Rcpp: multidimensional integration via sourceCpp and cubature

cubature.cpp

// [[Rcpp::depends(cubature)]]
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <exp_cubature.h>


using namespace Rcpp;

const double log2pi = std::log(2.0 * M_PI);
arma::vec mu;
arma::mat Sigma;
arma::mat rooti;
double rootisum;

// [[Rcpp::export]]
double dmvn(arma::vec x) {
  int xdim = x.n_elem;
  arma::vec z = rooti * (x - mu) ;
  double out = -(static_cast<double>(xdim)/2.0) * log2pi - 0.5 * arma::sum(z%z) + rootisum;
  return(exp(out));
}

int dmvn_wrapper(unsigned ndim, const double *x, void *fdata, unsigned fdim, double *fval){
  fval[0] = dmvn(arma::vec(x, ndim));
  return 0;
}

// [[Rcpp::export]]
List cubature() {
  double vecmu[] = {0.0, 0.0, 0.0};
  double vecsigma[] = {1.0, 3.0/5.0, 1.0/3.0, 3.0/5.0, 1.0, 11.0/15.0, 1.0/3.0, 11/15.0, 1.0};
  mu = arma::vec(vecmu, 3);
  Sigma = arma::mat(vecsigma, 3, 3);
  rooti = arma::trans(arma::inv(trimatu(arma::chol(Sigma))));
  rootisum = arma::sum(log(rooti.diag()));

  int retcode;
  double xmin[] = {-0.5, -0.5, -0.5};
  double xmax[] = {1.0, 4.0, 2.0};
  double val;
  double error;

  retcode = hcubature(1, &dmvn_wrapper, NULL, 3, xmin, xmax, 0, 0, 1e-5, ERROR_INDIVIDUAL, &val, &error);
  return Rcpp::List::create(_["value"] = val, _["error"] = error, _["retcode"] = retcode);
}


/*** R
cubature()
***/