这是一个快速比较mvtnorm::rmvnorm
and an Rcpp
给出的实施here作者:艾哈迈杜·迪科。所显示的时间是从维度范围从 500 到 2500 的多元正态分布中抽取 100 次的时间。从下图中,您可能可以推断出维度为 10000 所需的时间。时间包括生成随机数的开销mu
向量和diag
矩阵,但这些在不同方法中是一致的,并且对于所讨论的维度来说是微不足道的(例如 0.2 秒diag(10000)
).
library(Rcpp)
library(RcppArmadillo)
library(inline)
library(mvtnorm)
code <- '
using namespace Rcpp;
int n = as<int>(n_);
arma::vec mu = as<arma::vec>(mu_);
arma::mat sigma = as<arma::mat>(sigma_);
int ncols = sigma.n_cols;
arma::mat Y = arma::randn(n, ncols);
return wrap(arma::repmat(mu, 1, n).t() + Y * arma::chol(sigma));
'
rmvnorm.rcpp <-
cxxfunction(signature(n_="integer", mu_="numeric", sigma_="matrix"), code,
plugin="RcppArmadillo", verbose=TRUE)
rcpp.time <- sapply(seq(500, 5000, 500), function(x) {
system.time(rmvnorm.rcpp(100, rnorm(x), diag(x)))[3]
})
mvtnorm.time <- sapply(seq(500, 2500, 500), function(x) {
system.time(rmvnorm(100, rnorm(x), diag(x)))[3]
})
plot(seq(500, 5000, 500), rcpp.time, type='o', xlim=c(0, 5000),
ylim=c(0, max(mvtnorm.time)), xlab='dimension', ylab='time (s)')
points(seq(500, 2500, 500), mvtnorm.time, type='o', col=2)
legend('topleft', legend=c('rcpp', 'mvtnorm'), lty=1, col=1:2, bty='n')