Verschiedene Ergebnisse der Monte-Carlo-Integration aufgrund der Ausgabe

Question

Ich bin gerade auf einen C++ Bug/Feature gestoßen, den ich nicht vollständig verstehen kann und hoffte, dass jemand mit besseren Kenntnissen von C++ hier in die richtige Richtung zeigen könnte .

Nachstehend finden Sie meinen Versuch, die Fläche unter der Gauß-Kurve mit Monte-Carlo-Integration zu ermitteln. Das Rezept ist:

1. Generieren Sie eine große Anzahl von normalverteilten Zufallsvariablen (mit Mittelwert von Null und Standardabweichung von eins).
2. Platz diese Zahlen.
3. Nehmen Sie den Durchschnitt aller Quadrate. Der Durchschnitt wird eine sehr genaue Schätzung der Fläche unter der Kurve sein (im Falle des Gaußschen ist es 1,0).

Der folgende Code besteht aus zwei einfachen Funktionen: rand_uni, die eine zufällige Variable gibt, gleichmäßig verteilt zwischen null und eins, und rand_norm, die eine (eher schlecht, aber ‚gut genug für Regierungsarbeit‘) ist Annäherung einer normalverteilten Zufallsvariablen.

main läuft eine milliardenfach durch eine Schleife und ruft jedes Mal rand_norm auf, quadriert es mit pow und addiert zur akkumulierenden Variable. Nach dieser Schleife wird das akkumulierte Ergebnis einfach durch die Anzahl der Durchläufe geteilt und als Result=<SOME NUMBER> an das Terminal gedruckt.

Das Problem liegt in sehr quirky Verhalten des Codes folgt: wenn jede erzeugte Zufallsvariable cout gedruckt wird (ja, eine Milliarde mal), dann ist das Endergebnis korrekt, unabhängig von der verwendeten Compiler (1,0015, was ziemlich ist nah, was ich will). Wenn ich die Zufallsvariable nicht in jeder Schleifeniteration drucke, bekomme ich inf unter gcc und 448314 unter clang.

Ehrlich gesagt, dies ist nur der Kopf und da es meine erste (-ish) Begegnung mit C++ ist, weiß ich nicht wirklich, was das Problem sein könnte: ist es etwas mit pow? Ist cout seltsam?

Jeder Hinweis würde sehr geschätzt werden!

Die Quelle

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

} 

Das 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

Die result in double rand_norm() nicht auf 0 initialisiert, bevor Sie + = Zufallswerte zu starten.

Alle Couts verwenden und setzen einen Teil des Stack-Speichers zurück, der später beim Aufruf von rand_norm verwendet wird, um Ihr Problem zu beheben, wenn Sie die Couts ausführen.

Ändern Sie die erste Zeile in double result = 0.0;, um es zu beheben.

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

Auch Sie Ihre Zufallszahlengenerator Samen sollten, werden Sie genau die gleiche Sequenz von ohne es Zufallszahlen Pseudo bekommen immer, fügen Sie ein

srand(time(NULL)); 

vor dem ersten Aufruf von rand(), oder schauen Sie sich einige bessere Zufallszahlengeneratoren an, z Boost.Random

Answer 2

In Ihrer rand_norm() - Funktion ist das Ergebnis nicht initialisiert, das ist Ihr Fehler. Nur für den Fall, dass Sie sich wundern, warum es beim Drucken der Werte funktionierte, weil eine unitialized Variable den Wert hat, den die Speicherregion hat, in Ihrem Code ist es der Stapel, also cout < < Ihre Stapelwerte geändert. Hier ist die feste Version Ihrer Funktion:

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