使用 RcppGSL 的狄利克雷分布

Question

我有一个目前用 R 编写的 Gibbs 采样器，我正在尝试使用包 Rcpp 和 RcppGSL 使其更快。现在给我带来问题的是，我似乎无法将随机变量生成器用于狄利克雷分布。这是一个不能在我的电脑上运行的简短脚本：

#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...

为了它的价值，我最近还在 Python 中实现了相同的 Gibbs 采样器，使用 GSL RNG（在同一台计算机上），一切似乎都有效。

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() 作为第三个参数插入到对 Rcpp 一无所知的 GSL 函数中。那不行。

重新仔细阅读 RcppGSL 文档和示例。这是可以修复的，但需要不同的方法。

编辑： 即使按照 GSL 标准，该界面也有些奇怪。但是通过用

替换导致错误的行，您确实得到了一个比较干净的解决方案

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

之后它会构建并运行——不过您仍然需要为 RNG 引擎播种。但这在 GSL 参考手册中有所介绍...

使用 RcppGSL 的狄利克雷分布

Dirichlet distribution with RcppGSL

c++

r

gsl

rcpp