c++ - 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 源代码，最好的解决方案是用原生C++ 编写自定义的sample()。 Rcpp::sample() 的核心组件是 unif_rand()。在现代版本的 Fisher–Yates Shuffle 中集成 unif_rand()。问题解决了。

注意 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() is not recommended for serious random-number generation needs. It is recommended to use C++11's random number generation facilities to replace 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 的信息不足以为 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;
}