将 rmultinom 与 Rcpp 一起使用

sin*_*wav 3 c++ r bayesian mcmc rcpp

我想在 C++ 代码中使用 R 函数rmultinom与 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;                                                                                                                                         
}
Run Code Online (Sandbox Code Playgroud)

我是 C++ 新手,所以我意识到很多代码可能都是新手且效率低下。我愿意接受有关如何改进我的 C++ 代码的建议,但我的首要任务是了解 rmultinom 业务。谢谢!

顺便说一句,对于与此线程的相似之处,我深表歉意,但是

  1. 答案对我的目的不起作用
  2. 这种差异可能足以引发不同的问题(你认为是这样吗?)
  3. 该问题一年前发布并回答。

Dir*_*tel 5

下面是 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);
}
Run Code Online (Sandbox Code Playgroud)