Rcpp如何在Rcpp中生成随机多元法线向量? [英] Rcpp How to generate random multivariate normal vector in Rcpp?

问题描述

我想生成一些较大的随机多元变量(超过6个维度)普通样本.在R中，许多软件包都可以执行此操作，例如rmnorm，rmvn ...但是问题在于速度！因此，我尝试通过Rcpp编写一些C代码.我在线上看过一些教程，但似乎STL库中都没有用于多变量分发的糖".

I would like to generate some large random multivariate (more than 6 dimensions) normal samples. In R, many packages can do this such as rmnorm, rmvn... But the problem is the speed! So I tried to write some C code through Rcpp. I went through some tutorial online but it seems there is no "sugar" for multivariate distribution, neither in STL library.

感谢您的帮助！

谢谢！

推荐答案

我不确定Rcpp是否会有所帮助，除非您找到一个好的算法来近似您的多元变量(cholesky，svd等)并使用Eigen对其进行编程( RccpEigen)或犰狳(使用RcppArmadillo).

I'm not sure that Rcpp will help unless you find a good algorithm to approximate your multivariate (cholesky, svd, etc.) and program it using Eigen (RccpEigen) or Armadillo (using RcppArmadillo).

这是使用Cholesky分解和(Rcpp)Armadillo的一种方法

Here is one approach using the Cholesky decomposition and (Rcpp)Armadillo

#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语言实现天真的实现

Now a naive implementation in pure R

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

您还可以检查一切是否正常

You can also check if everythings work

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

现在让我们对其进行基准测试

Now let's benchmark it

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

在此示例中，我使用具有单位方差和零均值的正态分布，但是您可以轻松地将其推广到具有自定义均值和方差的高斯分布.

In this example I used a normal distribution with unit variance and null mean but you could easily generalize to gaussian distribution with custom mean and variance.

希望这会有所帮助

这篇关于Rcpp如何在Rcpp中生成随机多元法线向量?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！