0

ブースト乱数ジェネレーターをアダプター クラスでラップして、モンテカルロ ルーチンを実装しています。クラスのメンバー関数で単体テストを作成するとき、.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) の結果のシーケンスは、破棄のない同じシーケンスと 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 の動作についての私の解釈は間違っていますか? 2 つのシーケンスの最初の 3 つの出力が異なり、その後は同一になるのはなぜですか?

(このコードは、msvc 19.00.23918 および cygwin の g++ 4.9.2 でコンパイルされ、同じ結果が得られました)。

4

2 に答える 2

0

前の回答のコメントにあった重要な詳細を追加するだけです。@NathanOliver が述べたように、.discard はジェネレーターをインクリメントし、ユニフォームを正規分布に送信し、ユニフォームを正規分布に変換します。boost::normal_distributionは、「承認/拒否」タイプのアルゴリズムであるZiggurat アルゴリズムを使用します。ランダムなユニフォームを描画し、それを操作してから、それが目的の分布にあるかどうかを確認します。そうでない場合、それは拒否され、新しいランダムなユニフォームになります。 

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が再びトリガーされることです。これは、基礎となるジェネレーターがインクリメントされる回数がわからないことを意味します。

したがって、通常のシーケンスは、サブシーケンスの拒否の合計が同じになるポイントまで、異なる数を持ち、その時点で残りの均一なシーケンスは同一になります。これは、質問に投稿された例の 3 番目の位置で発生しました。(実際には、基になるジェネレーターをジッグラト アルゴリズムで 1 回または 2 回呼び出すことができるため、もう少し微妙ですが、本質は同じです。シーケンスが同期すると、異なる変量を生成することはできません)。

于 2016-08-17T19:24:30.667 に答える