r-在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))));
}

编辑：我已经模拟了所有三种方法并比较了时间。

// [[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（矢量化）来自R/Rcpp，参数化为rate = 1/scale——我在前面的问题中向您解释了这一点？

编辑：这与上次的答案基本相同：与使用"1/over"参数化的R版本相比，对rgamma的C++版本的简单C++调用。

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>