Rcpp如何在Rcpp中生成随机多变量法向量？

Question

我想生成一些大的随机多元（超过6维）正常样本。在R中，许多软件包可以做到这一点，比如rmnorm，rmvn ......但问题是速度！所以我试图通过Rcpp编写一些C代码。我在网上浏览了一些教程，但似乎在多变量分布中没有“糖”，STL库中也没有。

任何帮助表示赞赏！

谢谢！

Answer 1

我不确定Rcpp会有帮助，除非你找到一个很好的算法来逼近你的多变量（cholesky，svd等）并使用Eigen（RccpEigen）或Armadillo（使用RcppArmadillo）对它进行编程。

下面是使用Cholesky分解和一种方法（RCPP）犰狳

#include <RcppArmadillo.h> 

// [[Rcpp::depends(RcppArmadillo)]] 

// [[Rcpp::export]] 

using namespace arma; 
using namespace Rcpp; 

mat mvrnormArma(int n, mat sigma) { 
    int ncols = sigma.n_cols; 
    mat Y = randn(n, ncols); 
    return Y * chol(sigma); 
} 

现在，在纯的R

mvrnormR <- function(n, sigma) { 
    ncols <- ncol(sigma) 
    matrix(rnorm(n * ncols), ncol = ncols) %*% chol(sigma) 
} 

一个天真的实现也可以检查是否万物工作

sigma <- matrix(c(1, 0.9, -0.3, 0.9, 1, -0.4, -0.3, -0.4, 1), ncol = 3) 
cor(mvrnormR(100, sigma)) 
cor(MASS::mvrnorm(100, mu = rep(0, 3), sigma)) 
cor(mvrnormArma(100, sigma)) 

现在我们以它为基准吧

require(bencharmk) 
benchmark(mvrnormR(1e4, sigma), 
      MASS::mvrnorm(1e4, mu = rep(0, 3), sigma), 
      mvrnormArma(1e4, sigma), 
      columns=c('test', 'replications', 'relative', 'elapsed')) 


## 2 MASS::mvrnorm(10000, mu = rep(0, 3), sigma)   100 
## 3     mvrnormArma(10000, sigma)   100 
## 1      mvrnormR(10000, sigma)   100 
## relative elapsed 
## 2 3.135 2.295 
## 3 1.000 0.732 
## 1 1.807 1.323 

在这个例子中，我使用单位方差和零均值的正态分布，但你可以很容易地推广到高斯分布定制均值和方差。

希望这会有所帮助