RcppArmadillo Gamma 分布在具有相同种子的平台之间不同

Question

我正在研究 a package ，它使用来自 RcppArmadillo 的随机数。该软件包运行 MCMC 算法，为了获得精确的再现性，用户应该能够设置随机数种子。执行此操作时，似乎用于从 Gamma 分布生成随机数的 arma::randg() 函数会跨平台返回不同的值。 arma::randu() 或 arma::randn() 不是这种情况。这可能与 arma::randg() 需要 C++11 这一事实有关吗？

这是我在运行 R3.5.2 的 macOS 10.13.6 上得到的:

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 package 创建于 2019-02-22 (v0.2.1)}

这是我在运行 R3.5.2 的 Windows 10 上得到的结果:

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 package 创建于 2019-02-22 (v0.2.1)}

可以看出，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 package 创建于 2019-02-22 (v0.2.1)}

但是，如警告和结果所示，种子不会以这种方式设置。

在使用 arma::randg() 时，是否有办法在平台之间获得可重现的结果，或者我是否需要使用 RcppArmadillo 中提供的一些其他随机数生成器来实现 Gamma 分布？

更新

正如评论中所指出的，使用 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 package 创建于 2019-02-22 (v0.2.1)}

这为我解决了这个问题。但是，我不确定问题是否已解决，因为这似乎是意外行为，因此请保持打开状态。

Answer 1

总结评论中的讨论:

• 对于 gamma 分布，Armadillo 使用 C++11 的 random header 中的 std::gamma_distribution，c.f. https://gitlab.com/conradsnicta/armadillo-code/blob/9.300.x/include/armadillo_bits/arma_rng_cxx11.hpp#L165
• 用 C++ 生成标准随机数分布的算法是 implementation-defined .
• 如果您需要跨平台可重现的结果，最简单的解决方案是使用通过 R::rgamma 或 Rcpp::rgamma 在 R 中实现的 Gamma 分布。