Monte-Carlo Simulati..

Monte-Carlo Simulation
Sehr häufig hängen wichtige Ergebnisse von unbekannten Werten wesentlich ab, für die man
allerhöchstens statistische Daten hat oder für die man ein Modell der Wahrscheinlichkeitsrechnung
annimmt.
So sind z.B. Betonfestigkeiten, Bodenfestigkeiten NIE exakt vorhersagbar, dennoch hängt die
Sicherheit der Konstruktion stark von diesen Parametern ab.
Manchmal ist auch das mathematische Modell zu komplex, um selbst bei bekannten Daten das
Ergebnis zu errechnen.
Für solche Fälle bietet sich die Monte Carlo Simulation an. Man modelliert das ganze System als
Computermodell und lässt den Computer extrem viele Beispiele abarbeiten. Aus der statischen
Erfassung der Ergebnisse lassen sich so Rückschlüsse auf die gesuchten Werte machen:
z.B. ein Würfelspiel:
A hat 2 Würfel, B hat einen Würfel. A schlägt B, wenn einer seiner Würfel mehr Augen anzeigt als der
Würfel von B, andernfalls verliert er. Wie gut sind seine Chancen?
Dieses Problem lässt sich mit der Wahrscheinlichkeitsrechnung bzw. Kombinatorik exakt lösen. Man
kann aber ebenso leicht den Computer einige 1000 Mal gegen sich selbst würfeln lassen.
Ist p die Gewinnwahrscheinlichkeit von A und q jene von B, so gilt p+q = 1. Das Gesetz der großen
Zahlen sagt, dass A unter N Spielen mit hoher Wahrscheinlichkeit (p N) gewinnen wird. Hat er also
in der Simulation genau a gewonnen, so folgt
a ~ Np oder p ~ a/N
Der Computer würfelt also N mal, zählt die Anzahl a der Siege von A und dividiert a/N.
Das Problem ist nur noch: Wie würfelt man in C?
Zufallszahlen:
Das Würfeln ist ein Zufallsexperiment mit genau 6 möglichen Ergebnissen. Im Idealfall wird jede
Augenzahl mit derselben Wahrscheinlichkeit gewürfelt werden, d.h. wir haben eine
"Gleichverteilung" der Zahlen 1 bis 6 zu simulieren.
Ein getürkter Würfel würde z.B. die Zahl 6 wesentlich öfter würfeln und hätte daher eine andere
Verteilung.
Andere Experimente liefern "kontinuierliche" Zufallszahlen, d.h. es sind alle Zahlen aus Intervallen
oder gar alle reellen Zahlen zugelassen.
Gleichverteilung in [0,1]
Das ist die wichtigste Verteilung in der Simulation. Sie simuliert ein Experiment, in dem jede Zahl aus
[0,1] mit gleicher Wahrscheinlichkeit herauskommt. So eine Verteilung lässt sich nur mittels der
Verteilungsfunktion F(x) oder deren Ableitung f(x) (der Dichtefunktion) beschreiben:
F(x) = Wahrscheinlichkeit, dass das Ergebnis kleiner (oder gleich) x ist. Bei der Gleichverteilung in
[0,1]:
F(x) = 0 falls x <= 0, F(x) = 1 falls x >= 1, und F(x) = x für 0 <= x <= 1
f(x) = 1, falls 0 < x < 1 und 0 sonst.
Simulation einer Gleichverteilung, Zufallszahlen
Zufallszahlen kann man in C so erzeugen:
int rand(void) liefert eine Gleichverteilung der ganzen Zahlen 0,1,2, ..., RAND_MAX
also liefert
double X = (double) rand()/RAND_MAX)
(ACHTUNG!!! das double ist nötig, da es sich sonst um eine int-Division handeln würde und 0
herauskommt) eine Approximation für eine Gleichverteilung in [0,1]. Manchmal wird auch
double X = rand()/(RAND_MAX + 1.0)
oder auch etwas ähnliches wie
double X = rand()/(RAND_MAX + 0.000001)
empfohlen, wodurch man sich das double spart und eine Gleichverteilung in [0,1) erhält. Ich
verwende eigentlich fast immer die letzte Methode.
RAND_MAX ist ebenso wie rand() in stdlib.h deklariert.
Würfeln in C
wir brauchen eine Gleichverteilung von 1,2,3,4,5,6 was man am leichtesten so erzielt: Ist X die
Gleichverteilung in [0,1):
so "würfeln" wir 1, falls 0 <= X < 1/6 gilt. 0 <= 6X < 1,
1 <= 6X+1 < 2
so "würfeln" wir 2, falls 1/6 <= X < 2/6 gilt 1 <= 6X < 2,
2 <= 6X+1 < 3
...
Der Ausdruck 6X+1 liegt aber im Intervall [1,7) und daher ist diese Vorschrift
sehr leicht durch die Zuweisung an eine int-Variable erledigt:
int w = 6*X+1;
Ganz allgemein: Sucht man eine Gleichverteilung der Zahlen n, n+1, n+2, ..., n + (k-1) (k Zahlen), so
bekommt man das mit der Formel
int w = n + k*X;
zustande. So lassen sich auch Kopf/Zahl Spiele oder Roulette etc. simulieren.
Nicht gleichverteilte Wahrscheinlichkeiten gehen etwas komplizierter und erfordern eine nicht
triviale Modifikation obiger Formeln. Um eine gegebene diskrete Verteilung Y der n Zahlen 0, 1, 2, …,
(n-1) mit vorgegebenen Wahrscheinlichkeiten P(Y = i) = pi , deren Summe natürlich 1 sein muss, zu
generieren, geht man so vor:
1) berechne die Teilsummen s0 = p0, s1 = p0+p1, … , sn-1 = p0+p1+…+pn-1 = 1
2) Generiere Stichproben x einer Gleichverteilung X in [0,1) s.o.
3) Finde den Index i, sodass si ≤ x < si+1
4) Gib y = i zurück
Hier muss man im Allgemeinen eine Menge Vergleiche abarbeiten, ehe man das Resultat y = i erhält.
Für den einfachen Fall n = 2 (z.B. Münzwurf oder Erfolg/Misserfolg) mit gegebenen P(Y = 0) = p und
P(Y = 1) = q, p+q = 1 ist obiges Verfahren einfach zu implementieren. Die Teilsummen sind nämlich
nur p und 1 und wir erhalten:
1) Generiere Stichproben x einer Gleichverteilung X in [0,1) s.o.
2) Gilt x < p, dann gib 0 zurück, sonst 1
Simulation von kontinuierlichen Zufallsgrößen
Sehr oft braucht man auch Stichproben von kontinuierlichen Zufallsgrößen mit bekannter Verteilung.
Der C-Zufallszahlen-Generator ist mit Einschränkungen auch dazu zu gebrauchen. Mit
double X = rand() / (RAND_MAX + 0.000001);
erzeugt man eine annähernde Gleichverteilung in [0,1) mit derselben Verteilungsfunktion. Leider sind
solche Gleichverteilungen in der Praxis eher selten. Hier braucht man normalverteilte Zufallszahlen
oder in der Warteschlangentheorie meist exponentialverteilte Zufallsgrößen. Erstere sind relativ
schwierig zu erzeugen, letztere gehen einfach in C:
Mathematischer Background: Ist Y eine Zufallsgröße mit Verteilungsfunktion F(x), d.h.
P (Y < a) = F(a)
dann ist X = F(Y) gleichverteilt in [0,1):
P(X < a) = P( F(Y) < a) = P ( Y < F-1(a) ) = F( F-1(a) ) = a
Umgekehrt erhält man durch Anwendung der Umkehrfunktion von F aus einer [0,1)-Gleichverteilung
X die gewünschte Verteilung: Y = F-1(X)
Exponential-Verteilung:
Diese Verteilung wird fast immer für zufällige Wartezeiten in der Warteschlagentheorie verwendet
(z.B. Dauer einer Bearbeitung, Brenndauer einer Glühbirne, Funktionsdauer eines technischen
Bauteils usw.). Für die Exponential-Verteilung ist
F(x) = 1 - exp(- λ * x)
0
für x > 0
für x<= 0
(es können keine negativen Werte auftreten).
daher ist die Umkehrfunktion
F-1 (x) = - (1/ λ) log(1 - x)
log ist der natürliche Logarithmus (oft auch ln)
Die Exponential-Verteilung hat den Erwartungswert (1/λ) und die Varianz (1/λ)2. Wichtig ist, dass für
x nie der Wert 1 möglich ist, d.h. X muss eine Gleichverteilung in [0,1), und nicht [0,1], besitzen.
Beispiel: Erstelle einen zufälligen Zeitplan von Busabfahrten, die im Schnitt alle 5 Minuten stattfinden
sollen. Als Verteilung ist eine Exponentialverteilung anzunehmen:
Rechnet man in Minuten, so sucht man Stichproben einer Exponentialverteilung mit Erwartungswert
5 = 1/λ. Diese erzeugt man mittels:
double wz(double ew)
{
double x;
// zufällige Wartezeit mit Erwartungswert ew
x = rand()/(RAND_MAX + 0.000001);
return -ew*log(1-x);
// x [0,1)-gleichverteilt
}
und konstruiert einen Zeitplan wie folgt:
double start = 8*60; // 8 Uhr 00
double ende = 22*60; // um 22 Uhr schließt der Busservice
int busse = 0;
for
(double zeit = start; zeit <= ende; zeit += wz())
printf("Abfahrt Bus Nr. %d um %lf\n", ++busse, zeit);
Rechnet man in Sekunden, sind die Werte für start, ende und λ entsprechend zu ändern!
Andere Verteilungen
Leider sind die Verteilungsfunktionen F(x) von vielen anderen Verteilungen nicht so einfach zu
invertieren. Hier liefert die Numerik gute Näherungsverfahren.
Initialisierung des Zufallszahlengenerators
Der Zufallszahlengenerator (den wir durch den Aufruf von rand()) benutzen, liefert PseudoZufallszahlen (d.h. die Zahlen werden mit einem Algorithmus erstellt, verhalten sich aber wie echt
zufällige Werte), die aber stets dieselbe Folge von Zahlen durchlaufen. Will man mit anderen Werten
starten, sollte man vor dem 1. Aufruf den Zufallszahlengenerator neu „seed“en, d.h. mit anderen
Anfangsdaten belegen. Dazu dient die Funktion void srand(int seed);
Bei gleichem seed Wert wiederholt sich natürlich auch hier die Zufallsfolge. Das kann man
verhindern, indem man z.B. die time() Funktion aus <time.h> verwendet. Unter Linux kann
man auch den Seedwert aus /dev/random lesen usw.

Download Report