将 R 的种子从 Rcpp 更改为保证可重复性

Question

我正在尝试在 rcpp 中编写函数 r(d, n)。该函数从正态分布 N(0, d) 返回 n 次随机抽取。这个函数应该被很好地定义，因此只要 d 和 n 不改变它们的值，该函数应该返回相同的绘制。

如果 d 被限制为整数，这将不是问题，在这种情况下，我可以设置种子并完成工作

// set seed
// [[Rcpp::export]]
void set_seed(unsigned int seed) {
  Rcpp::Environment base_env("package:base");
  Rcpp::Function set_seed_r = base_env["set.seed"];
  set_seed_r(seed);  
}

// function r(d, n)
// [[Rcpp::export]]
vec randdraw(int d, int n){
  set_seed(d);
  vec out = randn(n);
  return out;
}

但显然我不想将 d 限制为整数。理想情况下，d 应该是两倍。有什么想法吗？谢谢！

Answer 1

我认为正在发生的问题是您试图驱散 Armadillo 提供的 randn，该randn 仅限于标准正常，例如N(0,1)，使其匹配 N(0,d)。有两种方法可以解决这个问题，因为这是一个标准的常态。

选项 1：使用统计属性

第一种方法只是将样本乘以d 的平方根，例如sqrt(d)*sample。这是可能的，因为方差和期望的随机变量属性给出了 sqrt(d)*N(0, 1) ~ N(0, sqrt(d)^2) ~ N(0, d)。

这里需要注意的更重要的事情之一是set_seed() 函数将在RcppArmadillo hooks into R's RNG library 的Armadillo 配置访问::Rf_runif 函数以生成随机值后工作。唯一需要关注的是，由于Section 6.3 of Writing R Extensions 中详述的 R/C++ 交互的限制，您不能使用arma::arma_rng::set_seed() 设置种子。如果你确实使用这个，那么you would get warned with：

当从 R 调用时，RNG 种子必须通过 set.seed() 在 R 级别设置

在第一个检测到的呼叫上。

话虽如此，下面是一个简短的代码示例，我们将其乘以 sqrt(d)。

代码：

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// set seed
// [[Rcpp::export]]
void set_seed(double seed) {
    Rcpp::Environment base_env("package:base");
    Rcpp::Function set_seed_r = base_env["set.seed"];
    set_seed_r(std::floor(std::fabs(seed)));
}

// function r(d, n)
// [[Rcpp::export]]
arma::vec randdraw(double d, int n){
    set_seed(d);              // Set a seed for R's RNG library
    // Call Armadillo's RNG procedure that references R's RNG capabilities
    // and change dispersion slightly.
    arma::vec out = std::sqrt(std::fabs(d))*arma::randn(n);
    return out;
}

输出：

> randdraw(3.5, 5L)
           [,1]
[1,] -0.8671559
[2,] -1.9507540
[3,]  2.9025090
[4,] -1.2953745
[5,]  2.0799176

注意：由于rnorm 过程与arma::randn 生成不同，因此没有直接等效项。

选项 2：依赖 R 的 RNG 函数

第二个，也是明显更好的解决方案，是明确地依赖 R 的 RNG 函数。之前，由于 RcppArmadillo 的配置，我们隐式使用了 R 的 RNG 库。我倾向于使用这种方法，因为您已经假设在使用 the set_seed() function 时代码是特定于 R 的（免责声明：我写了这篇文章）。如果您担心d 是integer 的限制，则可以使用std::floor(std::fabs(seed)) 将double 稍微强制转换为int。使用 Rcpp::r*() 或 R::r*() 生成值后，将使用 an advanced ctor 创建一个重用现有内存分配的犰狳向量。

代码：

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// set seed
// [[Rcpp::export]]
void set_seed(double seed) {
    Rcpp::Environment base_env("package:base");
    Rcpp::Function set_seed_r = base_env["set.seed"];
    set_seed_r(std::floor(std::fabs(seed)));
}

// function r(d, n)
// [[Rcpp::export]]
arma::vec randdraw(double d, int n){
    set_seed(d);                                      // Set a seed for R's RNG library
    Rcpp::NumericVector draws = Rcpp::rnorm(n, 0.0, d); // Hook into R's Library
    // Use Armadillo's advanced CTOR to re-use memory and cast as an armadillo object.
    arma::vec out = arma::vec(draws.begin(), n, false, true);
    return out;
}

输出：

> randdraw(3.21,10)
             [,1]
 [1,] -3.08780627
 [2,] -0.93900757
 [3,]  0.83071017
 [4,] -3.69834335
 [5,]  0.62846287
 [6,]  0.09669786
 [7,]  0.27419092
 [8,]  3.58431878
 [9,] -3.91253230
[10,]  4.06825360
> set.seed(3)
> rnorm(10, 0, 3.21)
 [1] -3.08780627 -0.93900757  0.83071017 -3.69834335  0.62846287  0.09669786  0.27419092  3.58431878 -3.91253230  4.06825360