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 des Verfahrens Zufallszahlen sind, wie man sie auch mit einem Roulette-Rad erzeugen könnte. Schon damals wurde eine ganze Reihe von grundlegenden Verfahren 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(AÈB)

P(A)+P(B)

P(A)+P(B) wenn AÇB=Æ

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 X₁, X₂, ... annehmen, so definiert man als Wahrscheinlichkeitsverteilung von x den Satz von Zahlen (Wahrscheinlichkeiten)

p_i = P(x=X_i)

für die dann gilt

p_i ³ 0

å
i

p_i = 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.

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(x) dx näherungsweise gleich der Wahrscheinlichkeit, daß die Zufallsvariable einen Wert im Intervall der Größe dx um x annimmt.

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 Variable 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,y) dx 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(x) h(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(x) x

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) = s² = á(x - áxñ)²ñ =

ó
õ

dx f(x)(x-áxñ)²

Durch Ausquadrieren des Integranden rechnet man leicht nach, daß auch gilt

s² = áx²ñ- áxñ²

Der hier verwendete Erwartungswert von x² ist ein Beispiel für die Momente einer Wahrscheinlichkeitsdichte

M_n = áxⁿñ =

ó
õ

dx f(x) xⁿ

Allgemein ist der Erwartungswert einer Funktion g(x) definert als

ág(x)ñ =

ó
õ

dx f(x) g(x)

Die Korrelation zweier Variablen x und y wird durch die Kovarianz

cov(x,y) = á(x-áxñ)(y-áyñ)ñ = áyñ- áxñáyñ

charakterisiert. Wenn x und y statistisch unabhängig sind, gilt, wie man ebenfalls aus den Definitionen leicht nachrechnet,

áxyñ = áxñáyñ

bzw. allgemeiner

ág(x)h(y)ñ = ág(x)ñáh(y)ñ

Parallel zur Kovarianz wird auch der Korrelationskoeffizient

r_xy =

cov(x,y)

s_xs_y

verwendet, für den wegen der Normierung auf die Varianzen von x und y gilt

-1 £ r_xy £ 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

C_k = r_{x_i,x_i+k} =

cov(x_i,x_i+k)

s_{x_i}s_{x_i+k}

áx_i x_i+kñ- áxñ²

áx²ñ- áxñ²

In der obigen Form ist C_k bei k=0 auf 1 normiert. Da nach einer gewissen Zeit x_i+k von x_i 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 an

f(x) =

ì
ï
í
ï
î

b-a

wenn a £ x £ b

sonst

[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) =

ì
ï
í
ï
î

e^-x/l

für x ³ 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) =

æ
è

2ps²

ö
ø

exp

é
ë

(x-m)²

2s²

ù
û

Hier kann die Zufallsvariable alle Werte zwischen -¥ und +¥ auf der reellen Achse annehmen. Der Erwartungswert von x ist m, die Varianz s². Den Standardfall m = 0, s = 1 bezeichet 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

B_m(n,p) =

æ
è

n
m

ö
ø

p^m(1-p)^n-m

Sie gibt an, mit welcher Wahrscheinlichkeit ein Ereignis in n unabhängigen Beobachtungen genau m-mal auftritt (m=0,1,2,¼,n). Dabei ist p die Wahrscheinlichkeit, daß das Ereignis in einer einzelnen Beobachtung auftritt.

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

{x₁,x₂,¼,x_n}

versteht man eine n-fache Wiederholung des Experiments, wobei die Ergebnisse in den einzelnen Experimenten einander nicht beeinflussen sollen. Mit anderen Worten, {x₁,x₂,¼,x_n} ist ein Satz von Zufallsvariablen mit denselben statistischen Eigenschaften wie x, und x_i und x_j sind für i ¹ j statistisch unabhängig (x_i 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

n
å
i=1

x_i

Man erwartet, daß

» á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

ñ = á

å
i

x_iñ =

å
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

(s²)¢ =

n
å
i=1

æ
è

x_i -

ö
ø

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

á(s²)¢ñ

å
i

æ
è

x_i -

å
j

x_j

ö
ø

æ
è

x_i -

å
k

x_k

ö
ø

å
i

x_i² -

å
ij

x_i x_j +

n²

å
ijk

x_j x_k > =

é
ë

å
i

áx_i²ñ-

å
ij

áx_i x_jñ

ù
û

é
ë

náx²ñ-

(náx²ñ+ n(n-1)áxñ²)

ù
û

n-1

(áx²ñ- áxñ²)

Wenn man also als Schätzwert

s² =

n-1

n
å
i=1

æ
è

x_i -

ö
ø

verwendet, dann ist diese Größe konstruktionsgemäß erwartungstreu, d.h. ás²ñ = s².

Eine vollkommen analoge Rechnung ergibt für die Varianz des Stichprobenmittelwerts,

Var(

) = <

æ
è

-áxñ

ö
ø

folgende Beziehung zur (wahren) Varianz von x

Var(

) =

Var(x)

Die Wurzel daraus bezeichnet man als Standardfehler des Mittelwerts. Da Var(x) für eine gegebene Zufallsvariable x eine feste Zahl ist, 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. In der Regel 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") i₀ < m vorgegeben und rekursiv die Folgen

i_n+1

(a·i_n + c) mod m

x_n+1

i_n+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 i₀ erhält man durch sukzessives Anwenden der obigen Vorschrift zunächst eine streng deterministische Folge von ganzen Zahlen i₀,i₁,i₂,¼ mit 0 £ i_n < m. Die eigentlichen "Funktionswerte" des Generators sind die x_n, für die nach Konstruktion 0 £ x_n < 1 gilt.

Auf Grund der Formeln ist klar, daß die i_n 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 x_n 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=2^r gesetzt, wo r die Wortlänge des Computers in Bits ist (bzw. m=2^r-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	i₀
16807	0	2³¹-1	ungerade
65539	0	2³¹	ungerade
69069	1	2³²
1664525	0	2³²	ungerade
25214903917	11	2⁴⁸

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 (x_n,x_n+1,x_n+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 X₁,X₂,¼,X_N mit Wahrscheinlichkeiten p₁,p₂,¼,p_N annimmt, so kann man

x₀=0, x₁=p₁, x₂=p₁+p₂, x₃=p₁+p₂+p₃, ¼, x_N=1

setzen und das Einheitsintervall in eine Folge von Teilintervallen

[x₀,x₁), [x₁,x₂), [x₂,x₃), ¼, [x_N-1,x_N)

zerlegen. Da die Wahrscheinlichkeit, in das Teilintervall [x_i-1,x_i) zu treffen, gleich seiner Länge x_i-x_i-1=p_i ist, "würfelt" man also gleichverteilte Zufallszahlen x aus [0,1) und setzt

x=X_i, wenn x_i-1 £ x < x_i

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 der 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 [Y₁,Y₂) liegt, ist

P(Y₁ £ y < Y₂) = P(X₁ £ x < X₂)

wobei X₁ und X₂ die Urbilder von Y₁ und Y₂ 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

ó
õ

X₂

X₁

dx f(x) =

ó
õ

Y₂

Y₁

dy g(y) =

ó
õ

x(Y₂)

x(Y₁)

g(y(x))

und im Limes X₂®X₁, Y₂®Y₁

f(x) =

g(y(x))

In der Form

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)

0 £ x < 1

g(y)

e^-y/l

0 £ y < ¥

also

l e^y/l

d(y/l) e^-y/l

-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

b-a

g(x)

so kann man den ersten Faktor des Integranden als (korrekt normierte) Wahrscheinlichkeitsdichte

f(x) =

ì
ï
í
ï
î

b-a

wenn a £ x £ b

sonst

einer auf [a,b] gleichverteilten Zufallsvariablen x auf fassen. 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(x) g(x) = (b-a) ág(x)ñ

Kann man das Integral bzw. den Erwartungswert nicht analytisch berechnen, so liegt es nahe, eine Stichprobe x₁,x₂,¼,x_n der Variablen x zu erzeugen und den Erwartungswert von g(x) durch den Stichprobenmittelwert zu approximieren

I » (b-a)

g(x)

= (b-a)

n
å
i=1

g(x_i)

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(x) g(x)

mit einer nicht-konstanten Funktion f(x) ³ 0 hat. Man kann dann f(x)/òdx f(x) als Wahrscheinlichkeitsdichte einer nicht-gleichverteilten Zufallsvariablen x auf fassen und das Integral

I =

é
ë

ó
õ

dx f(x)

ù
û

ó
õ

f(x)

ó
õ

dx¢ f(x¢)

g(x) =

é
ë

ó
õ

dx f(x)

ù
û

< g(x) >

wieder durch einen Stichprobenmittelwert annähern

I »

é
ë

ó
õ

dx f(x)

ù
û

n
å
i=1

g(x_i)

wobei die x_i 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

Var(

) =

Var(g)

und der Fehler des Integrals daher proportional zu 1/Ön ist

DI =

Ön

é
ë

ó
õ

dx f(x)

ù
û

ág²ñ- ágñ²

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.

File translated from T_EX by T_TH, version 3.59.
On 21 Jun 2005, 14:18.