Różne wyniki integracji Monte Carlo z powodu wyjścia

Question

Niedawno natknąłem się na błąd/funkcję C++, której nie mogę w pełni zrozumieć i miałem nadzieję, że ktoś o lepszej znajomości C++ może wskazać mi właściwy kierunek .

Poniżej znajduje się moja próba znalezienia obszaru pod krzywą Gaussa przy użyciu integracji Monte Carlo. Przepis to:

1. Generowanie dużej liczby normalnie rozłożonych zmiennych losowych (ze średnią zerową i odchyleniem standardowym jednego).
2. Kwadrat te liczby.
3. Weź średnią wszystkich kwadratów. Średnia będzie bardzo bliskim oszacowaniem powierzchni pod krzywą (w przypadku Gaussa jest to 1,0).

Poniższy kod składa się z dwóch prostych funkcji: rand_uni, która zwraca zmienną losową, równomiernie rozłożone między zero a jeden, a rand_norm, który jest (raczej słaba, ale „wystarczająco dobre dla pracy rządu”) przybliżeniem normalnie rozłożonej zmiennej losowej.

main przebiega przez pętlę miliard razy, wywołując rand_norm za każdym razem, wyrównując ją do pow i dodając do zmiennej akumulującej. Po tej pętli zakumulowany wynik jest po prostu dzielony przez liczbę przebiegów i drukowany na terminalu jako Result=<SOME NUMBER>.

Problem polega na bardzo dziwnym zachowaniu się kodu poniżej: kiedy każda wygenerowana zmienna losowa zostanie wydrukowana na cout (tak, miliard razy), wynik końcowy jest poprawny, niezależnie od użytego kompilatora (1.0015, co jest dość blisko, do tego, czego chcę). Jeśli nie wydrukuję zmiennej losowej w każdej iteracji pętli, otrzymam inf pod gcc i 448314 pod clang.

Szczerze mówiąc, to po prostu przeszkadza umysłowi i ponieważ jest to moje pierwsze (-jcze) spotkanie z C++, tak naprawdę nie wiem, jaki może być problem: czy jest to coś z pow? Czy cout zachowuje się dziwnie?

Każda wskazówka byłaby bardzo cenna!

Źródło

// 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; 

} 

Makefile

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 

Answer 1

result w double rand_norm() nie zostaje zainicjowany na 0 przed rozpoczęciem + = losowych wartości do niego.

Wszystkie użycie cout i resetowanie części pamięci stosu, która jest później używana przez wywołanie rand_norm, "naprawianie" twojego problemu, gdy robisz cout.

Zmień pierwszą linię na double result = 0.0;, aby to naprawić.

double rand_norm() { 
    double result = 0.0; 
    for(int i=12; i > 0; --i) { 
    result += rand_uni(); 
    } 
    return result - 6; 
}; 

Należy również podkreślić, nasion swój generator liczb losowych, zawsze możesz liczyć na dokładnie taką samą sekwencję pseudo losowych liczb bez niego, dodać

srand(time(NULL)); 

przed pierwszym wywołaniu rand() lub spojrzeć na lepsze generatory liczb losowych, np Boost.Random

Answer 2

W wyniku funkcji rand_norm() wynik nie jest inicjowany, to jest błąd. Na wypadek, gdybyś zastanawiał się, dlaczego zadziałał, gdy wydrukowałeś wartości, ponieważ zmienna zunifikowana będzie miała wartość, którą posiada region pamięci, w twoim kodzie jest stos, więc wywołanie cout < < zmieniło twoje wartości stosu. Oto poprawiona wersja twojej funkcji:

double rand_norm() { 
    double result = 0.0; 
    for(int i=12; i > 0; --i) { 
    result += rand_uni(); 
    } 
    return result - 6; 
};