在 Rcpp 中不正确使用随机生成器进行 Gamma 分布绘制

Question

我在 Rcpp 中编写了以下代码，用于生成 Gamma 分布随机变量，但是每次运行我都采用相同的输出。我读到，为了有不同的实现，我必须使用这里解释的随机生成器

我使用的代码如下，我做错了吗？？基本上，我觉得我使用的代码与发布的问题完全相同，但出于某种原因，我的情况不起作用。

  #include <RcppArmadilloExtensions/sample.h>
  #include <random>
  #include <iostream>
  #include <math.h>
  #include<Rmath.h>
  using namespace Rcpp;


// [[Rcpp::depends(RcppArmadillo)]]


// [[Rcpp::export]]

double Gama_Draw_1(double a , double b){
  
  std::default_random_engine generator;
  std::gamma_distribution<double> d(a,b);
  
  return d(generator);
}


// [[Rcpp::export]]


double Gama_Draw_2(double a , double b){
  
  
  std::mt19937 prng{ std::random_device{}() }; 
  std::gamma_distribution<double> d(a,b); 
  return d(prng);
  
  
}

// [[Rcpp::export]]


arma::vec Gama_Draw_3(double a , double b){
  
  arma::vec S(10);
  
  std::random_device rd; 
  std::mt19937 gen(rd()); 
  std::gamma_distribution<> distrib(a, b);
  
  for(int i=0; i<10; ++i){
    S[i] = distrib(gen);
  }
  
  return S;
}

例如，

Gama_Draw_1(1,1)
[1] 0.5719583
Gama_Draw_2(1,1)
[1] 1.065146
Gama_Draw_3(1,1)
           [,1]
 [1,] 1.0651459
 [2,] 0.8683230
 [3,] 2.7298295
 [4,] 0.7930546
 [5,] 0.3722135
 [6,] 0.7317269
 [7,] 0.2569218
 [8,] 2.0916565
 [9,] 0.9033717
[10,] 1.4198487

无论我运行多少次，我总是得到相同的结果。

Answer 1

Rcpp 包连接来自 R 的随机数生成器（和特殊函数），同时正确连接 R 本身内部的（高质量）RNG。

所以我建议你依靠它，因为它实际上很容易使用：

> Rcpp::cppFunction("NumericVector mygamma(int n, double a, double b) { return Rcpp::rgamma(n, a, b); }")
> set.seed(123)
> mygamma(3, 0.5, 0.5)
[1] 0.05796143 1.14034608 0.00742452
> mygamma(3, 0.5, 0.5)
[1] 0.0753645 0.2474057 0.0151682
> set.seed(123)
> mygamma(3, 0.5, 0.5)
[1] 0.05796143 1.14034608 0.00742452
>

重置种子让我们获得相同的抽奖，不重置让我们获得不同的抽奖。正如它应该。还有：

> set.seed(123)
> rgamma(3, 0.5, 2.0)
[1] 0.05796143 1.14034608 0.00742452
>

请注意，Gamma 的第二个参数用作 R 版本和编译版本之间的倒数。

编辑： 为了完整起见，下面是一个使用 C++11 RNG 的工作版本。显然我们不能将 R.

的值加倍

代码

#include <random>
#include <Rcpp.h>

std::mt19937 engine(123);

// [[Rcpp::export]]
Rcpp::NumericVector mygamma(int n, double a, double b) {
  Rcpp::NumericVector v(n);
  std::gamma_distribution<> gamma(a, b);
  for (auto i=0; i<n; i++) {
    v[i] = gamma(engine);
  }
  return v;
}

/*** R
mygamma(3, 0.5, 0.5)
mygamma(3, 0.5, 0.5)
*/

输出

> Rcpp::sourceCpp("~/git/Whosebug/68235476/answer.cpp")

> mygamma(3, 0.5, 0.5)
[1] 0.168918 1.254678 0.311767

> mygamma(3, 0.5, 0.5)
[1] 0.0451722 0.5485451 0.0443048
>

在 Rcpp 中不正确使用随机生成器进行 Gamma 分布绘制

Not proper use of Random Generator for Gamma distribution draws in Rcpp

c++

r

rcpp

代码

输出