logo标签列表

RcppArmadillo的gamma分布在相同种子的不同平台上存在差异

rrcpparmadillorcpp11
3
我正在处理一个软件包，它使用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

本文创建于2019年02月22日，使用reprex包（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 (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

本文创建于2019年2月22日，使用 reprex package (v0.2.1)。

然而，正如警告所述，并且结果显示的那样，这种方式并不能设置种子。

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

更新

正如评论中指出的那样，使用 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

本文创建于2019年2月22日，使用reprex package (v0.2.1)。

这对我解决了问题。但是，我不确定问题是否已经解决，因为这似乎是意外的行为，所以保持开放状态。

4
相比于使用Rcpp的rgamma函数，为什么要选择使用Armadillo？如果在使用Rcpp::rgamma()等效函数时问题没有持续发生，您可以在计算中使用该函数的输出，或者如果类型不同，则在计算过程中将其转换为所需的类型。 - Oliver
2
其实，我使用Armadillo的原因只是因为整个程序包大量使用了其他Armadillo函数。现在我已经更改了它，现在它可以工作了。 - Øystein S
3
我还没有检查过，但可能 Armadillo 使用了 C++11 的 <random> 头文件中的 std::gamma_distribution 。不幸的是，这些分布是"实现定义"的，这可能是 WRE 不建议使用该头文件的原因。 - Ralf Stubner
3
查看源代码后，你是完全正确的。如果使用C++11，则 <random> 库已经包含在内，而 arma::randg() 仅支持 C++ 11。 - Øystein S
1
依赖于特定的随机数序列通常是一个不好的想法。如果R开发人员决定稍微改变RNG，比如出于速度原因怎么办？实现可重复性的更健壮的方法是通过随机数序列的统计学来实现。例如，取平均值并查看其是否落在特定范围内。 - hbrerkere
当然，但是如果我想提交一篇基于该软件包分析的论文，期刊可能需要一个脚本来精确地重现每个十进制位的结果。虽然结论不应取决于所选的种子，但类似聚类标签之类的东西可能会发生变化。例如，在两次给定的运行之间，聚类1和聚类2可能会交换。此外，这是为了软件包开发，并且让用户能够设置种子似乎是合理的。 - Øystein S
1个回答
2
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接