8. Monte-Carlo-Methoden, Teil 2 – Monte-Carlo-Integration, Teil 2 und Random Walk¶

Im vorangegangenen Notebook haben wir uns mit Monte-Carlo-Methoden im Allgemeinen beschäftigt. Am Ende der Einheit haben wir außerdem das Konzept der Monte-Carlo-Integration gestreift, und zwar mit der Berechnung einer Näherung der Zahl $\pi$ über zufällig verteilte Punkte auf dem Einheitsquadrat. Das ist jedoch nicht die typische Art und Weise, Monte-Carlo-Integration zu betreiben.

Obwohl man mit einer Integration die Fläche unter der Kurve einer Funktion berechnet, ist das komplett analoge Beispiel mit der Berechnung von $\pi$ unter der Kurve des Viertelkreises, der im Einheitsquadrat eingeschrieben ist, tatsächlich mehr ein anschauliches Beispiel für eine MC-Simulation. Die MC-Integration verwendet allerdings meist ein anderes Prinzip.

8.1 Das Integral als Erwartungswert über eine Zufallsvariable¶

Konkret haben wir in der vergangenen Einheit den Wert des folgenden Integrals simuliert: $$4\int_0^1 dx \sqrt{1-x^2}$$ Das Resultat war eine Näherung der Zahl $\pi$, eine Fehlerabschätzung inklusive (weil wir mehrere Samples verwendet haben. Muss man das Integral allerdings tatsächlich so umständlich berechnen (also viele Punkte auf die ganze Fläche werfen, dann bestimmen welche unter dem Funktionswert liegen und welche darüber, alles abzählen und das Verhältnis bilden) oder geht das einfacher?

Es geht tatsächlich einfacher. Der Trick ist, das Integral, das wir gerade aufgeschrieben hatten, als Erwartungswert einer Funktion zu interpretieren. Um das noch etwas besser zu sehen, schreiben wir das Integral nochmals auf, und zwar zunächst einmal so: $$\int_{x_{min}}^{x_{max}} dx h(x)$$ Wir haben hier noch nicht viel gemacht, außer die Integrationsgrenzen mit zwei Variablen zu benennen, die Funktion allgemein als $h$ zu bezeichnen und den (irrelevanten) Vorfaktor wegzulassen. Der nächste Schritt ist allerdings der wesentliche. Wir schreiben das gleiche Integral jetzt als: $$\int_{x_{min}}^{x_{max}} dx f(x) g(x)$$ Auch das ist keine Hexerei, denn wir haben einfach nur die Funktion $h$ als ein Produkt von zwei anderen Funktionen $f$ und $g$ geschrieben. Warum ist das jetzt so viel besser als vorher?

Die Antwort auf diese Frage ist, dass uns nun die Interpretation als Erwartungswert leichter fällt. Sagen wir, $f$ sei eine Funktion, von der wir den Erwartungswert berechnen. Und $g$ ist eine Wahrscheinlichkeitsverteilung, und zwar der Zufallsvariablen $x$, deren Werte zwischen ${x_{min}}$ und ${x_{max}}$ liegen. Dann definiert dieses Integral tatsächlich den Erwartungswert der Funktion $f$ über die $g$-verteilte Zufallsvariable $x$.

Alles schön und gut, aber wie veranstalten wir damit nun eine MC-Integration?

8.2. Näherung des Erwartungswerts durch Sampling¶

Wir nähern das Integral, also den Erwartungswert durch Sampling der Zufallsvariablen $x$. Konkret müssen wir dazu folgendes tun:

• Den eigentlichen Integranden in ein Produkt von zwei Funktionen $f$ und $g$ zerlegen
• Eine der beiden (in unserem Beispiel $g$) als Wahrscheinlichkeitsverteilung interpretieren (und verwenden können)
• Wir sampeln $N$ Werte $x_i$ für $x$, und zwar aus der Verteilung $g$
• Wir bestimmen die Summe $$\frac{1}{N} \sum_{i=1}^N f(x_i)$$
• Wir betreiben wieder Statistik mit mehreren Samples und bestimmen Mittelwerte und Standardabweichung

Soweit, so gut. Sehen wir uns jetzt nochmals das Integral an, das zum Wert von $\pi$ führt: $$4\int_0^1 dx \sqrt{1-x^2}$$ Hier kann man z.B. einfach $f(x)=\sqrt{1-x^2}$ und $g(x)=1$ (also gleichverteilte Zufallszahlen) nehmen. Somit ergibt sich für die Näherung des Integrals die Summe $$\frac{1}{N} \sum_{i=1}^N \sqrt{1-x_i^2}$$ mit gleichverteilten $x_i$.

Probieren wir das mal aus.

8.3 Näherungswert für $\pi$, revisited¶

Zunächst lohnt es sich, der Übersichtlichkeit halber ein paar Schritte des Prozesses als Funktionen zu definieren. Die Schritte, die wir hier machen, sind:

• Die Funktion, um die es geht
• Die Summe (eigentlich der Mittelwert) als Näherung des Erwartungswertes
• Das Wiederholen des Aufrufs für mehrere Samples und die damit verbundene Statistik

Was wir hier nicht als Argument einer Funktion einbauen, ist die spezifische Wahrscheinlichkeitsverteilung, die wir verwenden wollen. Die Gleichverteilung wird stattdessen erstmal fix eingebaut.

Als nächstes können wir nun diese Funktionen zusammensetzen und damit eine Simulation bzw. Integration starten. Die Ergebnisse sehen wir uns in Abhängigkeit von der Anzahl der verwendeten Punkte pro Summe bzw. Erwartungswert an. Nachdem der Loop durchgelaufen ist, werden auch noch die Ergebnisse entsprechend ausgegeben.

8.4 Verwendung der Wahrscheinlichkeitsverteilung¶

Bisher haben wir noch nicht wirklich Gebrauch von der Möglichkeit gemacht, die Wahrscheinlichkeitsverteilung eigentlich ziemlich beliebig aus dem Integral zu separieren. Aber auch dafür wollen wir uns jetzt ein Beispiel ansehen. Dabei geht es um das Integral $$\int_0^\infty dx\; x^3 e^{-x^2}$$

Um es zu mit MC-Integration zu berechnen, nehmen wir es wie folgt auseinander: $f(x)=x^3$ und $g(x)=e^{-x^2}$. Dabei stellen wir fest, dass $g$ noch nicht ganz einer Wahrscheinlichkeitsverteilung entspricht. Genauer gesagt, hat eine Normalverteilung mit Mittelwert $0$ die folgende Form $$\frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{x^2}{2\sigma^2}}$$ Ein anderes Problem ist, dass wir nur eine Seite der Normalverteilung brauchen, weil das Integral von $0$ bis $\infty$ reicht. Wir nehmen also stattdessen die sogenannte Halbnormalverteilung der Form $$\frac{\sqrt{2}}{\sigma\sqrt{\pi}} e^{-\frac{x^2}{2\sigma^2}}$$

Damit das zusammenpasst, müssen wir also $\sigma=1/\sqrt{2}$ setzen, und es ergibt sich die Verteilung $$\frac{2}{\sqrt{\pi}} e^{-x^2}$$ Damit haben wir also als korrekte Zerlegung $f(x)=\frac{\sqrt{\pi}}{2} x^3$ und $g(x)=\frac{2}{\sqrt{\pi}} e^{-x^2}$ gefunden. Und damit ergibt sich für die Näherung des Integrals die Summe $$\frac{1}{N} \sum_{i=1}^N \sqrt{\pi}x_i^3 / 2$$ mit halbnormalverteilten $x_i$.

Zur Kontrolle sehen wir uns einmal das Integral an, wie SymPy es berechnet:

Als nächstes machen wir nun die numerische Abschätzung via MC-Integration. Dafür habe ich von oben nochmal die entsprechenden Zellen kopiert und die Funktion angepasst, sowie die Wahrscheinlichkeitsverteilung ausgetauscht.

Mit dieser Technik lassen sich im Prinzip beliebige Integrale numerisch nähern. Was wir hier nicht gemacht haben, was aber diesbezüglich wichtig ist: Die Methode funktioniert nicht nur genauso gut in mehr als einer Dimension, sie ist sogar besonders effizient, je höher die Anzahl der Dimensionen wird. Dazu noch ein Detail, das Sie sich gut merken können (auch, weil wir es in der vergangenen Einheit explizit gezeigt haben):

Der Fehler bei der Monte-Carlo-Integration skaliert wie der Kehrwert der Wurzel aus $N$, der Anzahl der verwendeten Punkte bzw. der Sample-Größe. Dieser Wert hängt nicht von der Dimension des Integrals ab, während das für andere numerische Integrationsmethoden sehr wohl der Fall ist (z.B. für Gauss-Quadratur-Verfahren muss man eine Anzahl Integrationspunkte in jeder Dimension festlegen). Wenn Sie sich also nur eine Eigenschaft bei Monte-Carlo-Methoden merken, dann merken Sie sich das:

Der Fehler skaliert wie $\frac{1}{\sqrt{N}}$.

8.5 Komplexeres Beispiel für Monte-Carlo-Simulation: Random Walk¶

Nachdem Sie nun mit der spezifischen Technik der MC-Integration vertraut sind, wenden wir uns wieder der allgemeinen Simulation zu. Insbesondere sehen wir uns ein Problem an, dem Sie vielleicht auch einmal im echten Leben begegnen werden: Dem Random Walk.

Damit ist folgendes gemeint: Wir stellen uns folgende Situation bzw. folgendes System vor:

• Ein Mensch/Hund/Roboter, nennen wir ihn “agent”, kann Schritte ausführen
• Die Schritte können in der Größe fixiert sein, oder durch eine Wahrscheinlichkeitsverteilung bestimmt werden
• Die Richtung eines Schrittes kann ebenfalls eingeschränkt sein (z.B. entlang der Koordinatenachsen auf einem Gitter) oder in beliebiger Richtung ausgeführt werden können
• Die Richtung wird für jeden Schritt zufällig über eine Wahrscheinlichkeitsverteilung gewählt
• Jeder Schritt wird von der Position ausgeführt, an der der vorige Schritt geendet hat
• Die Bewegung, die so entsteht, ist ein zusammenhängender Weg
• Simuliert man viele solcher Wege, dann erhält man verschiedene Informationen, wie immer über die statistische Auswertung der Ergebnisse

Somit haben wir das Grundgerüst, das eigentlich sehr übersichtlich und ziemlich geradlinig ist. Deshalb wollen wir uns das gleich einmal konkret ansehen.

Nachdem die Funktion für einen Schritt in mehreren Walks definiert ist, werden wir nun eine Runde von Walks laufen lassen. Dafür brauchen wir jeweils eine Anzahl von Schritten und Walks, und dann einen Loop, bei dem die Schritte ausgeführt und die Positionen mitgeschrieben werden.

Nachdem wir jetzt die Daten erzeugt und auch bereits visualisiert haben, wollen wir sie noch etwas auswerten. Dafür können wir einfach das NumPy-Array mit den Positionen verwenden und bestimmte Dinge berechnen (und dann ebenfalls grafisch darstellen). Interessant sind z.B.

• Mittlere Entfernung vom Ursprung in Abhängigkeit von der Anzahl der Schritte
• Maximale Entfernung vom Ursprung in Abhängigkeit von der Anzahl der Schritte
• Richtung am letzten Punkt im Walk

Sehen wir uns das alles einfach der Reihe nach an:

8.6 Übungsaufgabe: Experimentieren mit dem Random Walk¶

Für die Untersuchungen von Random Walks gibt es viele Möglichkeiten. Experimentieren Sie nun mit diesem Problem. Kopieren Sie den Code von oben und testen Sie je nach Zeit und Möglichkeit (und unter anderem, je nachdem, was Sie sonst noch alles verändern möchten) folgende Änderungen/Szenarien:

• Ändern Sie die Schrittweite auf eine andere Konstante. Was passiert?
• Wählen Sie die Schrittweite zufällig aus einer Wahrscheinlichkeitsverteilung Ihrer Wahl. Was passiert?
• Verpassen Sie dem Random Walk eine “Schlagseite” (“Biased Random Walk”), d.h. machen Sie z.B. die Schrittweite von der Richtung abhängig. Was passiert?
• Schränken Sie die Richtung auf ein Gitter ein (also die Richtung auf z.B. die 4 Winkel 0, 90, 180, 270 Grad). Was passiert?
• Erhöhen Sie die Anzahl der Schritte und beobachten Sie die Werte für mittlere und maximale Distanz vom Ursprung. Können Sie eine bestimmte Abhängigkeit dieser Größen von der Anzahl der Schritte feststellen?