c – 从Rcpp改变R的种子以保证再现性

我试图在rcpp中写一个函数r(d,n).该函数从正态分布N(0,d)返回n个随机抽取.应该很好地定义这个函数,因此只要d和n不改变它们的值,函数就应该返回相同的draw.

如果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应该是双倍的.有什么想法吗？谢谢！

解决方法:

我认为正在发生的问题是你试图驱散犰狳提供的限制为标准正常的兰特,例如N(0,1),使其与N(0,d)匹配.有两种方法可以解决这个问题,因为它是标准法线.

选项1：使用统计属性

第一种方法涉及将样本乘以d的平方根,例如SQRT(d)*样品.这是可能的,因为方差的随机变量性质和期望给出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函数来生成随机值.唯一值得关注的是你不能使用arma :: arma_rng :: set_seed()来设置种子,因为Section 6.3 of Writing R Extensions中详细介绍了R / C交互的限制.如果你使用它,那么you would get warned with：

When called from R, the RNG seed has to be set at the R level via set.seed()

在第一次检测到的电话上.

有了这个说,这是一个简短的代码示例,其中我们乘以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是一个整数的限制,可以使用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

标签：c,r,rcpp
来源： https://codeday.me/bug/20190929/1831898.html