Transcript
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 verwendet. Unter Linux kann man auch den Seedwert aus /dev/random lesen usw.
