R 复制（）函数不适用于 Rcpp 函数

Question

我在 R 中使用 replicate() 函数通过 Rcpp 函数生成随机数时遇到问题。考虑 R 中的以下函数：

trial <- function(){rnorm(1)}
replicate(10, trial())

它从高斯分布中生成 10 个随机数。它工作得很好，并产生如下结果：

 [1]  0.7609912 -0.2949613  1.8684363 -0.3358377 -1.6043926  0.2706250  0.5528813  1.0228125 -0.2419092 -1.4761937 

但是，我有一个 c++ 函数getRan()，它可以从高斯分布中生成一个随机数。我再次使用复制来调用这样的函数：

replicate(10,getRan())

它创建一个相同数字的向量，如下所示：

> replicate(10,getRan())
 [1] -1.374932 -1.374932 -1.374932 -1.374932 -1.374932 -1.374932 -1.374932 -1.374932 -1.374932 -1.374932
> replicate(10,getRan())
 [1] -0.3273785 -0.3273785 -0.3273785 -0.3273785 -0.3273785 -0.3273785 -0.3273785 -0.3273785 -0.3273785 -0.3273785
> replicate(10,getRan())
 [1] -0.7591953 -0.7591953 -0.7591953 -0.7591953 -0.7591953 -0.7591953 -0.7591953 -0.7591953 -0.7591953 -0.7591953
> replicate(10,getRan())
 [1] -1.698935 -1.698935 -1.698935 -1.698935 -1.698935 -1.698935 -1.698935 -1.698935 -1.698935 -1.698935

但是，如果我多次调用该函数，它工作正常：

 getRan()
[1] 1.345227
> getRan()
[1] 0.3555393
> getRan()
[1] 1.587241
> getRan()
[1] 0.5313518

那么这里的问题是什么？ replicate() 函数是否重复从 getRan() 返回的相同函数，而不是多次调用 getRan()？是bug吗？

PS：我知道我可以使用rnorm(n)来生成n个正常的随机数，但是，我想使用c++函数在生成随机数的基础上进行更复杂的计算

PPS：这是我的 C++ 代码：

double getRan(){
  unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
  std::default_random_engine generator(seed);
  std::normal_distribution<double> distribution (0.0,1.0);
  double epi = distribution(generator);
  return epi;
}

Answer 1

这是一个反例，表明它工作得很好：

代码

trialR <- function() { rnorm(1) }
Rcpp::cppFunction("double trialC() { return R::rnorm(0.0, 1.0); }")
Rcpp::cppFunction("Rcpp::NumericVector trialSugar() { return Rcpp::rnorm(1.0, 0.0, 1.0); }")

set.seed(123); replicate(3, trialR())
set.seed(123); replicate(3, trialC())
set.seed(123); replicate(3, trialSugar())

输出

通过Rscript 确保新的会话等pp：

edd@rob:/tmp$ Rscript so50543659.R 
[1] -0.560476 -0.230177  1.558708
[1] -0.560476 -0.230177  1.558708
[1] -0.560476 -0.230177  1.558708
edd@rob:/tmp$ 

Answer 2

德克斯的回答是正确的。您应该使用 R 的 RNG。如果你坚持在 C++ 中使用 RNG，你可以使用这样的东西：

#include <Rcpp.h>
// [[Rcpp::plugins(cpp11)]]
#include <random>

namespace {
  std::default_random_engine generator(std::random_device{}());
  std::normal_distribution<double> distribution (0.0,1.0);  
}

// [[Rcpp::export]]
double getRan(){
  return distribution(generator);
}

/*** R
replicate(10,getRan())
*/

这避免了在每次函数调用时创建std::default_random_engine（和std::normal_distribution）的新实例。这很重要，因为只能保证从一个 RNG 重复抽取 RNG 的属性。不适用于从不同 RNG 中重复抽取（希望是不同的）种子。

顺便说一句，在我的系统上，您的原始代码不会多次生成相同的数字。如果您在使用std::random_device 时遇到问题并且正在使用Windows，您可能会受到this mingw bug 的影响。在这种情况下，按时间播种是更好的选择。