提升随机数生成器的丢弃行为

Question

我正在用一个适配器类包装 boost 随机数生成器来实现蒙特卡洛例程。在对类的成员函数编写单元测试时，我假设 .discard(unsigned int N) 的行为是抽取 N 个随机数而不存储它们，从而推进 rng 的状态。升压代码是：

void discard(boost::uintmax_t z)
{
    if(z > BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD) {
        discard_many(z);
    } else {
        for(boost::uintmax_t j = 0; j < z; ++j) {
            (*this)();
        }
    }
}

这支持了我的假设。但是，我发现 .discard(1) 产生的序列与没有丢弃的相同序列不是一个数字。代码：

#include <iostream>
#include <iomanip>
#include <random>
#include <boost/random.hpp>

int main()
{
    boost::mt19937 uGenOne(1);
    boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > distOne(uGenOne, boost::normal_distribution<>());

    boost::mt19937 uGenTwo(1);
    boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > distTwo(uGenTwo, boost::normal_distribution<>());
    distTwo.engine().discard(1);

    unsigned int M = 10;
    std::vector<double> variatesOne(M);
    std::vector<double> variatesTwo(M);

    for (unsigned int m = 0; m < M; ++m) {
        variatesOne[m] = distOne();
        variatesTwo[m] = distTwo();
    }

    for (unsigned int m = 0; m < M; ++m)
        std::cout << std::left << std::setw(15) << variatesOne[m] << variatesTwo[m] << std::endl;

    return 0;
}

输出

2.28493        0.538758  
-0.668627      -0.0017866
0.00680682     0.619191  
0.26211        0.26211   
-0.806832      -0.806832 
0.751338       0.751338  
1.50612        1.50612   
-0.0631903     -0.0631903
0.785654       0.785654  
-0.923125      -0.923125

我对 .discard 运作方式的解释是否不正确？为什么前三个输出中的两个序列不同，然后是相同的？

（此代码在 cygwin 上的 msvc 19.00.23918 和 g++ 4.9.2 上编译，结果相同）。

Answer 1

这里的问题似乎是引擎没有被正确修改，或者发行版正在添加一些额外的工作。如果我们直接用引擎就好了

int main()
{
    boost::mt19937 uGenOne(1);

    boost::mt19937 uGenTwo(1);
    uGenTwo.discard(1);

    unsigned int M = 10;
    std::vector<double> variatesOne(M);
    std::vector<double> variatesTwo(M);

    for (unsigned int m = 0; m < M; ++m) {
        variatesOne[m] = uGenOne();
        variatesTwo[m] = uGenTwo();
    }

    for (unsigned int m = 0; m < M; ++m)
        std::cout << std::left << std::setw(15) << variatesOne[m] << variatesTwo[m] << std::endl;

    return 0;
}

它产生

1.7911e+09     4.28288e+09
4.28288e+09    3.09377e+09
3.09377e+09    4.0053e+09
4.0053e+09     491263
491263         5.5029e+08
5.5029e+08     1.29851e+09
1.29851e+09    4.29085e+09
4.29085e+09    6.30312e+08
6.30312e+08    1.01399e+09
1.01399e+09    3.96591e+08

正如我们预期的那样，这是一个 1 移位序列，因为我们丢弃了第一个输出。

所以你对丢弃的工作方式是正确的。我不确定为什么通过boost::variate_generator 进行操作时会出现差异。我不明白为什么前三个数字不同，但所有其余的输出都匹配。

Answer 2

只是添加上一个答案的 cmets 中的重要细节。正如@NathanOliver 提到的， .discard 增加了生成器，它将制服发送到正态分布，将制服转换为正态分布。 boost::normal_distribution 使用Ziggurat algorithm，这是一种“接受/拒绝”类型的算法。它绘制一个随机制服，对其进行操作，然后检查它是否在所需的分布中。如果不是，它被拒绝并且一个新的随机制服。

for(;;) {
        std::pair<RealType, int> vals = generate_int_float_pair<RealType, 8>(eng);
        int i = vals.second;
        int sign = (i & 1) * 2 - 1;
        i = i >> 1;
        RealType x = vals.first * RealType(table_x[i]);
        if(x < table_x[i + 1]) return x * sign;
        if(i == 0) return generate_tail(eng) * sign;
        RealType y = RealType(table_y[i]) + uniform_01<RealType>()(eng) * RealType(table_y[i + 1] - table_y[i]);
        if (y < f(x)) return x * sign;
    }

关键是，如果最后一个if 失败，那么for 循环将再次启动，并再次触发对generate_int_float_pair 的调用。这意味着底层生成器增加的次数是未知的。

因此，正常序列将具有不同的数字，直到子序列中的拒绝总和相同，此时剩余的统一序列相同。这发生在问题中发布的示例中的第三个位置。 （它实际上有点微妙，因为在 Ziggurat 算法中可以调用底层生成器一次或两次，但本质是相同的——一旦序列同步，它们就永远不会产生不同的变量）。