我正在尝试让 pmvnorm 的 Rcpp 版本至少与 R 中的 mvtnorm::pmvnorm 一样快。
我已经发现https://github.com/zhanxw/libMvtnorm并创建了一个带有相关源文件的Rcpp框架包。我添加了以下使用犰狳的函数(因为我在我编写的其他代码中使用它)。
//[[Rcpp::export]]
arma::vec triangl(const arma::mat& X){
arma::mat LL = arma::trimatl(X, -1); // omit the main diagonal
return LL.elem(arma::find(LL != 0));
}
//[[Rcpp::export]]
double pmvnorm_cpp(arma::vec& bound, arma::vec& lowtrivec){
double error;
int n = bound.n_elem;
double* boundptr = bound.memptr();
double* lowtrivecptr = lowtrivec.memptr();
double result = pmvnorm_P(n, boundptr, lowtrivecptr, &error);
return result;
}
构建包后,从 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)
rbenchmark::benchmark(pmvnorm_cpp(bounds, triang),
mvtnorm::pmvnorm(upper=bounds, corr = corrmat),
replications=1000)
这表明 pmvnorm_cpp 比 mvtnorm::pmvnorm 慢得多。结果不同。
> pmvnorm_cpp(bounds, triang)
[1] 0.04300643
> mvtnorm::pmvnorm(upper=bounds, corr = corrmat)
[1] 0.04895361
这让我很困惑,因为我认为基本的 Fortran 代码是相同的。我的代码中是否有某些东西使一切变得缓慢?或者我应该尝试直接移植 mvtnorm::pmvnorm 代码?我确实没有使用 Fortran 的经验。
感谢建议,否则请原谅我的无能。
编辑:为了与替代方案进行快速比较,如下:
//[[Rcpp::export]]
NumericVector pmvnorm_cpp(NumericVector bound, NumericMatrix cormat){
Environment stats("package:mvtnorm");
Function f = stats["pmvnorm"];
NumericVector lower(bound.length(), R_NegInf);
NumericVector mean(bound.length());
NumericVector res = f(lower, bound, mean, cormat);
return res;
}
与 R 调用具有本质上相同的性能(以下在 40 维 mvnormal 上):
> rbenchmark::benchmark(pmvnorm_cpp(bounds, corrmat),
+ mvtnorm::pmvnorm(upper=bounds, corr = corrmat),
+ replications=100)
test replications elapsed relative user.self sys.self
2 mvtnorm::pmvnorm(upper = bounds, corr = corrmat) 100 16.86 1.032 16.60 0.00
1 pmvnorm_cpp(bounds, corrmat) 100 16.34 1.000 16.26 0.01
所以在我看来,前面的代码中一定发生了一些事情。要么是我如何处理犰狳的事情,要么是其他事情是如何联系起来的。我认为与最后的实现相比应该有性能增益。