Monte-Carlo-Simulation

Die Kreiszahl Pi wird mit der Monte-Carlo-Methode angenähert bestimmt durch das Vierfache der Wahrscheinlichkeit, mit der ein innerhalb des Quadrats zufällig gewählter Punkt in den Kreis fällt. Aufgrund des Gesetzes der großen Zahlen sinkt mit steigender Anzahl von Experimenten die Varianz des Ergebnisses.

Monte-Carlo-Simulation (auch MC-Simulation oder Monte-Carlo-Studie) ist ein Verfahren aus der Stochastik bzw. Wahrscheinlichkeitstheorie, bei dem wiederholt Zufallsstichproben einer Verteilung mithilfe von Zufallsexperimenten gezogen werden^[1].

Ziel ist es analytisch nicht (oder nur aufwendig) lösbare Probleme mithilfe der gezogenen Stichproben numerisch zu lösen. Als Grundlage ist vor allem das Gesetz der großen Zahlen zu sehen. Die Zufallsexperimente können entweder – etwa durch Würfeln – real durchgeführt werden oder in Computerberechnungen mittels Monte-Carlo-Algorithmen. Bei Monte-Carlo-Algorithmen werden zur Simulation von zufälligen Ereignissen Zufallszahlen oder auch Pseudozufallszahlen benutzt.

Zu den Pionieren der Monte-Carlo-Methode in den 1940er Jahren gehören Stanislaw Ulam, Nicholas Metropolis und John von Neumann. Als grundlegende Veröffentlichung gilt eine Arbeit von Metropolis, Edward Teller, Augusta H. Teller, Marshall Rosenbluth und Arianna W. Rosenbluth von 1953.

Geschichte

Das 1733 von Georges-Louis Leclerc de Buffon vor der Pariser Akademie der Wissenschaften vorgestellte Nadelproblem^[2], das mit Hilfe des Zufalls die näherungsweise Bestimmung der Kreiszahl Pi ermöglicht, war eine der ersten Anwendungen einer Monte-Carlo-Simulation. (→Probabilistische Bestimmung der Zahl Pi)

Enrico Fermi hatte in den 1930er Jahren die eigene Ideen zu Monte-Carlo-Simulationen mittels elektronischer Rechenmaschinen. Ausgeführt wurden diese 1946 von Stanislaw Ulam und dem von ihm deshalb kontaktierten John von Neumann.^[3] Dies geschah zur Zeit des 2. Weltkriegs während der Arbeit an einem damals geheimen Projekt am Los Alamos Scientific Laboratory, für das ein Codename nötig war. Es ging im Rahmen der Entwicklung der ersten Atombombe um die Neutronendiffusion in nuklearen Materialien.^[4] Auch die mathematische Methode der Simulation musste geheim gehalten werden. Der Name Monte-Carlo wurde von Nicholas Metropolis geprägt und hängt wie folgt mit der Methode zusammen: Stan Ulam hatte einen Onkel, der sich zum Spielen immer Geld von Verwandten geliehen hatte, denn „er musste nach Monte Carlo gehen“.^[5] Dies ist natürlich eine Anspielung auf die Spielbank Monte-Carlo im gleichnamigen Stadtteil des Stadtstaates Monaco.^[6][7][8]

Als grundlegende Veröffentlichung gilt eine Arbeit von Nicholas Metropolis, Marshall N. Rosenbluth und dessen Ehefrau Arianna W. Rosenbluth, Edward Teller und dessen Ehefrau Augusta H. Teller, veröffentlicht 1953 im Journal of Chemical Physics.^[9] Ziel war die Berechnung der Zustandsgleichung eines zweidimensionalen Systems starrer Kugeln als Modelle einer Flüssigkeit. Simuliert wurde mit 224 Teilchen und periodischen Randbedingungen. Jede Simulation bestand aus bis zu 48 Zyklen, in denen jeweils jedes Teilchen einen Bewegungsschritt ausführte. Ein Zyklus benötigte drei Minuten auf dem MANIAC I Computer des Los Alamos National Laboratory. Verwendet wurde eine Sampling-Methode mit Wichtung über den Boltzmannfaktor, das Herzstück des MC-Verfahrens im Metropolis-Algorithmus, wobei die Idee nach Marshall Rosenbluth von Teller gekommen sein soll.^[10][11] Nach Rosenbluth leisteten er und seine Frau die Hauptarbeit für den Artikel (Metropolis hätte hauptsächlich Computerzeit zur Verfügung gestellt) und sie waren die einzigen der Autoren, die das Verfahren in anschließenden Publikationen weiterverfolgten,^[12][13] sie wandten sich aber selbst ebenfalls bald darauf anderen Forschungsthemen (Plasmaphysik) zu.

Mathematik

Mathematisch ist das System ein wahrscheinlichkeitsgewichteter Weg im Phasenraum (allgemein Zustandsraum). Monte-Carlo-Simulationen sind besonders geeignet, um statistische Mittelwerte einer Größe $ {\mathcal {A}} $,

$ \left\langle {\mathcal {A}}\right\rangle =\sum _{x\in \Omega }P(x)\,{\mathcal {A}}(x), $

oder hochdimensionale Integrale (Monte-Carlo-Integration) wie

$ \int _{x\in \Omega }\!\!P(x)\,{\mathcal {A}}(x)\;\mathrm {d} ^{n}x $

zu berechnen. $ P(x) $ soll in diesem Zusammenhang ein normiertes statistisches Gewicht (etwa ein Boltzmanngewicht) sein. $ {\mathcal {A}}(x) $ ist der Wert der Größe $ {\mathcal {A}} $ im Zustand $ x $. Die Summation bzw. Integration verläuft hier über einen Raum $ \Omega $, also der Phasenraum der Teilchen im System.

Häufig ist der Raum $ \Omega $ so groß, dass die Summation nicht vollständig durchgeführt werden kann. Stattdessen erzeugt man nun eine Markow-Kette $ x_{1},x_{2},x_{3},\ldots $ von Zuständen in $ \Omega $, deren Häufigkeit wie das vorgegebene Gewicht $ P(x) $ verteilt ist. Bereiche des Raumes $ \Omega $ mit hohem Gewicht sollen also häufiger in der Markow-Kette vertreten sein als Bereiche mit niedrigem Gewicht. Man spricht hier von Importance Sampling. Gelingt dies, so lassen sich die Erwartungswerte einfach als arithmetisches Mittel der Größe $ {\mathcal {A}} $ zu diesen Zuständen der Markow-Kette berechnen, also als

Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): \left\langle\mathcal{A}\right\rangle \approx \frac{1}{N}\sum_{i=1}^N \mathcal{A}(x_i).

Dieser Zusammenhang basiert auf dem Gesetz der großen Zahlen. Je nach physikalischem System kann es schwierig sein, diese Markow-Kette zu erzeugen. Insbesondere ist sicherzustellen, dass die Markow-Kette tatsächlich den gesamten Raum $ \Omega $ bedeckt und nicht nur einen Teil des Raumes abtastet. Man sagt: der Algorithmus muss ergodisch sein.

Quasi-Monte-Carlo-Simulationen verwenden keine Pseudozufallszahlen, sondern eine Sequenz mit geringer Diskrepanz (zum Beispiel eine Sobol-Sequenz).

Anwendungen

Mit der Monte-Carlo-Methode können Probleme mit statistischem Verhalten simuliert werden. Diese Methode hat deshalb besonders in der Physik wichtige Anwendungen gefunden.

Heutige Supercomputer (HPC) basieren auf massivem Multiprocessing mit vielen tausend Einzelprozessoren, die parallel arbeiten. Diese Gegebenheiten lassen sich besonders gut mit solchen probabilistischen Lösungsverfahren ausnutzen.^[14]

Anwendungen der Monte-Carlo-Simulation sind beispielsweise folgende.

Berechnung von Integralen

• Bestimmung von (höherdimensionalen) Integralen, siehe Beispiel

Resampling

Untersuchung der Verteilungseigenschaften von Zufallsvariablen unbekannten Verteilungstyps,

• Beispiel: die Ermittlung der nichtzentralen Verteilung des Korrelationskoeffizienten. Mit Hilfe von Zufallszahlen wird die Realisierung beliebig vieler Korrelationskoeffizienten simuliert. Eine Zusammenfassung der Koeffizienten in eine Häufigkeitstabelle ergibt eine empirische Verteilungsfunktion oder ein Histogramm.
• die Eigenschaften von Schätzfunktionen bei Vorliegen von Ausreißern in Daten. Mit Hilfe der Simulation kann gezeigt werden, dass das arithmetische Mittel nicht mehr ein bester Schätzer für den Erwartungswert ist.
• die Schätzung von Verteilungsparametern.

Nachbildung von komplexen Prozessen

• Produktionsprozesse in einem Fertigungsunternehmen, um Engpässe und Opportunitäten in der Produktion aufzudecken
• Klimamodelle
• Rekonstruktionsverfahren in der Nuklearmedizin.
• Risikoaggregation zur Bestimmung des Gesamtrisikoumfangs eines Unternehmens im Risikomanagement^[15]
• Ableitung von Bewertungen in der Wertermittlung, z. B. bei der Unternehmensbewertung oder Immobilienwirtschaft^[16][17]. Bepreisung komplexer Finanzkontrakte wie „exotische“ Optionen, bei denen keine analytische Formel für die Bewertung eines Finanzproduktes bekannt ist.
• räumliche Verteilung des energieabhängigen Neutronenflusses in einem heterogenen Medium, etwa im Blanket eines Kernfusionsreaktors.
• Supercomputer und MC Methoden werden u. a. für die Simulation der alternden Nuklearwaffen (siehe auch Stockpile stewardship) der USA benutzt.^{[18][19][20][21]}
• Wege eines einzelnen Regentropfens simulieren, der mit zufällig verteilten anderen Tropfen kollidiert. Nach der Simulation mehrerer konkreter Tropfen sind Aussagen über die durchschnittliche Tropfengröße möglich oder auch zu Temperatur und Tröpfchendichte, bei denen Schnee oder Hagel entstehen.
• Verteilung der Kugeln auf die Fächer beim Galtonbrett.

Beispiele

Probabilistische Bestimmung der Zahl Pi

Man wählt hierzu zufällige Punkte $ \left(x,y|x\in \left[-1..1\right]\wedge y\in \left[-1..1\right]\right) $ aus und überprüft (durch Anwendung des Satzes von Pythagoras), ob diese innerhalb des Einheitskreises liegen:

$ x^{2}+y^{2}\leq 1 $.

Über das Verhältnis der Anzahl der Punkte innerhalb und außerhalb des Kreises kann dann folgendermaßen $ \ \pi $ bestimmt werden:

$ {\frac {{\text{Kreisfl}}\mathrm {\ddot {a}} {\text{che}}}{{\text{Quadratfl}}\mathrm {\ddot {a}} {\text{che}}}}={\frac {r^{2}\cdot \pi }{(2\cdot r)^{2}}}={\frac {\pi }{4}}={\frac {{\text{Treffer in Kreisfl}}\mathrm {\ddot {a}} {\text{che}}}{\text{generierte Punkte im Quadrat}}}=P\left({\text{Im Kreis}}\right) $

Numerische Integration

Das obige Beispiel zur Bestimmung von Pi bildet praktisch das Flächenintegral einer Viertelkreisfläche. Entsprechend kann man das Flächenintegral allgemeiner, auch höherdimensionaler Funktionen nach dieser Methode berechnen. Soll das Integral

Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): S(f)=\int_0^1 f(x) \, dx

einer Funktion $ f $ berechnet werden, dann wählt man Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): m unabhängige im Intervall Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): [0,1] gleichverteilte Punkte Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): x_1,\dots, x_m und approximiert Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): S(f) durch

$ S_{m}(f)={\frac {1}{m}}\sum _{i=1}^{m}f(x_{i}). $

Im allgemeineren Fall von höherdimensionalen Funktionen ist das Vorgehen ähnlich. Sei Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): K\subset\mathbb{R}^n eine beliebige Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): n -dimensionale Menge und $ f\colon K\rightarrow \mathbb {R} $ eine integrierbare Funktion. Um den Wert

$ S(f)=\int _{K}f(x)\,dx $

näherungsweise zu berechnen, wählt man zufällig in der Menge Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): K gleichverteilte Punkte $ x_{i} $ für $ i=1,\dots ,m $. Dann approximiert

$ S_{m}(f)={\frac {\mathrm {vol} (K)}{m}}\sum _{i=1}^{m}f(x_{i}) $

den Wert $ S(f) $ in Abhängigkeit von $ m $ beliebig genau. Um wie oben vorgestellt Pi zu berechnen, muss man $ K:=[-1,1]^{2} $ und $ f:=\chi _{\mathrm {Kreis} } $ als charakteristische Funktion des Einheitskreises wählen. Hier ist $ S(\chi _{\mathrm {Kreis} })=\pi $ gerade die Fläche des Einheitskreises.

In der Praxis werden Monte-Carlo Verfahren vor allem für die Berechnung hochdimensionaler Integrale verwendet. Hier sind klassische Integrationsalgorithmen stark vom Fluch der Dimensionalität betroffen und nicht mehr anwendbar. Allerdings sind speziell hochdimensionale Integranden meist stark lokalisiert^[22]. In diesen Fällen erlauben insbesondere MCMC-Verfahren die Erzeugung von Stichproben mit einer Verteilung die eine effiziente Berechnung solcher hochdimensionaler Integrale erlaubt.

Miller-Rabin-Primzahltest

Ein Beispiel für eine Monte-Carlo-Simulation ist der Miller-Rabin-Test, bei dem probabilistisch bestimmt wird, ob eine natürliche Zahl prim ist oder nicht. Die Ausgabe des Tests lautet entweder „sicher zusammengesetzt“ oder „wahrscheinlich prim“. Die Wahrscheinlichkeit, dass eine zusammengesetzte Zahl als „wahrscheinlich prim“ klassifiziert wird, liegt pro Durchgang unter 25 % und kann durch mehrfache Ausführung weiter gesenkt werden. Der Miller-Rabin-Test liefert keine Aussage über die Faktoren einer zusammengesetzten Zahl, ist also kein Faktorisierungsverfahren.

Programmpakete (Auswahl)

• MCNP, der Monte-Carlo N-Particle Transport Code, ist ein prototypisches, weltweit verbreitetes reaktorphysikalisches Programm, das sehr häufig angewendet wird, auch in der Kerntechnik und der Kernfusionstechnik.^[23] Die aktuelle Version ist MCNP6.2. Auf der MCNP-Webseite^[24] sind Handbücher und Release Notes als Internetdokumente zu finden, zum Beispiel der Band I des MCNP-Handbuchs Overview and Theory.^[25]

• PYTHIA ist ein Simulationsprogramm für die Teilchenphysik und simuliert Kollisionen und dabei entstehende Teilchen.

• SPICE ist ein Simulationsprogramm für analoge, digitale und gemischte elektronische Schaltungen. Mit der Monte-Carlo-Simulation ist es möglich, die Auswirkungen der Streuung der Bauteilewerte innerhalb der angegebenen Toleranz zu berechnen.

Siehe auch

• Binder-Kumulante

Literatur

• Kurt Binder: Monte Carlo methods in statistical physics. Springer, Berlin [u. a.] 1979, ISBN 3-540-09018-5.
• Kurt Binder: Applications of the Monte Carlo method in statistical physics. Springer, Berlin 1984, ISBN 3-540-12764-X.
• Paul Glasserman: Monte Carlo Methods in Financial Engineering. Springer, Berlin 2003, ISBN 978-1-4419-1822-2.
• David P. Landau, Kurt Binder: A guide to Monte Carlo simulations in statistical physics. Cambridge University Press, Cambridge, 2014, ISBN 978-1107074026

Weblinks

Commons: Monte-Carlo-Simulation – Sammlung von Bildern, Videos und Audiodateien

• Interaktive Animation der Monte-Carlo-Simulation zur Abschätzung der Kreiszahl π.
• Die Monte-Carlo-Simulation für die Teilchenphysik. Film des Karlsruher Instituts für Technologie zur Monte-Carlo-Simulation in der Teilchenphysik auf YouTube

Einzelnachweise