与Rcpp一起使用Rcpp

Question

我希望在与Rcpp一起使用的代码中使用R函数-Rcpp-。我遇到了一个关于没有足够的参数的错误--我不知道这些参数应该是什么，因为它们与R中函数所使用的参数不匹配，我也没有幸运地使用"::Rf_foo“语法从Rcpp代码中访问R函数。

下面是我的代码的简化版本(是的，我正在编写一个gibbs采样器)。

#include <Rcpp.h>                                                                                                                                     
using namespace Rcpp;                                                                                                                                 

// C++ implementation of the R which() function.                                                                                                      
int whichC(NumericVector x, double val) {                                                                                                             
  int ind = -1;                                                                                                                                       
  int n = x.size();                                                                                                                                   
  for (int i = 0; i < n; ++i) {                                                                                                                       
    if (x[i] == val) {                                                                                                                                
      if (ind == -1) {                                                                                                                                
        ind = i;                                                                                                                                      
      } else {                                                                                                                                        
        throw std::invalid_argument( "value appears multiple times." );                                                                               
      }                                                                                                                                               
    } // end if                                                                                                                                       
  } // end for                                                                                                                                        
  if (ind != -1) {                                                                                                                                    
    return ind;                                                                                                                                       
  } else {                                                                                                                                            
    throw std::invalid_argument( "value doesn't appear here!" );                                                                                      
    return -1;                                                                                                                                        
  }                                                                                                                                                   
}                                                                                                                                                     

// [[Rcpp::export]]                                                                                                                                   
int multSample(double p1, double p2, double p3) {                                                                                                     
  NumericVector params(3);                                                                                                                            
  params(0) = p1;                                                                                                                                     
  params(1) = p2;                                                                                                                                     
  params(2) = p3;                                                                                                                                     

  // HERE'S THE PROBLEM.                                                                                                                              
  RObject sampled = rmultinom(1, 1, params);                                                                                                          
  int out = whichC(as<NumericVector>(sampled), 1);                                                                                                    
  return out;                                                                                                                                         
}

我是c++新手，所以我意识到这段代码的很多地方可能是nooby和低效的。对于如何改进我的c++代码，我是开放的，但是我的首要任务是了解rmultinom业务。谢谢!

顺便说一句，我很抱歉和这条线有相似之处，但是

1. 这个答案不符合我的目的
2. 这种差异可能足以引发另一个问题(你认为是这样吗？)
3. 这个问题一年前发表并得到了答复。

Answer 1

下面是user95215修改后的答案，以便进行编译，另一个版本则更多地采用Rcpp风格：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerVector oneMultinomC(NumericVector probs) {
    int k = probs.size();
    SEXP ans;
    PROTECT(ans = Rf_allocVector(INTSXP, k));
    probs = Rf_coerceVector(probs, REALSXP);
    rmultinom(1, REAL(probs), k, &INTEGER(ans)[0]);
    UNPROTECT(1);
    return(ans);
}

// [[Rcpp::export]]
IntegerVector oneMultinomCalt(NumericVector probs) {
    int k = probs.size();
    IntegerVector ans(k);
    rmultinom(1, probs.begin(), k, ans.begin());
    return(ans);
}

Answer 2

如果我试图编译您的代码，就会得到编译器错误：

> Rcpp::sourceCpp('~/scratch/multSample.cpp')
multSample.cpp:33:21: error: no matching function for call to 'rmultinom'
  RObject sampled = rmultinom(1, 1, params);
                    ^~~~~~~~~
/Library/Frameworks/R.framework/Resources/include/Rmath.h:449:6: note: candidate function not viable: requires 4 arguments, but 3 were provided
void    rmultinom(int, double*, int, int*);
        ^
1 error generated.

正如它所暗示的，您没有正确地指定参数。注意，与其他函数相比，rmultinom接口有点尴尬:它填充了*rn所指向的内存，而不是返回一个新对象(有它自己新分配的内存)。

如果您查看R源，您将看到接口，您也可以看到使用它的示例这里 (实际上，stats提供了一个包装器函数，它可以进行更多的参数检查等等)。但是注意它是如何在这里使用的：

rmultinom(size, REAL(prob), k, &INTEGER(ans)[ik]);

换句话说，它通过将该内存的地址传递给INTSXP函数来填充一个名为ans的rmultinom。

因此，如果您想使用Rcpp中的这个函数，您将不得不做一些类似的事情--但也许这应该得到类似的糖矢量化处理，以避免接口的丑陋。

你可以试着做这样的事情：

IntegerMatrix sampled(nrow, ncol);
rmultinom(1, 1, params, sampled.begin());

或者其他类似的东西。

Answer 3

凯文提供的例子和链接使我找到了一个有效的答案。有一些类型的争论。我编写了一个函数，允许您从多项式分布中抽取一个向量。代码在下面。

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector oneMultinomC(NumericVector probs) {
    int k = probs.size();
    SEXP ans;
    PROTECT(ans = RF_allocVector(INTSXP, k));
    probs = RF_coerceVector(probs, REALSXP);
    rmultinom(1, REAL(probs), k, &INTEGER(ans)[0]);
    UNPROTECT(1);
    return ans;
}

我不明白这里发生了什么。特别是，我不明白第四种观点。我知道这是一个指向存储输出的内存位置的指针，但我不理解“”位。尽管如此，它还是有效的。吉布斯样本走了，男孩和女孩。