-4
我想生成一些大的随机多元(超过6维)正常样本。在R中,许多软件包可以做到这一点,比如rmnorm,rmvn ......但问题是速度!所以我试图通过Rcpp编写一些C代码。我在网上浏览了一些教程,但似乎在多变量分布中没有“糖”,STL库中也没有。Rcpp如何在Rcpp中生成随机多变量法向量?
任何帮助表示赞赏!
谢谢!
我想生成一些大的随机多元(超过6维)正常样本。在R中,许多软件包可以做到这一点,比如rmnorm,rmvn ......但问题是速度!所以我试图通过Rcpp编写一些C代码。我在网上浏览了一些教程,但似乎在多变量分布中没有“糖”,STL库中也没有。Rcpp如何在Rcpp中生成随机多变量法向量?
任何帮助表示赞赏!
谢谢!
我不确定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
在这个例子中,我使用单位方差和零均值的正态分布,但你可以很容易地推广到高斯分布定制均值和方差。
希望这会有所帮助
不错的答案+1。你想为Rcpp Gallery(gallery.rcpp.org)写这个吗? – 2013-03-10 19:54:38
@DirkEddelbuettel非常感谢..很高兴能为Rcpp画廊作出贡献。我将分叉github回购并提交一个pull请求。 – dickoa 2013-03-10 21:08:25
简单的文件作为.Rmd或.cpp与标记也很好,但你当然也去正式和长途... – 2013-03-11 00:30:42