Реализация mvtnorm :: pmvnorm в rcpp медленнее, чем исходная функция R

Вместо того, чтобы пытаться использовать для этого дополнительную библиотеку, я бы попытался использовать C API, экспортированный mvtnorm, c.f. https://github.com/cran/mvtnorm/blob/master/inst/NEWS#L44-L48. При этом я обнаружил три причины, по которым результаты различаются. Один из них также отвечает за разницу в производительности:

1. mvtnorm использует RNG, хотя он был удален из библиотеки, которую вы используете, c.f. https://github.com/zhanxw/libMvtnorm/blob/master/libMvtnorm/randomF77.c.

2. Ваша triangl функция неверна. Он возвращает нижнюю треугольную матрицу в порядке столбцов. Однако базовый код fortran ожидает этого в строковом порядке, ср. https://github.com/cran/mvtnorm/blob/master/src/mvt.f#L36-L39 и https://github.com/zhanxw/libMvtnorm/blob/master/libMvtnorm/mvtnorm.cpp#L60

3. libMvtnorm использует 1e-6 вместо 1e-3 в качестве относительной точности, ср. https://github.com/zhanxw/libMvtnorm/blob/master/libMvtnorm/mvtnorm.cpp#L65. Это также отвечает за разницу в производительности.

Мы можем проверить это, используя следующий код:

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::depends(mvtnorm)]]
#include <mvtnormAPI.h>

//[[Rcpp::export]]
arma::vec triangl(const arma::mat& X){
  int n = X.n_cols;
  arma::vec res(n * (n-1) / 2);
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < i; ++j) {
      res(j + i * (i-1) / 2) = X(i, j);
    }
  }
  return res;
}

// [[Rcpp::export]]
double pmvnorm_cpp(arma::vec& bound,
           arma::vec& lowertrivec,
           double abseps = 1e-3){

  int n = bound.n_elem;
  int nu = 0;
  int maxpts = 25000;     // default in mvtnorm: 25000
  double releps = 0;      // default in mvtnorm: 0
  int rnd = 1;            // Get/PutRNGstate

  double* bound_ = bound.memptr();
  double* correlationMatrix = lowertrivec.memptr();
  double* lower = new double[n];
  int* infin = new int[n];
  double* delta = new double[n];

  for (int i = 0; i < n; ++i) {
    infin[i] = 0; // (-inf, bound]
    lower[i] = 0.0;
    delta[i] = 0.0;
  }

  // return values
  double error;
  double value;
  int inform;

  mvtnorm_C_mvtdst(&n, &nu, lower, bound_,
           infin, correlationMatrix, delta,
           &maxpts, &abseps, &releps,
           &error, &value, &inform, &rnd);
  delete[] (lower);
  delete[] (infin);
  delete[] (delta);

  return value;
}

/*** R
set.seed(1)
covar <- rWishart(1, 10, diag(5))[,,1]
sds <- diag(covar) ^-.5
corrmat <- diag(sds) %*% covar %*% diag(sds)
triang <- triangl(corrmat)
bounds <- c(0.5, 0.9, 1, 4, -1)
set.seed(1)
system.time(cat(mvtnorm::pmvnorm(upper=bounds, corr = corrmat), "\n"))
set.seed(1)
system.time(cat(pmvnorm_cpp(bounds, triang, 1e-6), "\n"))
set.seed(1)
system.time(cat(pmvnorm_cpp(bounds, triang, 0.001), "\n"))
 */

Полученные результаты:

> system.time(cat(mvtnorm::pmvnorm(upper=bounds, corr = corrmat), "\n"))
0.04896221 
   user  system elapsed 
  0.000   0.003   0.003 

> system.time(cat(pmvnorm_cpp(bounds, triang, 1e-6), "\n"))
0.04895756 
   user  system elapsed 
  0.035   0.000   0.035 

> system.time(cat(pmvnorm_cpp(bounds, triang, 0.001), "\n"))
0.04896221 
   user  system elapsed 
  0.004   0.000   0.004 

При одинаковом ГСЧ (и состоянии ГСЧ), правильной нижней треугольной корреляционной матрице и той же относительной точности результаты идентичны, а производительность сопоставима. При более высокой точности страдает производительность.

Все это для отдельного файла, использующего Rcpp::sourceCpp. Чтобы использовать это в пакете, вам нужно добавить LinkingTo: mvtnorm в ваш DESCRIPTION файл.