Генерировать случайную гамма в Rcpp

Question

У меня есть параметры формы и масштаба вектора для генерации чисел из случайного распределения гаммы. Я не мог найти способ в R;

lambda<-matrix(rgamma(p,ak,scale=1/bk),p,1) 

Так что я пытаюсь создать в цикле и, кажется, я теряю время, которое является основной целью использования Rcpp.

for(int i = 0; i<p;i++){ 
lambda(i)=arma::conv_to<double>::from(arma::randg<arma::mat (1,1,arma::distr_param(ak(i),pow(bk(i),-1)))); 
} 

Редактировать: Я моделировал все 3 метода и сравнивал время.

// [[Rcpp::plugins(cpp11)]] 
// [[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 
// [[Rcpp::export]] 

arma::mat rg(const arma::mat& w,double scale,int num) { 
int p = w.n_rows; 
arma::mat lambda(num,p,arma::fill::zeros); 


for (int j=0;j<num;j++){ 
for(int i = 0; i<p;i++){ 
    lambda(j,i)=arma::conv_to<double>::from(arma::randg<arma::mat>(1,1,arma::distr_param(w(i),scale))); 
} 
} 
return mean(lambda.rows(0,num-1),0);  
} 

Это коды R;

fun<-function(w,sca,num){ 
p=nrow(w) 
lambda=matrix(0,num,p) 
for (jj in 1:num){ 
    for (i in 1:p){ 
     lambda[jj,i]=rgamma(1,w[i],scale=sca) 
    } 
} 
return(colMeans(lambda)) 
} 

fun2<-function(w,sca,num){ 
p=nrow(w) 
lambda=matrix(0,num,p) 
for (jj in 1:num){  
     lambda[jj,]=rgamma(p,w,scale=sca)  
} 
return(colMeans(lambda)) 
} 

Вот результат;

a=matrix(c(1,2,3)) 


    microbenchmark(rg(a,2,100000),fun(a,2,1e5),fun2(a,2,1e5)) 
    Unit: milliseconds 
       expr  min  lq  mean median  uq  max neval 
    rg(a, 2, 1e+05) 891.6197 900.0092 952.3833 915.6609 991.7918 1271.9117 100 
    fun(a, 2, 1e+05) 1018.6645 1059.3994 1171.9160 1130.7481 1200.3192 1754.3445 100 
    fun2(a, 2, 1e+05) 321.8309 339.2670 373.5317 355.5914 395.1515 604.4365 100 

Answer 1

Вы понимаете, что arma::randg от пространства имен Armadillo?

В то время как R::rgamma (скалярной) или Rcpp::rgamma (vectorised) взяты из R/Rcpp с параметрированием rate = 1/scale - что я объяснил вам в this earlier question of yours?

Edit: Это в основном тот же ответ as the last time: простой C++ позвонить в версии C++ в rgamma по сравнению с R версии с помощью «1/над» параметризацию.

R> set.seed(42); stem(rgamma(100, 0.5, 2.0)) # calling from R 

    The decimal point is 1 digit(s) to the left of the | 

    0 | 00000000000001111111112222222333345566677777778888990111111122345666 
    2 | 002257268 
    4 | 235788 
    6 | 269457 
    8 | 70345 
    10 | 28 
    12 | 1 
    14 | 
    16 | 0 

R> cppFunction("NumericVector callrgamma(int n, double shape, double scale) { 
+          return(rgamma(n, shape, scale)); }") 
R> set.seed(42); stem(callrgamma(100, 0.5, 1.0/2.0)) # calling from C++ 

    The decimal point is 1 digit(s) to the left of the | 

    0 | 00000000000001111111112222222333345566677777778888990111111122345666 
    2 | 002257268 
    4 | 235788 
    6 | 269457 
    8 | 70345 
    10 | 28 
    12 | 1 
    14 | 
    16 | 0 

R>