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 

在這個例子中，我使用單位方差和零均值的正態分佈，但你可以很容易地推廣到高斯分佈定製均值和方差。

希望這會有所幫助