STL random_shuffle生成高度相关的序列

Question

请注意接受的答案指出问题在于重新种植.重播不是原因.没有重新种植的测试在发布前产生了很高的相关性.见注1.

我在R中生成了1,000,000个统一随机数,对序列进行了排序,并调用std::random_shuffle()此序列的副本100次.100个置换序列非常相关.但是,如果我不首先对统一数字进行排序,那么100个置换序列或多或少是不相关的.以下是代码.

// [[Rcpp::export]]
IntegerVector testRandomShuffle(IntegerVector x, int rd) // rd is the seed
{
  IntegerVector y(x.begin(), x.end()); // copy
  std::srand(rd); // seeding
  std::random_shuffle(y.begin(), y.end());
  return y;
}


/***R
v = runif(1000000)
vSorted = sort(v)
sqc = 1L : length(v) # indexes
rd = sample.int(length(v), 100) # random seeds


# Compute correlation matrices
corMatForUnsorted = cor(as.data.frame(lapply(rd, function(x) 
  v[testRandomShuffle(sqc, x)])))
corMatForSorted = cor(as.data.frame(lapply(rd, function(x) 
  vSorted[testRandomShuffle(sqc, x)])))


# plot histograms
par(mfrow = c(1, 2)) 
hist(corMatForUnsorted[abs(corMatForUnsorted) < 1], breaks = 200, xlab = 
  "Correlation for unsorted")
hist(corMatForSorted[abs(corMatForSorted) < 1], breaks = 200, xlab = 
  "Correlation for sorted")
*/

我做错了什么吗？我只是希望改组排序和未排序的序列产生或多或少相同的相关分布.这些相关性有多小是另一个故事.R的本机sample.int()排列函数的相同实验在两种情况下产生低相关性.

谢谢!

注1:问题是我在Windows上使用g ++ 4.9.3附带的Rtools 3.4.此版本的C++库中的shuffle函数工作不正常.

注2:确认Rcpp::sample()在多线程中有效.一个小测试用例:

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


struct testSampleInPara: public RcppParallel::Worker
{
  IntegerVector tmp;
  List rst;


  void operator() (std::size_t st, std::size_t end)
  {
    if(st == 0)
    {
      // is tmp / rst a copy or a reference ?
      std::cout << std::to_string((std::size_t)&tmp[0]) + "\n";
      IntegerVector rst0 = Rcpp::sample(tmp, 5);
      rst[0] = rst0; // assume rst not a copy
    }
    else // if(st == 1)
    {
      std::cout << std::to_string((std::size_t)&tmp[0]) + "\n";
      IntegerVector rst1 = Rcpp::sample(tmp, 10);
      rst[1] = rst1;
    }
  }


  testSampleInPara(IntegerVector tmp, List rst):
    tmp(tmp), rst(rst)
  {
    RcppParallel::parallelFor(0, 2, *this);
  }
};


// [[Rcpp::export]]
List testIfSampleCopy(IntegerVector tmp)
{
  List rst(2);
  testSampleInPara(tmp, rst);
  return rst;
}

/***R
testIfSampleCopy(1L : 10L)
# printout:
# 356036792
# 356036792
# [[1]]
# [1] 10  5  9  7  8
# 
# [[2]]
#  [1] 10  3  7  6  2  1  8  4  9  5
*/

我对Rcpp容器的体验对于多线程的性能很差.我通常创建指向Rcpp容器的起始元素的指针或指针数组,在线程之间共享这些指针和容器的大小.注意Rcpp::sample()接受并返回Rcpp容器.

注3:通过阅读Rcpp源代码,最好的解决方案是使用sample()本机C++ 编程自定义.中心的组成部分Rcpp::sample()是unif_rand().融入unif_rand()Fisher-Yates Shuffle的现代版本.问题解决了.

注意4:unif_rand()在多线程环境中使用会大大降低线程的速度.我没有时间阅读Dirk Eddelbuettel建议的文档,但我猜R的源同步unif_rand()对我们来说是不可见的,就像malloc()在C.最终的解决方案是包含// [[Rcpp::plugins("cpp11")]]和使用std::random.

Answer 1

std::random_shuffle(begin, end)经常使用std::rand,已知是一个坏的随机数发生器.从cppreference:

rand()不建议用于严重的随机数生成需求.建议使用C++ 11的随机数生成工具进行替换rand().

请std::shuffle改用.

// Note the lack of `int rd`. `std::random_device` is better for
// seeding purposes, but it is non-deterministic.
IntegerVector testShuffle(IntegerVector x)
{
  IntegerVector y(x.begin(), x.end()); // copy

  // std::mt19937 is a rather heavy type. As such, it's often recommended
  // to make it a static variable. If you will be calling this function
  // from multiple threads, you'd want to make it `thread_local` instead
  // of `static` (or otherwise avoid the data race on `engine`).
  static std::mt19937 engine = [] {
    // Using the Immediately Invoked Lambda Expression (IILE) idiom to
    // initialize the static variable.

    // Seed the RNG.
    std::random_device rd;

    // Note that there are better ways to seed the mersenne twister.
    // This way is flawed, as it can't possibly initialize all of the
    // mersenne twister's state, but it's the simplest way for
    // demonstration purposes
    std::mt19937 engine(rd());

    return engine;
  }();

  // You should be able to just use y.begin(), y.end()
  std::shuffle(y.begin(), y.end(), engine);
  return y;
}

如果你想要一个确定性的种子,请注意单个int信息不足以完全播种a std::mt19937,但你仍然可以使用它:

IntegerVector testShuffle(IntegerVector x, int seed)
{
  IntegerVector y(x.begin(), x.end());

  static std::mt19937 engine;

  // Not thread-friendly, but simple.
  // Also, note that you'll get bad results if you seed a mersenne twister
  // (or a lot of RNGs) with 0, so avoid that
  engine.seed(seed);

  std::shuffle(y.begin(), y.end(), engine);
  return y;
}