r- rcpparmadillo伽马分布在具有相同种子的平台之间有所不同

Question

我正在处理一个软件包，该软件包使用了rcpparmadillo的随机数。该软件包运行MCMC算法，为了获得精确的可重复性，用户应该能够设置随机数种子。在执行此操作时，似乎arma::randg()的函数用于从伽马分布生成随机数的函数在平台上返回不同的值。arma::randu()或arma::randn()并非如此。它可能与arma::randg()需要C 11？

有关

这是我在MacOS 10.13.6上获得的，运行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

^{由Reprex软件包（v0.2.1）在2019-02-22创建}

这是我在Windows 10上获得的，运行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

^{由Reprex软件包（v0.2.1）在2019-02-22创建}

可以看出，使用arma::randg()生成的随机数在内部一致，但平台之间有所不同。

我试图使用Armadillo文档中的说明来设置种子：

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

^{由Reprex软件包（v0.2.1）在2019-02-22创建}

但是，正如警告所说的那样，结果表明，种子没有以这种方式设置。

使用arma::randg()时，是否有一种方法可以在平台之间获得可复制的结果，或者我需要使用rcpparmadillo中的其他一些随机数生成器来实现伽马分布？

update

如注释中指出的那样，使用R::rgamma()解决了此问题。以下代码在Mac和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

^{由Reprex软件包（v0.2.1）在2019-02-22创建}

这为我解决了。但是，我不确定问题是否已解决，因为这似乎是出乎意料的行为，所以将其打开。

Answer 1

在评论中总结讨论：

• 对于伽马发行版，Armadillo使用C 11的random标头std::gamma_distribution，C.F。https://gitlab.com/conradsnicta/armadillo-code/blob/9.300.x/include/armadillo_bits/arma_rng_cxxxx11.hpp#l165
• 在C 中生成标准随机数分布的算法是实现定义的。
• 如果您需要跨平台可重复的结果，最简单的解决方案是使用R在R中通过R::rgamma或Rcpp::rgamma实现的伽马分布。