从 Rcpp 中获取与基 R 相同的整数样本

Posted 2023-02-16

【中文标题】从 Rcpp 中获取与基 R 相同的整数样本

【英文标题】：Get the same sample of integers from Rcpp as base R

【发布时间】：2020-02-07 19:06:28

【问题描述】：

是否有可能从Rcpp 获得与从基本 R 的 sample 获得相同的整数 sample？

我曾尝试使用Rcpp::sample 和Rcpp::RcppArmadillo::sample，但它们返回的值不同——示例代码如下。此外，https://gallery.rcpp.org/articles/using-the-Rcpp-based-sample-implementation/ 帖子的 快速示例 部分返回了来自 Rcpp 和基础 R 的相同示例，但是，我无法重现这些结果（我在最后附上了这段代码）。

可以这样做吗/请问我做错了什么？

我的尝试：

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

// [[Rcpp::export]]
Rcpp::IntegerVector mysamp1( int n) 
  Rcpp::IntegerVector v = Rcpp::sample(n, n);
  return v;


// [[Rcpp::export]]
Rcpp::IntegerVector mysamp2(int n)   
  Rcpp::IntegerVector i = Rcpp::seq(1,n);
  Rcpp::IntegerVector v = wrap(Rcpp::RcppArmadillo::sample(i,n,false));
  return v;


// set seed https://***.com/questions/43221681/changing-rs-seed-from-rcpp-to-guarantee-reproducibility
// [[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)));


// [[Rcpp::export]]
Rcpp::IntegerVector mysamp3( int n, int seed) 
  set_seed(seed); 
  Rcpp::IntegerVector v = Rcpp::sample(n, n);
  return v;



/***R
set.seed(1)
sample(10)
#  [1]  9  4  7  1  2  5  3 10  6  8
set.seed(1)
mysamp1(10)
#  [1]  3  4  5  7  2  8  9  6 10  1
set.seed(1)
mysamp2(10)
#  [1]  3  4  5  7  2  8  9  6 10  1
mysamp3(10, 1)
#  [1]  3  4  5  7  2  8  9  6 10  1

*/

来自Using the RcppArmadillo-based Implementation of R's sample() 库帖子的代码，在我的系统上返回FALSE：

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadilloExtensions/sample.h>
using namespace Rcpp ;

// [[Rcpp::export]]
CharacterVector csample_char( CharacterVector x, 
                              int size,
                              bool replace, 
                              NumericVector prob = NumericVector::create()) 
  CharacterVector ret = RcppArmadillo::sample(x, size, replace, prob) ;
  return ret ;


/*** R
N <- 10
set.seed(7)
sample.r <- sample(letters, N, replace=T)

set.seed(7)
sample.c <- csample_char(letters, N, replace=T)

print(identical(sample.r, sample.c))
# [1] FALSE
*/

【问题讨论】：

或许设置RNGKind

本指南有帮助吗：gallery.rcpp.org/articles/random-number-generation

感谢 akrun 做到了。例如，设置RNGversion("1.6.2") 或RNGkind("Super") 会产生相同的结果。

不错的快速回答 :) 我认为我们有一个未解决的问题票可以将示例实现更新为 R 现在所做的。如果有人有时间贡献...

我忘记了细节，但它是 R 的 RNG 的一个有据可查的变化，因为有人注意到 IIRC 使用抽样（N 非常大）存在偏差。这就是为什么你要在 R 中打开一个选项来获得旧的（匹配的）行为。现在有时间的人应该将 Rcpp 端更新为 R 中的 new 行为，这样您就不需要将 R 转回来。

【参考方案1】：

将 cmets 编译为答案。 Akrun 指出，通过设置RNGkind 或RNGversion，我们可以复制结果。来自 DirkEddelbuettel； “R 的 RNG 发生了变化，因为有人注意到 IIRC 使用采样（N 非常大）存在偏差。这就是为什么你要在 R 中打开一个选项来获得旧的（匹配）行为. " RalfStubner 指出这是一个已知问题：https://github.com/RcppCore/RcppArmadillo/issues/250 和 https://github.com/RcppCore/Rcpp/issues/945

目前 R 使用不同的默认采样器，导致不同的结果

RNGkind(sample.kind = "Rejection")
set.seed(1)
sample(10)
# [1]  9  4  7  1  2  5  3 10  6  8
set.seed(1)
mysamp1(10)
# [1]  3  4  5  7  2  8  9  6 10  1

但是，可以使用更早的版本来使用

RNGkind(sample.kind = "Rounding")
#Warning message:
#  In RNGkind("Mersenne-Twister", "Inversion", "Rounding") : non-uniform 'Rounding' sampler used

set.seed(1)
sample(10)
# [1]  3  4  5  7  2  8  9  6 10  1
set.seed(1)
mysamp1(10)
# [1]  3  4  5  7  2  8  9  6 10  1