Поведение сброса для форсированных генераторов случайных чисел
Я обертываю импульсные генераторы случайных чисел классом адаптера для реализации процедуры Монте-Карло. При написании модульных тестов для функций-членов класса я предполагал, что поведение .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 работает неправильно? Почему две последовательности отличаются в первых трех выходах, а затем идентичны?
(Этот код был скомпилирован в msvc 19.00.23918 и g ++ 4.9.2 в cygwin с одинаковыми результатами).
Решение
Проблема здесь заключается в том, что движок не изменяется должным образом или дистрибутив добавляет дополнительную работу. Если мы используем двигатель напрямую, как
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, Я не понимаю, почему первые три числа отличаются, но все остальные выходные данные совпадают.
Другие решения
Просто чтобы добавить важную деталь, которая была в комментариях к предыдущему ответу. Как упоминалось в @NathanOliver, .discard увеличивает генератор, который отправляет униформу в нормальное распределение, которое преобразует униформу в нормальное распределение. boost :: normal_distribution использует Алгоритм зиккурата который является алгоритмом типа «принятие / отклонение». Он рисует случайную форму, манипулирует ею, а затем проверяет, находится ли она в нужном распределении. Если нет, то отклоняется и новая случайная форма.
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 будет запущен снова. Это означает, что количество приращений базового генератора неизвестно.
Следовательно, нормальные последовательности будут иметь разные номера до тех пор, пока сумма, в которой сумма отклонений в подпоследовательности не будет одинаковой, в этот момент оставшаяся однородная последовательность будет идентичной. Это произошло на третьей позиции в примере, размещенном в вопросе. (На самом деле это немного более тонко, потому что базовый генератор может быть вызван один или два раза в алгоритме Зиккурата, но суть одна и та же — после того, как последовательности синхронизируются, они никогда не смогут производить разные переменные).