Using the Stan Math C++ Library

General
1

I am trying to use the Stan function lkj_corr_rng inside of my c++ code. I can execute the c++ code available in the section “Using Higher-Order Functions in the StanHeaders Package” of the documentation in Using the Stan Math C++ Library. But I am failing when I try to adapt that code to use the lkj_corr_rng function. For example, when I execute the code bellow

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(StanHeaders)]]
#include <stan/math/prim/mat/prob/lkj_corr_rng.hpp>
#include <stan/math.hpp>
#include <Rcpp.h>
#include <RcppEigen.h>

// [[Rcpp::plugins(cpp14)]]

// [[Rcpp::export]]
Eigen::MatrixXd rlkgcorr_cpp(int K, double eta) {
  return stan::math::lkj_corr_rng(K, eta);
}

I have the following error: no matching function for call to ‘lkj_corr_rng(int&, double&)’. What am I doing wrong?

2

That error occurs because lkj_corr_rng needs to access a random number generator. I didn’t test this but the following should compile:

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(StanHeaders)]]
#include <boost/random/additive_combine.hpp>
#include <stan/math/prim/mat/prob/lkj_corr_rng.hpp>
#include <stan/math.hpp>
#include <Rcpp.h>
#include <RcppEigen.h>

// [[Rcpp::plugins(cpp14)]]

// [[Rcpp::export]]
Eigen::MatrixXd rlkgcorr_cpp(int K, double eta) {
  unsigned int seed = 1906563677;
  boost::ecuyer1988 rng = boost::ecuyer1988(seed);
  return stan::math::lkj_corr_rng(K, eta, rng);
}

This constructs a new RNG for each function call from a fixed seed so it produces the same pseudorandom matrix every time. For more random results you need to either pass a different seed every time or store the RNG between calls.

1 Like
3

Hi Niko, your code works for me! I just had to adapt it to

// [[Rcpp::export]]
Eigen::MatrixXd rlkjcorr_cpp(int D, double eta, unsigned int seed) {
  boost::ecuyer1988 rng = boost::ecuyer1988(seed);
  return stan::math::lkj_corr_rng(D, eta, rng);
}

and then use the function in the the c++ this way:

const int max_int = std::numeric_limits<int>::max();
unsigned int seed = Rcpp::sample(max_int,1)[0];
Eigen::MatrixXd R = rlkjcorr_cpp(D, eta, seed);

Thanks a lot!
Ricardo Pedroso