从多元正态分布中有效地随机抽取答案

【问题标题】：Efficiently randomly drawing from a multivariate normal distribution从多元正态分布中有效地随机抽取
【发布时间】：2017-09-25 19:40:41
【问题描述】：

只是想知道是否有人遇到过他/她需要从非常高维的多元正态分布（比如维度 = 10,000）中随机抽取的问题，因为 mvtnorm 包的 rmvnorm 函数对于那个。

我知道这个article 有一个Rcpp 实现mvtnorm 包的dmvnorm 函数，所以我想知道rmvnorm 是否存在等效的东西？

【问题讨论】：

标签： r rcpp

【解决方案1】：

这里是 mvtnorm::rmvnorm 和 Rcpp 实现的快速比较，由 Ahmadou Dicko 给出的 here。给出的时间是从维度范围从 500 到 2500 的多元正态分布中抽取 100 次的时间。从下图中，您可能可以推断出维度为 10000 所需的时间。时间包括生成随机 mu 向量的开销和diag 矩阵，但这些在不同的方法中是一致的，并且对于所讨论的维度来说是微不足道的（例如，diag(10000) 为 0.2 秒）。

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')

【讨论】：

非常感谢 jbaums。我需要在此代码中添加什么才能使用 openmp（多核）？
我没有使用 openmp 的经验，但这里的瓶颈是 Cholesky 分解，所以this 可能会有所帮助。 This 也很有用。
jbaums 我最初得到的性能比你上面显示的还要好。但是，当我将采样的数量从 100 增加到 10000 时，Rcpp 大大减慢（而 mvtnorm 实际上没有改变）。我无法弄清楚发生了什么。内存管理？你有什么想法吗？