Box-Müller-Methode in C++

[SeppJ]

tobias kempkens schrieb:

@SeppJ:
sry, da habe ich dich dann falsch verstanden.
Denn der letzte Satz war
"Tritt der Fehler immer noch bei der selben Pfadtiefe auf?"
und deine antwort war "NEIN". also nahm ich an, dein "nein" bezog sich auf den letzten satz und nicht auf seinen vorschlag.

Ahh, ok. Tschuldige dass das missverständlich war. Ich wollte Exhumed darauf hinweisen, dass dies ganz offensichtlich kein Genauigkeitsproblem ist, aber auch nicht seitenlang über Numerik schreiben um ihm zu erklären, warum.

edit: Falls möglich, wäre auch die rand-Implementierung deiner Standardbibliothek hilfreich. Bei meiner (libstdc++ auf 64-Bit Debian) tritt der Fehler nämlich nicht auf, selbst bei vielen Millionen Pfaden. Da spielt dann eventuell die Wahrscheinlichkeit der 0 eine Rolle, weil ich vielleicht Glück hatte und einfach keine bekommen habe.

[Blupp]

Schau dir mal den Wert von RAND_MAX an. Dieser sagt dir wie viele verschiedene Zahlen zwischen 0 und 1 du höchstens haben kannst. Wahrscheinlich gibt es da nur 2^16 Zahlen. Was die Genauigkeit von double bei weitem unterschreitet.

[Tobias Kempkens]

also das könnte natürlich der grund sein. ich bin natürlich absoluter informatik-laie und dachte, dass mit dem typen double dann auch eben wirklich 15 kommastellen zur verfügung ständen und alle werte im intervall [0;1] entsprechend gleichverteilt wären.

da ich leider nicht weiss wo ich sowas für den bereich RAND_MAX nachgucken kann oder evtl ändern hab ich einfach mal folgendes gemacht:

die gleichverteilte ZV (welche später logarithmiert wird) so bestimmt, dass wenn sie =0 ist, stattdessen den wert 0,000001 annimmt.
und ihr hattet recht. selbst mit 10 Millionen Simulation bleibt das Ergebnis stabil bzw. errechenbar. (so wie es eigentlich seien sollte.)

also wäre das problem an der stelle gelöst. (warum es ZU instabil war).

noch eine lästige frage:
könntet ihr jdm der (quasi) NULL PLAN von der ganzen materie hat (ausser dem zusammenbasteln) erklären wie ich den befehl
rand()/static_cast<double>(RAND_MAX)
dann entsprechend umwandeln könnte, damit er in etwa ZV im double format zieht?
gibts einen anderen befehl in der standard-bib? kann man den befehl iwie erweitern? muss ich mir erst ne andere bibliothek-datei besorgen?

der GRUND warum ich eigentlich die suboptimal Box-Müller-Methode genommen habe anstatt die bessere Polar-Methode ist, dass ich gehört habe, dass BM-Methode für C++ instabil wäre bei bereits zahlen NAHE null!!! und ich gerne rausfinden möchte inwiefern.

auf auf jeden fall schonmal DANKE für eure mühe!

[JaykopX]

Um Zufallszahlen zu bekommen die größer sind als RAND_MAX hab ich mal ein Beispiel gebastelt für mein RAND_MAX = 32767 = 2^15 -1:

#include <iostream>

using namespace std;

#define MAX_UINT64 18446744073709551615 //2^64-1
double RandomNumber() {
	unsigned long long retVal=0;//Rückgabewert
	retVal = rand();			//Ersten 15 bits füllen, da RAND_MAX bei mir 2^15
	for(int i=0; i<3; ++i) {		//3*15 bits nachfüllen
		retVal <<=15;			//15 Bits verschieben
		retVal |= rand();		//nächsten 15 Bits füllen
	}
	return (double)retVal/(double)(MAX_UINT64>>4);	//Wir nutzen nur 4*15=60 bits, also 4 bits "abziehen"
}

int main() {
	cout<<"Wert von RAND_MAX:"<<RAND_MAX<<'\n';
	for(int i=0; i<100; ++i) {
		cout<<RandomNumber()<<" ## ";
	}
}

Das ganze ohne Gewähr! Da könnte noch was faul dran sein, hab es nicht ausgiebig getestet. Das ganze sollte nun eine 60 bit Zufallszahl liefern, ich habe aber nicht nachgeprüft ob da durch Rundungsfehler was verloren geht. Die Idee dahinter ist einfach das man eine Zufallszahl aus 4 Zufallszahlen zusammenbastelt.

[Michael E.]

Das verkürzt die Periode auf ein Viertel und hat wahrscheinlich bescheidene Eigenschaften bzgl. Zufälligkeit. Nimm lieber nen PRNG, der für größere Zahlen ausgelegt ist.

[SeppJ]

Also für Monte Carlo ist rand() sowieso nicht gut, weil man nie weiß was dahinter steckt. Normalerweise ist das ein ganz einfacher LCG und der ist für ernsthaftes MC absolut ungeeignet. Such dir am besten einen guten, schnellen Zufallsgenerator (multiply and carry, TAUS, ...) und implementier den selbst oder nimm eine fertige Implementierung, z.B. aus der GSL. Da ist der Zahlenraum dann auch mindestens 32 Bit, du findest bestimmt auch welche mit 64-Bit.

[Tobias K.]

o.k., danke für den tipp / die tipps.
(wobei rand() für diese monte carlo simulation noch auszureichen schien. denn bei 50Mio. Simulationen war das numerische Ergebnis identisch mit dem analytischen aus dem Black-Scholes-Modell zur Optionspreisbewertung.)

kurze frage nochmal des informatik-laien. hab in der library hinsichtlich des befehls RAND_MAX folgendes gefunden:

#define RAND_MAX 0x7FFF

was genau sagt jetzt dieses 0x7FFF aus?

[SeppJ]

Tobias K. schrieb:

was genau sagt jetzt dieses 0x7FFF aus?

32767

Und das erklärt auch, warum bei dir die 0 so wahrscheinlich ist und warum ich auf dem Computer an dem ich getestet habe (RAND_MAX = 2147483647) viele Millionen Schritte ohne 0 machen konnte.

(Und es sagt, dass ich mit meiner Antwort ganz zu Anfang zu 100% Recht hatte )

[Tobias K.]

Das ist auf jeden Fall schonmal interessant zu wissen.

was genau müsste ich denn da statt 0x7FFF eingeben um auch 2147483647 mögliche zahlen zu erhalten? oder kann ich da nicht einfach so rumfummeln?

[SeppJ]

Tobias K. schrieb:

oder kann ich da nicht einfach so rumfummeln?

Dies. Der Wert ist dazu da, dich darüber zu informieren, wie rand() implementiert ist. Der Wert wird von der rand()-Implementierung selbst aber nicht benutzt. D.h. du kannst gerne was ändern, aber das Verhalten von rand() ändert sich nicht. Aber die Programme die du compilierst würden dann falsche Annahmen über RAND_MAX machen.

[Michael E.]

Tobias K. schrieb:

Das ist auf jeden Fall schonmal interessant zu wissen.

was genau müsste ich denn da statt 0x7FFF eingeben um auch 2147483647 mögliche zahlen zu erhalten? oder kann ich da nicht einfach so rumfummeln?

Du solltest da nicht rumfummeln. PRNGs zu ändern ohne sich sehr gut auszukennen, ist meistens eine sehr schlechte Idee. Schnapp dir nen anderen Generator.

Edit: Jaja, wenn der Tab zehn Minuten geöffnet ist...

[Michael E.]

SeppJ schrieb:

Der Wert wird von der rand()-Implementierung selbst aber nicht benutzt.

Die LCGs, die ich bisher gesehen haben, haben alle vor dem return ein &(1<<x) durchgeführt... Ob dabei allerdings auf die Konstante zugegriffen wird, weiß ich nicht.

[MrBesserwisser]

da ich hier keinen C++-Compiler zur Verfügung habe.

Bei mir funktioniert es astrein.

Was mir aber aufgefallen ist, dass Du für innere und äussere Schleife die Zählvariable i benutzt. Versuch es mal mit j als innerer!!!!

[SeppJ]

Michael E. schrieb:

SeppJ schrieb:

Der Wert wird von der rand()-Implementierung selbst aber nicht benutzt.

Die LCGs, die ich bisher gesehen haben, haben alle vor dem return ein &(1<<x) durchgeführt... Ob dabei allerdings auf die Konstante zugegriffen wird, weiß ich nicht.

Aber das rand() wird nicht neu compiliert, da wird der fertige Objektcode aus der Standardbibliothek gelinkt, der kennt das Compilezeitmakro gar nicht. Das RAND_MAX ist wirklich nur dazu da, um zur Compilezeit Informationen über die Laufzeitumgebung zu haben. Genauso wie man durch ändern von CHAR_BIT nicht die Größe eine Bytes der Maschine ändern kann.

[Tobias K.]

SeppJ schrieb:

Such dir am besten einen guten, schnellen Zufallsgenerator (multiply and carry, TAUS, ...) und implementier den selbst oder nimm eine fertige Implementierung, z.B. aus der GSL. Da ist der Zahlenraum dann auch mindestens 32 Bit, du findest bestimmt auch welche mit 64-Bit.

=> o.k., also wie schon gesagt, ich bin absoluter anfänger, was programmiersprachen etc. angeht.
ich hab mir jetzt GSL1.13 runtergeladen und den ordner zu den anderen bibliotheken kopiert. aber was genau mache ich nun?
ich habe versucht eine header-datei von dem ordner to "include", aber das funktioniert nicht. wie genau komme ich nun an meine "bessere" ZV ran?^^

sorry, wenn ich euch hier mit anfänger-fragen bugge.