Rcpp Как создать случайный многомерный нормальный вектор в Rcpp?

Question

Я хотел бы создать несколько случайных многомерных (более 6 измерений) обычных образцов. В R многие пакеты могут делать это, например, rmnorm, rmvn ... Но проблема в скорости! Поэтому я попытался написать код C через Rcpp. Я просмотрел несколько обучающих онлайн, но, похоже, нет никакого «сахара» для многомерного распространения, ни в библиотеке STL.

Любая помощь приветствуется!

Спасибо!

Answer 1

Я не уверен, что Rcpp поможет, если вы не найдете хороший алгоритм для аппроксимации своего многомерного (cholesky, svd и т. Д.) И запрограммируйте его с помощью Eigen (RccpEigen) или Armadillo (используя RcppArmadillo).

Вот один подход, использующий разложение Холецкого и (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

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

Вы также можете проверить, если Everythings работать

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 

В этом примере я использовал нормальное распределение с единичной дисперсией и нулевым средним значением, но вы можете легко обобщить на гауссовский дистрибутив с пользовательским средним значением и дисперсией.

Надеюсь, что это поможет