Monte-Carlo-Simulation

Grundprinzip und mathematische Grundlage

Eine Monte-Carlo-Simulation (auch MC-Simulation oder Monte-Carlo-Studie) ist ein Verfahren der Stochastik bzw. Wahrscheinlichkeitstheorie. Dabei werden wiederholt Zufallsstichproben aus einer Verteilung mithilfe von Zufallsexperimenten gezogen. Ziel ist es, Probleme numerisch zu lösen, die analytisch nicht oder nur mit großem Aufwand lösbar sind.

Die Zufallsexperimente können real durchgeführt werden, etwa durch Würfeln, oder in einem Computer mit Monte-Carlo-Algorithmen. Diese verwenden Zufallszahlen oder Pseudozufallszahlen, um zufällige Ereignisse nachzubilden. Grundlage ist vor allem das Gesetz der großen Zahlen: Bei wachsender Zahl von Experimenten nähert sich ein empirischer Mittelwert dem entsprechenden Erwartungswert an.

Monte-Carlo-Simulationen eignen sich besonders zur Berechnung des Erwartungswertes einer Funktion f. Für eine diskrete Zufallsvariable gilt E_X[f(X)] = ∑{x∈Ω} P(x) f(x), für eine kontinuierliche Zufallsvariable E_X[f(X)] = ∫{x∈Ω} P(x) f(x) dⁿx.

Dabei ist P(x) im diskreten Fall eine Wahrscheinlichkeit und im kontinuierlichen Fall eine Wahrscheinlichkeitsdichte. In der statistischen Mechanik kann P(x) beispielsweise ein Boltzmanngewicht sein. f(x) bezeichnet den Wert der Funktion bei der Realisierung x. Ω ist der Raum der möglichen Ergebnisse, etwa der Phasenraum der Teilchen eines physikalischen Systems.

Der Erwartungswert wird durch Stichproben geschätzt. Je nach zugrunde liegender Verteilung kann der Schätzer eine größere oder kleinere Varianz besitzen. Varianz bezeichnet hier die Streuung der Schätzwerte. Ziel ist häufig, mit möglichst wenigen Stichproben einen Schätzer mit möglichst kleiner Varianz zu erhalten. Dazu dienen Verfahren zur Varianzreduktion. Stichproben aus einer Wahrscheinlichkeitsverteilung können auch mit Markov-Chain-Monte-Carlo-Verfahren erzeugt werden.

Konvergenz und wichtige Anwendungsfelder

Die Konvergenz einer Monte-Carlo-Simulation kann mit der Gelman-Rubin-Statistik bewertet werden. Mit steigender Zahl der Experimente sinkt die Varianz des Ergebnisses; die Schätzung wird dadurch im Allgemeinen zuverlässiger, bleibt aber eine Näherung.

Monte-Carlo-Methoden werden überall eingesetzt, wo Prozesse statistisches Verhalten zeigen oder viele mögliche Zufallsverläufe untersucht werden sollen. Wichtige Anwendungsgebiete sind:

• Physik, insbesondere Statistische Physik und Quantenmechanik. Observablen eines Systems im Gleichgewicht werden dort als Erwartungswerte beschrieben.
• Integralgeometrie und stochastische Geometrie.
• Mathematische Optimierung, zum Beispiel in Form evolutionärer Algorithmen.
• Numerische Schätzung der Fehlerfortpflanzung, besonders bei nichtlinearen Funktionen.
• Bayessche Statistik. Komplizierte Bayessche Modelle können an Daten angepasst und mithilfe der Posterior predictive distribution für Vorhersagen verwendet werden; ein Beispiel ist Approximate Bayesian Computation.

Monte-Carlo-Simulationen dienen außerdem der Nachbildung komplexer Prozesse. Beispiele sind Produktionsprozesse zur Aufdeckung von Engpässen und Opportunitäten, Klimamodelle, Rekonstruktionsverfahren in der Nuklearmedizin und die Risikoaggregation im Risikomanagement. In der Wertermittlung werden sie etwa zur Unternehmens- oder Immobilienbewertung sowie zur Preisbestimmung komplexer Finanzkontrakte wie exotischer Optionen eingesetzt, wenn keine analytische Bewertungsformel bekannt ist.

Weitere Beispiele sind die räumliche Verteilung des energieabhängigen Neutronenflusses in einem heterogenen Medium, etwa im Blanket eines Kernfusionsreaktors, die Simulation alternder Nuklearwaffen, die Modellierung von Regentropfenkollisionen und die Untersuchung der Kugelverteilung in den Fächern eines Galtonbretts.

Numerische Integration und Resampling

Eine zentrale Anwendung ist die Monte-Carlo-Integration, besonders bei hochdimensionalen Integralen. Klassische Integrationsalgorithmen leiden dort stark unter dem Fluch der Dimensionalität und sind häufig nicht mehr anwendbar. Hochdimensionale Integranden sind jedoch oft stark lokalisiert. MCMC-Verfahren können dann Stichproben aus einer geeigneten Verteilung erzeugen und so eine effiziente Berechnung ermöglichen.

Sei Ω ⊂ ℝⁿ eine beliebig dimensionale Menge, f eine integrierbare Funktion und V das Volumen von Ω. Wird ein gleichverteilter Zufallsvektor X mit der Wahrscheinlichkeitsdichte p(x) = 1/V verwendet, gilt S(f) := ∫_{x∈Ω} f(x) dⁿx = V E_X[f(X)].

Für gleichverteilte Stichproben X_i ∼ p lautet der Schätzer: Ŝ(f) = V · (1/n) · ∑_{i=1}ⁿ f(x_i).

Mit steigender Zahl an Stichproben nähert sich dieser Schätzer dem Integralwert beliebig genau an.

Beim Resampling werden Verteilungseigenschaften von Zufallsvariablen untersucht, deren Verteilungstyp unbekannt ist. So kann beispielsweise die nichtzentrale Verteilung eines Korrelationskoeffizienten ermittelt werden: Viele Realisierungen werden mit Zufallszahlen simuliert und anschließend in einer Häufigkeitstabelle, einer empirischen Verteilungsfunktion oder einem Histogramm zusammengefasst. Ebenso lassen sich die Eigenschaften von Schätzfunktionen bei Ausreißern untersuchen. Eine Simulation kann zeigen, dass das arithmetische Mittel dann nicht mehr der beste Schätzer für den Erwartungswert ist. Auch Verteilungsparameter können auf diese Weise geschätzt werden.

Typische Beispiele

Zur probabilistischen Bestimmung der Kreiszahl π werden zufällige Punkte gleichverteilt in [-1,+1]² gewählt. Für jeden Punkt (x,y) wird geprüft, ob er im Einheitskreis liegt, also ob x² + y² ≤ 1 gilt. Mit der Indikatorfunktion z = I(x²+y²≤1) = 1, falls x²+y²≤1, und 0 sonst, entsteht eine Bernoulli-verteilte Zufallsvariable Z ∼ Ber(p). Ihr Parameter ist das Verhältnis von Kreisfläche zu Quadratfläche: p = (r²·π)/(2r)² = π/4, wobei r = 1 ist.

Für n stochastisch unabhängige und identisch verteilte Realisierungen z₁,…,zₙ ist die relative Häufigkeit p̂ = (Anzahl der Punkte in der Kreisfläche)/(Anzahl der erzeugten Punkte im Quadrat) = (∑ᵢ₌₁ⁿ zᵢ)/n der übliche Schätzwert für p. Da E[Z] = P(Z=1) = p, kann p auch durch das arithmetische Mittel geschätzt werden. Somit ist π̂ = 4p̂ ein Monte-Carlo-Schätzwert für π. Für diese Schätzung können mithilfe des Standardfehlers Konfidenzintervalle angegeben werden, die ihre Unsicherheit ausdrücken.

Dasselbe Beispiel ist als numerische Integration formulierbar. Man setzt Ω = [-1,1]² und verwendet als Funktion die Indikatorfunktion I_Kreis des Einheitskreises. Dann gilt S(f) = V E_{X,Y}[I_Kreis(X,Y)] = π, weil V = 2·2 = 4 die Fläche des Quadrats ist. Der Schätzer lautet Ŝ(f) = V · (1/N) · ∑ᵢ₌₁ᴺ I_Kreis(xᵢ,yᵢ).

Weitere Beispiele sind die näherungsweise Auswertung der Crofton-Formel zur Bestimmung der Bogenlänge einer Kurve und der Miller-Rabin-Primzahltest. Dieser entscheidet probabilistisch, ob eine natürliche Zahl prim ist. Seine Ausgabe lautet entweder „sicher zusammengesetzt“ oder „wahrscheinlich prim“. Die Wahrscheinlichkeit, dass eine zusammengesetzte Zahl pro Durchgang als „wahrscheinlich prim“ eingestuft wird, liegt unter 25 % und kann durch wiederholte Ausführung weiter gesenkt werden. Der Test liefert keine Faktoren einer zusammengesetzten Zahl und ist daher kein Faktorisierungsverfahren.

Entstehung und technische Umsetzung

Eine frühe Anwendung war das 1733 von Georges-Louis Leclerc de Buffon vorgestellte Nadelproblem, mit dem die Kreiszahl π mithilfe des Zufalls näherungsweise bestimmt werden kann. In den 1930er Jahren entwickelte Enrico Fermi erste Ideen zu Monte-Carlo-Simulationen mit elektronischen Rechenmaschinen. Während des Manhattan-Projekts wurden solche Methoden im Los Alamos Scientific Laboratory zur Untersuchung des Neutronentransports in nuklearen Materialien eingesetzt; die Methode musste zunächst geheim gehalten werden. Weitere wichtige Beiträge leisteten unter anderem Stanisław Marcin Ulam und John von Neumann.

Als grundlegende Veröffentlichung gilt die 1953 im Journal of Chemical Physics erschienene Arbeit Equation of State Calculations by Fast Computing Machines von Nicholas Metropolis, Marshall N. Rosenbluth, Arianna W. Rosenbluth, Edward Teller und Augusta H. Teller. Untersucht wurde die Zustandsgleichung eines zweidimensionalen Systems starrer Kugeln als Modell einer Flüssigkeit. Die Simulation verwendete 224 Teilchen und periodische Randbedingungen. Es gab bis zu 48 Zyklen; in jedem Zyklus führte jedes Teilchen einen Bewegungsschritt aus. Ein Zyklus dauerte auf dem MANIAC I drei Minuten. Verwendet wurde eine Stichprobenmethode mit Wichtung durch den Boltzmannfaktor, die den Kern des MC-Verfahrens im Metropolis-Algorithmus bildet.

Der Name Monte-Carlo wurde von Nicholas Metropolis geprägt und spielt auf die Spielbank im Stadtteil Monte-Carlo des Stadtstaates Monaco an.

Bekannte Programme und Codes sind Geant4, ein vom CERN entwickeltes Simulationsprogramm, MCNP für den Teilchentransport in Reaktor- und Kernfusionstechnik, OpenMC für den Transport von Neutronen und Photonen, PYTHIA für Kollisionen in der Teilchenphysik, Serpent sowie SPICE für analoge, digitale und gemischte elektronische Schaltungen. Bei SPICE kann eine Monte-Carlo-Analyse die Auswirkungen der Streuung von Bauteilwerten innerhalb ihrer Toleranzen berechnen. Die aufgeführten Codes bilden nur eine Auswahl; daneben gibt es kommerzielle und für spezielle Prozesse entwickelte Programme.