带RcppGSL的Dirichlet分布

Question

我有一个吉布斯采样器目前在R，我试图使它更快使用包Rcpp和RcppGSL。现在给我带来问题的是，我似乎无法使用随机变量生成器来处理dirichlet分布。这里有一个在我的电脑上不起作用的短脚本：

#include <Rcpp.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>

#include <RcppGSL.h>

using namespace Rcpp;

// [[Rcpp::depends(RcppGSL)]]
// [[Rcpp::export]]

NumericVector rdirichlet_cpp(NumericVector alpha) {
  int n = alpha.size();
  NumericVector results(n);

  // Allocate random number generator
  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);

  gsl_ran_dirichlet(r, n, alpha, results);

  // Release random number generator
  gsl_rng_free(r);

  return(results);
}

当我尝试使用sourceCpp()来获取它时，我会收到一条错误消息，表示存在no matching function for call to 'gsl_ran_dirichlet'。我几乎没有使用C/C++的经验，所以我可能犯了一个愚蠢的错误(我仍然不太确定Rcpp如何管理内存)。但是也许问题实际上是在RcppGSL包中，不知怎么的，它与一个不包含Dirichlet随机变量生成器的较旧版本的GSL相关联。

值得注意的是，我最近还使用GSL (在同一台计算机上)用Python实现了相同的Gibbs采样器，而且一切似乎都正常。

Answer 1

您报告的错误与我看到的错误不同：

gsldiri.cpp: In function ‘Rcpp::NumericVector rdirichlet_cpp(Rcpp::NumericVector)’:
gsldiri.cpp:20:41: error: cannot convert ‘Rcpp::NumericVector {aka Rcpp::Vector<14, Rcpp::PreserveStorage>}’ to ‘const double*’ for argument ‘3’ to ‘void gsl_ran_dirichlet(const gsl_rng*, size_t, const double*, double*)’
   gsl_ran_dirichlet(r, n, alpha, results);
                                         ^
make: *** [gsldiri.o] Error 1

这个非常有意义:您将Rcpp::NumericVector()作为第三个参数随机插入到一个GSL函数中，该函数对Rcpp一无所知。那是行不通的。

更仔细地阅读RcppGSL文档和示例。这是可以修复的，但它需要一种不同的方法。

编辑:即使按照GSL标准，界面也有些奇怪。但是，通过将引起错误的行替换为

gsl_ran_dirichlet(r, n, alpha.begin(), results.begin());

在此之后，它构建并运行-您仍然需要在RNG引擎上播种。但在GSL参考手册中.