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

Question

I wrote the following code in Rcpp, for generating Gamma Distribution random variables, however each time that I run it I take the same output. I read that in order to have different realizations I have to use a random generator as explained here Why is this random generator always output the same number

The code that I use is the following, am I doing something wrong?? Basically, I feel that I use the exact same code as in the question posted, but for some reason in my case doesn't work.

  #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;
}

For example,

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

No matter how many times I run it I always get the same result.

I can't reproduce your result. If I create a reproducible example, with all includes and export statements then different values are generated. Can you edit your question to show an Rcpp example that compiles in case something else is going on. https://stackoverflow.com/questions/54822702/rcpparmadillo-gamma-distribution-differs-between-platforms-with-same-seed?answertab=votes#tab-top may be worth reading over `std::gamma_distribution` — user2957945, Jul 03 '21 at 11:34
@user2957945 I think that it can be reproduced, I included all the code in the .cpp file — Jonathan1234, Jul 03 '21 at 11:52
your gamma function code is different to what you initially showed -- the first version worked. — user2957945, Jul 03 '21 at 12:02
I think it is the same the only difference I think is that I forgot to write the Gamma with double m :/ — Jonathan1234, Jul 03 '21 at 12:03
no its not the same, your code was `std::mt19937 prng{ std::random_device{}() }; std::gamma_distribution d(a,b); return d(prng);` — user2957945, Jul 03 '21 at 12:04
Either way you must _pass a seed to the RNG_ or for example pass the `random_device` to your RNG constructor: https://en.cppreference.com/w/cpp/numeric/random/uniform_int_distribution But you work with Rcpp so you probably should use the R RNGs anyway... — Dirk Eddelbuettel, Jul 03 '21 at 12:06
@DirkEddelbuettel I think I did that on the ```Gama_Draw_3``` but it still gives the same output every time. I also, checked if the set. seed in R is the problem, so I said to Null and still the same. — Jonathan1234, Jul 03 '21 at 12:15
You _must_ keep the state of the RNG persistent and outside the call of the function. Else you _recreate it each time_ getting you the same values. — Dirk Eddelbuettel, Jul 03 '21 at 12:44

Dirk Eddelbuettel · Accepted Answer · 2021-07-03T12:46:32.217

The Rcpp package interfaces the random number generators (and special functions) from R, while properly interfacing the (high-quality) RNG inside R itself.

So I recommend you rely on that as it is in fact easy to use:

> 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
>

Resetting the seed gets us the same draw, not resetting gets us different draws. As it should. Also:

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

Note that the second parameter for the Gamma is used as a reciprocal value between the R version and the compiled version.

Edit: For completeness, a working version with C++11 RNG follows. Obviously we cannot double the values from R.

Code

#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)
*/

Output

> Rcpp::sourceCpp("~/git/stackoverflow/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
>

Thank you for reply!!And the question now is, is the ```Rcpp::rgamma``` faster than a built-in function in c++?? That's why I was seeking to make it work in c++ — Jonathan1234, Jul 03 '21 at 12:40
Sometimes I feel really stupid... So, the whole problem was exactly what you said I was using the seed inside the function -.- — Jonathan1234, Jul 03 '21 at 12:51

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

1 Answers1

Code

Output