Search code examples
rrcpparmadillorcpp11

RcppArmadillo gamma distribution differs between platforms with same seed

I am working on a package, which uses random numbers from RcppArmadillo. The package runs MCMC algorithms, and to obtain exact reproducibility, the user should be able to set a random number seed. When doing this, it seems like the arma::randg() function for generating random numbers from the gamma distribution returns different values across platforms. This is not the case for arma::randu() or arma::randn(). Could it be related to the fact that arma::randg() requires C++11?

Here is what I get on macOS 10.13.6, running R3.5.2:

library(Rcpp)
library(RcppArmadillo)

sourceCpp(code = '
#include <RcppArmadillo.h>
using namespace Rcpp;

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

 // [[Rcpp::export]]
 double random_gamma() {
 return arma::randg();
 }

 // [[Rcpp::export]]
 double random_uniform() {
 return arma::randu();
 }

 // [[Rcpp::export]]
 double random_normal() {
 return arma::randn();
 }
 '
)

replicate(2, {set.seed(1); random_gamma()})
#> [1] 1.507675 1.507675
replicate(2, {set.seed(432); random_gamma()})
#> [1] 0.02234341 0.02234341
replicate(2, {set.seed(1); random_uniform()})
#> [1] 0.2655087 0.2655087
replicate(2, {set.seed(1); random_normal()})
#> [1] -1.390378 -1.390378

Created on 2019-02-22 by the reprex package (v0.2.1)

Here is what I get on Windows 10, running R3.5.2:

library(Rcpp)
library(RcppArmadillo)

sourceCpp(code = '
          #include <RcppArmadillo.h>
          using namespace Rcpp;

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

          // [[Rcpp::export]]
          double random_gamma() {
          return arma::randg();
          }

          // [[Rcpp::export]]
          double random_uniform() {
          return arma::randu();
          }

          // [[Rcpp::export]]
          double random_normal() {
          return arma::randn();
          }
          '
)

replicate(2, {set.seed(1); random_gamma()})
#> [1] 0.2549381 0.2549381
replicate(2, {set.seed(432); random_gamma()})
#> [1] 0.2648896 0.2648896
replicate(2, {set.seed(1); random_uniform()})
#> [1] 0.2655087 0.2655087
replicate(2, {set.seed(1); random_normal()})
#> [1] -1.390378 -1.390378

Created on 2019-02-22 by the reprex package (v0.2.1)

As can be seen, the random numbers generated with arma::randg() are internally consistent, but differ between the platforms.

I tried to set the seed using the instructions in the Armadillo documentation:

library(Rcpp)
library(RcppArmadillo)

sourceCpp(code = '
          #include <RcppArmadillo.h>
          using namespace Rcpp;

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

          // [[Rcpp::export]]
          double random_gamma(int seed) {
          arma::arma_rng::set_seed(seed);
          return arma::randg();
          }
          '
)

replicate(4, random_gamma(1))
#> Warning in random_gamma(1): When called from R, the RNG seed has to be set
#> at the R level via set.seed()
#> [1] 1.3659195 0.6447221 1.1771862 0.9099034

Created on 2019-02-22 by the reprex package (v0.2.1)

However, as the warning tells, and the result shows, the seed does not get set this way.

Is there a way to get reproducible results between platforms when using arma::randg(), or do I need to implement a gamma distribution using some of the other random number generators available in RcppArmadillo?

Update

As pointed out in the comment, using R::rgamma() solves this problem. The following code returns the same numbers on both Mac and Windows:

library(Rcpp)

sourceCpp(code = '
          #include <Rcpp.h>
          using namespace Rcpp;

          // [[Rcpp::export]]
          double random_gamma() {
            return R::rgamma(1.0, 1.0);
          }
          '
)

replicate(2, {set.seed(1); random_gamma()})
#> [1] 0.1551414 0.1551414

Created on 2019-02-22 by the reprex package (v0.2.1)

This solves it for me. However, I am not sure if the issue is solved, since this seems like unexpected behavior, so leaving it open.

Solution

  • Summarizing the discussion in the comments: