Monte Carlo-Simulation

Monte Carlo-Methoden

Der Begriff "Monte Carlo-Methode" entstand in den 1940er Jahren, als man im Zusammenhang mit dem Bau der Atombombe die Simulation von Zufallsprozessen erstmals in größerem Stil einsetzte, um die Wechselwirkung von Neutronen mit Materie theoretisch vorherzusagen. Die Bezeichnung ist eine Anspielung auf den für Glücksspiele bekannten Ort, da die Grundlage dieser Verfahren Zufallszahlen sind, wie man sie auch mit einem Roulette-Rad erzeugen könnte. Schon damals wurde eine ganze Reihe von grundlegenden Algorithmen entwickelt, und heute zählen Monte Carlo (MC)-Methoden zu den wichtigsten numerischen (und auch nicht-numerischen) Verfahren, die sich auf viele naturwissenschaftliche, technische und medizinische Probleme mit großem Erfolg anwenden lassen. Dabei ist es gleichgültig, ob das Problem ursprünglich statistischer Natur war oder nicht, sondern man wendet die Bezeichnung auf alle Verfahren an, bei denen die Verwendung von Zufallszahlen eine entscheidende Rolle spielt. Ein Beispiel dafür - und eine der wichtigsten Anwendungen schlechthin - ist die MC-Integration, d.h. die numerische Berechnung hochdimensionaler Integrale, die ja mit Zufall überhaupt nichts zu tun haben müssen.

Wahrscheinlichkeitstheoretische Grundbegriffe

Zufallsvariablen

Ausgangspunkt für wahrscheinlichkeitstheoretische Betrachtungen sind die Begriffe "(zufälliges) Ereignis" und "Wahrscheinlichkeit".
Unter einem Elementarereignis w versteht man den möglichen Ausgang einer Messung, Beobachtung usw. (In den Naturwissenschaften verwenden wir dafür gern den Begriff "Experiment".) Die Gesamtheit der Elementarereignisse bildet den Ereignisraum W. Teilmengen A von W, die i.a. mehr als ein Elementarereignis enthalten können, werden als Ereignisse bezeichnet. So entsprechen beim Würfeln mit einem Würfel die sechs möglichen Elementarereignisse dem Auftreten einer bestimmten Augenzahl, während z.B. das Ereignis "gerade Augenzahl" die Tatsache bezeichnet, 2, 4 oder 6 Augen geworfen zu haben. Spezielle Ereignisse sind die leere Menge (das unmögliche Ereignis) und der ganze Ereignisraum (das sichere Ereignis).
Entsprechend der relativen Häufigkeit, mit der sie bei einer großen Anzahl von Wiederholungen des Experiments auftreten, definiert man für Ereignisse Wahrscheinlichkeiten P mit den Eigenschaften
0 P(A) 1
P()=0, P(W)=1
P(AB)
P(A)+P(B)
=
P(A)+P(B)    wenn     AB=
Eine Zufallsvariable ist eine Abbildung
x: W R
die jedem Elementarereignis einen reellen Zahlenwert x(w) zuordnet, z.B. die Körpergröße einer zufällig herausgegriffenen Person. Damit ist auch die Wahrscheinlichkeit definiert, daß x einen bestimmten Wert X annimmt, nämlich mittels Rückführung auf die ursprünglichen Wahrscheinlichkeiten über dem Ereignsraum W
P(x=X) = P({w|x(w)=X})
Kann die Zufallsvariable nur abzählbar viele diskrete Werte X1, X2, ... annehmen, so definiert man als Wahrscheinlichkeitsverteilung von x den Satz von Zahlen (Wahrscheinlichkeiten)
pi = P(x=Xi)
für die dann gilt
pi 0


i 
pi = 1
Ein anderer wichtiger Fall ist der, daß die Zufallsvariable kontinuierliche Werte in einem Intervall oder auf der ganzen reellen Achse annehmen kann. Die kumulative Verteilungsfunktion F(X) ist dann die Wahrscheinlichkeit, daß x einen Wert kleiner oder gleich X annimmt
F(X) = P(x X)
mit
F(-)=0    und    F(+)=1
Wenn x nur Werte in einem endlichen Intervall annimmt, dann werden diese Grenzwerte natürlich schon am Rand des Intervalls erreicht. Dazwischen ist F(X) eine monoton wachsende Funktion.

 cdf.png 
Wenn F differenzierbar ist, dann läßt sich diese Funktion als Integral
F(X) = P(x X) =
X

- 
dx f(x)
einer Wahrscheinlichkeitsdichte f(x) darstellen, die
f(x) 0



- 
dx f(x) = 1
erfüllt. Für kleine dx ist f(xdx näherungsweise gleich der Wahrscheinlichkeit, daß die Zufallsvariable einen Wert im Intervall der Größe dx um x annimmt.

 pdf.png 
Da sich die Behandlung von diskreten und kontinuierlichen Zufallsvariablen hauptsächlich dadurch unterscheidet, ob man über Wahrscheinlichkeiten summiert oder über Wahrscheinlichkeitsdichten integriert, werden im folgenden nur kontinuierliche Variablen explizit betrachtet.
Einem Elementarereignis können auch zwei oder mehrere Zufallsvariablen, x, y, ..., zugeordnet sein (z.B. Orts- und Geschwindigkeitskomponenten eines Gasmoleküls). Im Fall von kontinuierlichen Variablen wird dann das gleichzeitige Auftreten bestimmter Werte der einzelnen Variablen durch eine multivariate Wahrscheinlichkeitsdichte f(x,y,) beschrieben. Speziell ist bei zwei Variablen f(x,ydx dy wieder die Wahrscheinlichkeit, daß die erste Variable einen Wert im Intervall dx um x annimmt und gleichzeitig die zweite Variable einen Wert im Intervall dy um y.
Zwei (oder mehrere) Zufallsvariablen heißen statistisch unabhängig, wenn sich ihre gemeinsamen Wahrscheinlichkeiten bzw. Wahrscheinlichkeitsdichten nach dem Muster
f(x,y) = g(xh(y)
faktorisieren lassen. Dabei sind g(x) und h(y) die Wahrscheinlichkeitsdichten für das Auftreten von x bzw. y allein (d.h. jeweils ohne Rücksicht auf die andere Variable).

Charakterisierung von Zufallsvariablen

Die Wahrscheinlichkeitsdichte (im diskreten Fall die Wahrscheinlichkeitsverteilung) enthält die gesamte Information über eine Zufallsvariable. Oft ist die Wahrscheinlichkeitsdichte aber modellmäßig nicht bekannt oder nicht im Detail meßbar, und man versucht daher, sie durch abgeleitete, einfacher bestimmbare Größen zu charakterisieren.
Der Mittelwert (oder Erwartungswert) einer Zufallsvariablen x mit Wahrscheinlichkeitsdichte f(x) ist definiert als
x =
dx f(xx
Das ist der Wert, der sich ergeben würde, wenn man x sehr oft mißt und die Summe der Meßwerte durch die Anzahl der Messungen dividiert. Für den Mittelwert verwendet man gern das Symbol m,
m = x
Die Varianz ist ein Maß für die mögliche Abweichung (Streuung) einer Einzelbeobachtung vom Mittelwert
Var(x) = s2 = (x - x)2 =
dx f(x)(x-x)2
Durch Ausquadrieren des Integranden rechnet man leicht nach, daß auch gilt
s2 = x2- x2
Der hier verwendete Erwartungswert von x2 ist ein Beispiel für die Momente einer Wahrscheinlichkeitsdichte
Mn = xn =
dx f(xxn
Allgemein ist der Erwartungswert einer Funktion g(x) definert als
g(x) =
dx f(xg(x)
Die Korrelation zweier Variablen x und y wird durch die Kovarianz
cov(x,y) = (x-x)(y-y) = xy- xy
charakterisiert. Wenn x und y statistisch unabhängig sind, gilt, wie man ebenfalls aus den Definitionen leicht nachrechnet,
xy = xy
bzw. allgemeiner
g(x)h(y) = g(x)h(y)
Parallel zur Kovarianz wird auch der Korrelationskoeffizient
rxy = cov(x,y)

sxsy
verwendet, für den wegen der Normierung auf die Varianzen von x und y gilt
-1 rxy 1
Oft betrachtet man auch den Korrelationskoeffizienten einer Zufallsvariablen zu verschiedenen Zeitpunkten, z.B. zur Zeit i und eine Zeit k später. Wenn die statistischen Eigenschaften von x nicht von der Zeit abhängen, dann hängt diese Zeitkorrelationsfunktion nur von der Zeitdifferenz k, aber nicht von i selbst ab
Ck = rxi,xi+k = cov(xi,xi+k)

sxisxi+k
= xi xi+k- x2

x2- x2
In der obigen Form ist Ck bei k=0 auf 1 normiert. Da in den meisten Fällen xi+k nach einer gewissen Zeit von xi statistisch unabhängig ist, gehen Zeitkorrelationsfunktionen in der Regel für k nach Null.

Beispiele

Der einfachste Fall einer kontinuierlichen Zufallsvariablen ist die Gleichverteilung über einem Intervall [a,b], d.h. x nimmt alle Werte aus dem Intervall mit der gleichen Wahrscheinlichkeit
f(x) =



1

b-a
wenn    a x b
0
sonst
an. [Der Wert von f(x) im Inneren des Intervalls ergibt sich aus der Normierungsbedingung dx f(x)=1.]
Eine ganze Reihe von physikalischen Prozessen, darunter die Lebensdauern der Kerne beim radioaktiven Zerfall, werden durch die Exponentialverteilung
f(x) =



1

l
 e-x/l
für    x 0
0
sonst
beschrieben. Die Zufallsvariable x kann hier alle Werte auf der positiven rellen Achse annehmen; allerdings sind große Werte viel seltener als kleine. Der Erwartungswert von x ist durch den Parameter l gegeben.
Das wichtigste Modell einer kontinuierlichen Zufallsvariablen ist jedoch die Gauß- oder Normalverteilung
f(x) =
1

2ps2

 exp
- (x-m)2

2s2

Hier kann die Zufallsvariable alle Werte zwischen - und + auf der reellen Achse annehmen. Der Erwartungswert von x ist m, die Varianz s2. Den Standardfall m = 0, s = 1 bezeichnet man als N(0,1)-Verteilung. Die Normalverteilung spielt nicht nur in der Statistik eine zentrale Rolle, auch in der Physik sind z.B. die kartesischen Komponenten der Geschwindigkeit von Gasmolekülen normalverteilt.
Ein Beispiel für eine diskrete Wahrscheinlichkeitsverteilung ist die Binomialverteilung
Bk(n,p) =
n
k

 pk(1-p)n-k
Sie gibt an, mit welcher Wahrscheinlichkeit ein Ereignis in n unabhängigen Beobachtungen genau k-mal auftritt (k=0,1,2,,n), wenn p die Wahrscheinlichkeit ist, daß das Ereignis bei einer einzelnen Beobachtung auftritt [und daher (1-p) die Wahrscheinlichkeit, daß es bei einer einzelnen Beobachtung nicht auftritt]. Mittelwert und Varianz von k sind hier durch < k > =pn bzw. s2=p(1-p)n gegeben.
Im Limes p0, n, sodaß pnl konstant bleibt, geht die Binomialverteilung in die Poissonverteilung
Pl(k) = lk

k!
 e-l
für k=0,1,2, , mit Mittelwert < k > =l und Varianz s2=l über. Sie beschreibt das Auftreten unabhängiger "seltener" Ereignisse, wie z.B. der Anzahl der radioaktiven Zerfälle oder Molekülkollisionen in einem gegebenen Zeitintervall, wird aber auch in der Warteschlangentheorie zur Modellierung von Ereignissen wie dem Eintreffen von Anrufen in einer Telephonzentrale verwendet.

Stichproben

Sind die Wahrscheinlichkeitsverteilungen oder -dichten von Zufallsvariablen bekannt, so lassen sich daraus im Prinzip alle statistischen Kenngrößen wie Mittelwerte, Varianzen, Korrelationen usw.  berechnen - sofern man nämlich in der Lage ist, die entsprechenden mathematischen Ausdrücke analytisch oder numerisch auszuwerten. In der Praxis tritt aber häufig der Fall auf, daß man zwar ein theoretisches Modell für die Verteilung hat, aber die Parameter nicht kennt und daher aus einer Messung bestimmen will, oder man hat überhaupt noch keine Vorstellung von der Form der Verteilung. In diesem Fall kann man die statistischen Kenngrößen auch mit Hilfe einer Stichprobe schätzen.
Wir betrachten ein Zufallsexperiment, dessen Ausgang durch eine Zufallsvariable x beschrieben wird. Unter einer (unabhängigen) Stichprobe vom Umfang n
{x1,x2,,xn}
versteht man eine n-fache Wiederholung des Experiments, wobei die Ergebnisse in den einzelnen Experimenten einander nicht beeinflussen sollen. Mit anderen Worten, {x1,x2,,xn} ist ein Satz von Zufallsvariablen mit denselben statistischen Eigenschaften wie x, und xi und xj sind für i j statistisch unabhängig (xi ist einfach der Wert von x im i-ten Versuch). Die Erzeugung der Stichprobe kann selbst wieder als Zufallsexperiment aufgefaßt werden.
Erwartungswerte und andere statistische Kenngrößen werden nun durch Mittelung über die Stichprobe geschätzt. So ist z.B. der Stichprobenmittelwert definiert als
-
x
 
= 1

n
  n

i=1 
xi
Man erwartet, daß
-
x
 
x
gilt, wobei x den "wahren" (aus der eventuell unbekannten theoretischen Verteilung berechneten) Mittelwert von x bezeichnet. Der Stichprobenmittelwert ist tatsächlich erwartungstreu, denn der Erwartungswert für das Zufallsexperiment "Stichprobenmittelwert bilden" ist
-
x
 
= 1

n
 

i 
xi = 1

n
 

i 
x = x
Dies ist zunächst aber nur eine Aussage für den Fall, daß man unendlich viele Stichproben erzeugen und daraus jeweils den Stichprobenmittelwert bestimmen könnte.
Der analoge Ansatz für die Varianz
(s2) = 1

n
  n

i=1 

xi -
-
x
 

2
 
unterschätzt aber den wahren Wert um einen Faktor (n-1)/n, denn
(s2)
=
1

n
  <

i 

xi - 1

n


j 
xj

xi - 1

n


k 
xk
>
=
1

n
  <

i 
xi2 - 2

n


ij 
xi xj + 1

n2


ijk 
xj xk > = 1

n
 


i 
xi2- 1

n


ij 
xi xj
=
1

n
 
nx2- 1

n
(nx2+ n(n-1)x2)
= n-1

n
 (x2- x2)
Wenn man also als Schätzwert
s2 = 1

n-1
  n

i=1 

xi -
-
x
 

2
 
verwendet, dann ist diese Größe konstruktionsgemäß erwartungstreu, d.h. s2 = s2.
Eine vollkommen analoge Rechnung ergibt für die Varianz des Stichprobenmittelwerts,
Var(
-
x
 
) = <
-
x
 
-x
2
 
>
folgende Beziehung zur (wahren) Varianz von x
Var(
-
x
 
) = 1

n
 Var(x)
Die Wurzel daraus bezeichnet man als Standardfehler des Mittelwerts. Da Var(x) für eine gegebene Zufallsvariable x eine feste Zahl ist,1 zeigt die letzte Formel, daß auch bei Durchführung nur einer einzigen Stichprobe der Stichprobenmittelwert nahe am wahren Mittelwert liegen wird, wenn nur der Umfang der Stichprobe hinreichend groß ist.

Erzeugung von Zufallszahlen am Computer

Pseudo-Zufallszahlengeneratoren

Sowohl zur Simulation von Zufallsprozessen im eigentlichen Sinn als auch zur Lösung vieler anderer mathematischen Aufgaben ist es wünschenswert, am Computer Stichproben von Zufallsvariablen erzeugen zu können. Im Prinzip könnte man dazu natürlich "echte" Zufallsprozesse wie den Zerfall einer radioaktiven Quelle oder thermisches Rauschen verwenden. Abgesehen von instrumentellen Problemen (und wer hätte schon gern eine Strahlungsquelle auf dem Schreibtisch), würden sich auf diese Weise aber kaum Zufallszahlen in hinreichender Menge und mit der Geschwindigkeit erzeugen lassen, wie sie von vielen Anwendungen verbraucht werden. Außerdem muß es für Testzwecke möglich sein, ein und dieselbe Folge von Zufallszahlen beliebig oft zu reproduzieren.
Aus diesem Grund verwendet man in Simulationen fast ausschließlich Pseudozufallszahlen: Das sind Folgen von Zahlen, die zwar durch einen streng deterministischen Algorithmus erzeugt werden, aber "zufällig aussehen", d.h. gewisse statistische Tests auf Zufälligkeit erfüllen. (Daraus folgt auch, daß es für jeden Pseudo-Zufallszahlengenerator eine Anwendung geben kann, bei der er sich als nicht zufällig genug erweist.)
In den meisten Compilern, Interpretern und manchen Script-Sprachen ist wengstens ein (Pseudo)-Zufallszahlengenerator implementiert, der gleichverteilte Zufallszahlen aus dem Intervall [0,1) liefert. Häufig ist das ein linearer Kongruenzgenerator (Modulo-Generator), der im einfachsten Fall nach dem folgenden Prinzip arbeitet:
Es werden ganze Zahlen a, c und m sowie ein Startwert ("Seed") i0 < m vorgegeben und rekursiv die Folgen
in+1
=
(a·in + cmod m
xn+1
=
in+1/m
gebildet. Die Zahlen a, c und m sind fest und bestimmen die Eigenschaften des Zufallszahlengenerators. Nach Vorgabe eines (mit gewissen Einschränkungen) beliebigen Startwerts i0 erhält man durch sukzessives Anwenden der obigen Vorschrift zunächst eine streng deterministische Folge von ganzen Zahlen i0,i1,i2, mit 0 in < m. Die eigentlichen "Funktionswerte" des Generators sind die xn, für die nach Konstruktion 0 xn < 1 gilt.
Auf Grund der Formeln ist klar, daß die in höchstens m verschiedene Werte annehmen können; dann muß sich die Folge wiederholen. Unter bestimmten Voraussetzungen an die Parameter a, c und m kann man erreichen, daß die Periode des Generators den Maximalwert m erreicht. Dann sind die xn in ihrer Gesamtheit auf jeden Fall gleichverteilt auf dem Intervall [0,1). Je nach Wahl von a, c und m erfüllen sie auch gewisse Tests auf Zufälligkeit (statistische Unabhängigkeit). Auf keinen Fall dürfen jedoch für die Parameter des Generators x-beliebige Zahlen eingesetzt werden.
Um sich die Berechnung der Modulo-Funktion zu ersparen, wird häufig m=2r gesetzt, wo r die Wortlänge des Computers in Bits ist (bzw. m=2r-1, wenn nicht Zweierkomplement-Darstellung verwendet, sondern das Vorzeichen separat gespeichert wird). In diesem Fall werden vom Ergebnis einer Operation mit ganzen Zahlen, das einen "Overflow", also mehr als r Binärstellen, produziert, automatisch nur die niedrigsten r Bits gespeichert. Die folgende Tabelle enthält die Parameter für einige einfache Generatoren.
a c m i0
 16807   0   231-1   ungerade 
65539 0 231 ungerade
69069 1 232  
1664525 0 232 ungerade
25214903917 11 248  
Der Eintrag in der zweiten Zeile ("IBM-Generator") war seinerzeit auf Großrechnern weit verbreitet, ist aber eigentlich ein Gegenbeispiel, da die resultierenden Zufallszahlen stark korreliert sind. Faßt man nämlich Tripel von aufeinanderfolgenden Funktionswerten zu Koordinaten (xn,xn+1,xn+2) im dreidimensionalen Einheitswürfel zusammen, so kommen alle diese Punkte auf einer geringen Anzahl von Ebenen zu liegen. Solche Korrelationen treten zwar prinzipiell immer auf, bei guten Generatoren aber erst in viel höheren Dimensionen.
Einfache (in Compiler "eingebaute") Generatoren wie die obigen sind zwar für viele Zwecke ausreichend, bei Anwendungen, die kritisch von der Qualität der Zufallszahlen abhängen, empfiehlt es sich jedoch, einen gut getesteten Generator aus einer Programmbibliothek zu verwenden.

Erzeugung nicht-gleichverteilter Zufallszahlen

Während man davon ausgehen kann, daß am Computer - entweder über den eingebauten Generator oder eine Bibliotheksfunktion - standardmäßig gleichverteilte Zufallszahlen aus dem Intervall [0,1) zur Verfügung stehen, müssen in der Praxis häufig Zufallsvariablen mit einer anderen Verteilung simuliert werden.
Ist z.B. x eine diskrete Zufallsvariable, die die Werte X1,X2,,XN mit Wahrscheinlichkeiten p1,p2,,pN annimmt, so kann man
x0=0,    x1=p1,    x2=p1+p2,    x3=p1+p2+p3,    ,    xN=1
setzen und das Einheitsintervall in eine Folge von Teilintervallen
[x0,x1),    [x1,x2),    [x2,x3),   ,    [xN-1,xN)
zerlegen. Da die Wahrscheinlichkeit, in das Teilintervall [xi-1,xi) zu treffen, gleich seiner Länge xi-xi-1=pi ist, "würfelt" man also gleichverteilte Zufallszahlen x aus [0,1) und setzt
x=Xi,       wenn       xi-1 x < xi
Diese sogenannte inverse Transformationsmethode läßt sich im Prinzip auch anwenden, wenn x abzählbar unendlich viele Werte annimmt. Man muß dann nur das Intervall [0,1) in unendlich viele (immer kleiner werdende) Teilabschnitte zerlegen.
Ein weiterer einfacher Fall ist die Simulation einer gleichverteilten (kontinuierlichen) Zufallsvariablen x über dem Intervall [a,b). Hier leistet offenbar die Transformation
x  x=a+x(b-a)
wobei wieder x [0,1) gleichverteilt ist, das Gewünschte.
Die Verallgemeinerung dieser Transformationsmethode für kontinuierliche Zufallsvariablen erhalten wir, indem wir zunächst eine beliebige Transformation x y(x) betrachten, die streng monoton zunehmend (oder streng monoton abnehmend) und differenzierbar sein soll. y ist dann ebenfalls eine Zufallsvariable, und die Wahrscheinlichkeit, daß y in einem Intervall [Y1,Y2) liegt, ist
P(Y1 y < Y2) = P(X1 x < X2)
wobei X1 und X2 die Urbilder von Y1 und Y2 sind und der streng monoton zunehmende Fall angenommen wurde. Sind f(x) und g(y) die Wahrscheinlichkeitsdichten von x und y, so ergibt sich daraus

X2

X1 
dx f(x) =
Y2

Y1 
dy g(y) =
x(Y2)

x(Y1) 
dx  dy

dx
 g(y(x))
und im Limes X2X1, Y2Y1
f(x) = dy

dx
 g(y(x))
In der Form
dy

dx
= f(x)

g(y)
kann man dies als Differentialgleichung für jene Transformation ansehen, die aus x eine Zufallsvariable y mit vorgegebener Wahrscheinlichkeitsdichte g(y) macht. (Im Fall einer streng monoton abnehmenden Transformation ist die linke Seite der letzten Gleichung durch -dy/dx zu ersetzen.)
Soll beispielsweise die Exponentialverteilung mit Hife des eingebauten Zufallszahlengenerators simuliert werden, so sind
f(x)
=
1
       0 x < 1
g(y)
=
1

l
 e-y/l
       0 y <
also
dy

dx
=
l ey/l
d(y/le-y/l
=
dx
-e-y/l
=
x+c
Die Integrationskonstante ergibt sich aus der Forderung y(0)=0 zu c=-1. Die Transformation
x  y = - l ln(1-x)
macht also aus der auf [0,1) gleichverteilten Zufallsvariablen x eine exponentialverteilte Variable y auf [0,).
Der größte Nachteil der Transformationsmethode, die im wesentlichen darauf hinausläuft, die Gleichung
G(y) = F(x)
für gegebenes x nach y aufzulösen, besteht darin, daß man häufig entweder schon die Stammfunktionen zu f und g, F und G, nicht analytisch angeben oder die Gleichung G(y)=F(x) nicht nach y auflösen kann. Es gibt aber neben der Transformationsmethode einige andere allgemein anwendbare Verfahren sowie eine Unzahl von Methoden zur Simulation von Stichproben aus speziellen Verteilungen.

Monte Carlo-Integration

Eines der wichtigsten - wenn nicht das wichtigste - Monte Carlo-Verfahren ist eine Methode zur numerischen Berechnung von Integralen, insbesondere in hohen Dimensionen: die sogenannte Monte Carlo-Integration.
Angenommen, es sei das Integral einer Funktion g(x) auf dem Intervall [a,b] zu berechnen
I =
b

a 
dx g(x)
Spaltet man vom Integral die Länge des Integrationsgebietes, b-a, ab
I = (b-a)
b

a 
dx  1

b-a
 g(x)
so kann man den ersten Faktor des Integranden als (korrekt normierte) Wahrscheinlichkeitsdichte
f(x) =



1

b-a
wenn    a x b
0
sonst
einer auf [a,b] gleichverteilten Zufallsvariablen x auffassen. Damit läßt sich aber, abgesehen von einem Vorfaktor, das Integral als Erwartungswert der Funktion g der Zufallsvariablen x interpretieren
I = (b-a)
dx f(xg(x) = (b-ag(x)
Kann man das Integral bzw. den Erwartungswert nicht analytisch berechnen, so liegt es nahe, eine Stichprobe x1,x2,,xn der Variablen x zu erzeugen und den Erwartungswert von g(x) durch den Stichprobenmittelwert zu approximieren
I (b-a

g(x)
 
= (b-a 1

n
  n

i=1 
 g(xi)
Dieser Ausdruck hat große Ähnlichkeit mit konventionellen Formeln für numerische Integration, allerdings werden hier die Stützpunkte zufällig im Integrationsgebiet gewählt und nicht äquidistant oder in einer anderen regelmäßigen Anordnung.
Beim obigen, naive MC-Integration genannten Verfahren würfelt man die Stützstellen gleichverteilt im Integrationgebiet. Ein allgemeinerer Fall liegt vor, wenn das Integral von vornherein die Form
I =
dx f(xg(x)
hat, mit einer nicht-konstanten Funktion f(x) 0, für die dx f(x) < ist. Man kann dann f(x)/dx f(x) als Wahrscheinlichkeitsdichte einer nicht-gleichverteilten Zufallsvariablen x auffassen und das Integral
I =

dx f(x)
 
dx  f(x)


dx f(x)
 g(x) =

dx f(x)
  < g(x) >
wieder durch einen Stichprobenmittelwert annähern
I

dx f(x)
  1

n
  n

i=1 
 g(xi)
wobei die xi jetzt allerdings eine Stichprobe aus der Wahrscheinlichkeitsdichte f(x)/dx f(x) sein müssen. Dieses Verfahren nennt man Importance Sampling.
Der große Vorteil von MC-Integrationsverfahren besteht darin, daß bei unabhängigen Stichproben vom Umfang n die Varianz des Stichprobenmittelwerts durch
Var(
-
g
 
) = 1

n
 Var(g)
gegeben und der Fehler des Integrals daher proportional zu 1/n ist:
DI = 1

n
 

dx f(x)
 

 

g2- g2
 
Da diese Beziehung unabhängig von der Dimensionalität des Integrationsgebiets ist, sind MC-Verfahren zwar bei niedrig-dimensionalen Problemen konventionallen Integrationsverfahren unterlegen, bei hochdimensionalen Integralen stellen sie jedoch die Methode der Wahl dar.

Footnotes:

1Das setzt natürlich voraus, daß sx2 < ist.


File translated from TEX by TTH, version 3.80.
On 15 Jul 2008, 08:16.