c - Rcpp Rcppでランダムな多変量法線ベクトルを生成するには?

Question

大きなランダム多変量 (6 次元以上) の正規サンプルをいくつか生成したいと思います。R では、rmnorm や rmvn など、多くのパッケージでこれを実行できますが、問題は速度です。そこで、Rcpp を使って C コードを書いてみました。オンラインでいくつかのチュートリアルを実行しましたが、多変量分布の「砂糖」はSTLライブラリにもないようです。

どんな助けでも大歓迎です！

ありがとう！

Answer 1

多変量 (コレスキー、svd など) を近似するための優れたアルゴリズムを見つけて、Eigen (RccpEigen) または Armadillo (RcppArmadillo を使用) を使用してプログラムしない限り、Rcpp が役立つかどうかはわかりません。

コレスキー分解と (Rcpp)Armadillo を使用したアプローチの 1 つを次に示します。

#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

この例では、単位分散とゼロ平均の正規分布を使用しましたが、カスタムの平均と分散を使用してガウス分布に簡単に一般化できます。

お役に立てれば