Dirichlet distribution with RcppGSL

I have a Gibbs sampler currently written in R, and I'm trying to make it faster using the package Rcpp and RcppGSL. What is causing me problems right now is that I can't seem to be able to use the random variate generator for the dirichlet distribution. Here's a short script that doesn't work on my computer:

#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);
}

When I try to source this using sourceCpp(), I get an error message, saying that there is no matching function for call to 'gsl_ran_dirichlet'. I have close to no experience with C/C++, so I may be making a stupid mistake (I'm still not quite sure how Rcpp manages memory). But perhaps the issue is actually with the RcppGSL package, that somehow it is linked to an older version of GSL which does not contain the Dirichlet random variate generator...

For what it's worth, I also recently implemented that same Gibbs sampler in Python, using the GSL RNG (on the same computer), and everything seems to work.

Solution

You report a different error from the one I see:

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

and this one makes perfect sense: You somewhat randomly inserted the Rcpp::NumericVector() as a third argument to a GSL function which knows nothing about Rcpp. That cannot work.

Re-read the RcppGSL documentation and examples more carefully. This is fixable, but it needs a different approach.

Edit: Even by GSL standards that interface is somewhat odd. But you do get a somewhat clean solution by replacing your error-causing line with

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

after which it builds and runs---you still need to work on seeding the RNG engine though. But that is covered in the GSL Reference manual...