Гамма-распределение 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

^{Создано 2019-02-22 пакетом представлением (v0.2.1)}

Вот что я получаю на 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

^{Создано 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

^{Создано 2019-02-22 пакетом представлением (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-02-22 пакетом представлением (v0.2.1)}

Это решает это для меня. Однако я не уверен, что проблема решена, поскольку это выглядит как неожиданное поведение, поэтому оставляю его открытым.

Answer 1

Подводя итоги обсуждения в комментариях:

• Для гамма-распределения Armadillo использует std::gamma_distribution из C++11 random заголовок, см. https://gitlab.com/conradsnicta/armadillo-code/blob/9.300.x/include/armadillo_bits/arma_rng_cxx11.hpp#L165
• Алгоритмы создания стандартных распределений случайных чисел в C++ определяются реализацией.
• Если вам нужны кроссплатформенные воспроизводимые результаты, самое простое решение - использовать гамма-распределение, реализованное в R через R::rgamma или же Rcpp::rgamma,