Различные результаты интеграции Монте-Карло в связи с выходом
Я только недавно наткнулся на ошибку / функцию C ++, которую я не могу полностью понять, и надеялся, что кто-то с лучшим знанием C ++ здесь может указать мне правильное направление.
Ниже вы найдете мою попытку найти область под кривой Гаусса с использованием интеграции Монте-Карло. Рецепт:
- Генерация большого количества нормально распределенных случайных величин (со средним значением нуля и стандартным отклонением, равным единице).
- Квадрат эти числа.
- Возьмите среднее значение всех квадратов. Среднее будет очень близко оценить площадь под кривой (в случае Гаусса, это 1,0).
Код ниже состоит из двух простых функций: rand_uni, который возвращает случайную величину, равномерно распределенную между нулем и единицей, и rand_norm, что является (довольно плохим, но «достаточно хорошим для работы в правительстве») приближением нормально распределенной случайной величины.
main проходит через цикл миллиард раз, вызывая rand_norm каждый раз, возводя его в квадрат pow и добавление к накопительной переменной. После этого цикла накопленный результат просто делится на количество прогонов и выводится на терминал как Result=<SOME NUMBER>,
Проблема заключается в очень причудливом поведении кода ниже: когда каждая сгенерированная случайная величина печатается в cout (да, один миллиард раз), тогда конечный результат будет правильным независимо от используемого компилятора (1.0015, что довольно близко к тому, что я хочу). Если я не печатаю случайную переменную в каждой итерации цикла, я получаю inf под gcc и 448314 под clang,
Честно говоря, это просто поражает воображение, и, поскольку это мое первое (-ое) знакомство с C ++, я не совсем понимаю, в чем может быть проблема: pow? Есть ли cout вести себя странно?
Любая подсказка будет высоко ценится!
// Monte Carlo integration of the Gaussian curve
#include <iostream>
#include <cstdlib>
#include <cmath>using namespace std;enum {
no_of_runs = 1000000
};// uniform random variable
double rand_uni() {
return ((double) rand() / (RAND_MAX));
};// approximation of a normaly distributed random variable
double rand_norm() {
double result;
for(int i=12; i > 0; --i) {
result += rand_uni();
}
return result - 6;
};int main(const int argc,
const char** argv) {
double result = 0;
double x;
for (long i=no_of_runs; i > 0; --i) {
x = pow(rand_norm(), 2);
#ifdef DO_THE_WEIRD_THING
cout << x << endl; // MAGIC?!
#endif
result += x;
}
// Prints the end result
cout << "Result="<< result / no_of_runs
<< endl << endl;
}
CLANG=clang++
GCC=g++
OUT=normal_mc
default: *.cpp
$(CLANG) -o $(OUT).clang.a *.cpp
$(CLANG) -o $(OUT).clang.b -DDO_THE_WEIRD_THING *.cpp
$(GCC) -o $(OUT).gcc.a *.cpp
$(GCC) -o $(OUT).gcc.b -DDO_THE_WEIRD_THING *.cpp
Решение
result в double rand_norm() не инициализируется в 0, прежде чем начать + = случайные значения для него.
Все cout используют и сбрасывают некоторую часть стековой памяти, которая позже используется вызовом rand_norm, «исправляя» вашу проблему, когда вы выполняете cout.
Измените первую строку на double result = 0.0; починить это.
double rand_norm() {
double result = 0.0;
for(int i=12; i > 0; --i) {
result += rand_uni();
}
return result - 6;
};
Кроме того, вы должны заполнить свой генератор случайных чисел, вы всегда получите точно такую же последовательность псевдослучайных чисел без него, добавьте
srand(time(NULL));
перед вашим первым вызовом rand (), или посмотрите на некоторые лучшие генераторы случайных чисел, например Boost.Random
Другие решения
В вашей функции rand_norm () результат не инициализирован, это ваша ошибка. На тот случай, если вам интересно, почему это работает, когда вы печатаете значения, это потому, что унитарная переменная будет иметь значение, которое у нее есть в области памяти, в вашем коде это стек, поэтому вызывается cout<< изменил ваши значения стека. Вот исправленная версия вашей функции:
double rand_norm() {
double result = 0.0;
for(int i=12; i > 0; --i) {
result += rand_uni();
}
return result - 6;
};