Allgemeine Verfahren

Der n-dimensionale Gleichverteilungs- Generator.
In diesem Abschnitt werden wir eine systematische Darlegung der Verfahren zur Simulation von zufälligen Veränderlichen bringen. In den Kapiteln 1 und 3 haben wir schon des öfteren Simulationsmethoden von zufälligen Veränderlichen diskutiert, die alle auf dem folgenden Prinzip beruhten: die Definition einer zufälligen Veränderlichen wurde Wort für Wort in ein Rechnerprogramm übersetzt, wobei wir den experimentellen Versuch durch den Aufruf eines Generators für Pseudozufallszahlen ersetzt hatten. Wir erinnern in diesem Zusammenhang an die Simulation der Bernoulli- Verteilung sowie der Multinominalverteilung in Kap.1.7. Auch die Simulation der $\chi^{2}$- Verteilung in Kap.3.2 kann in diesem Zusammenhang erwähnt werden. In allen diesen Beispielen mußten wir für die Erzeugung einer Zufallszahl der gesuchten Verteilung im allgemeinen viele Aufrufe des Pseudozufallszahlen Generators durchführen. Eine Verbesserung der Generatoren liegt daher auf der Hand, nämlich die Erzeugung von $N$ Zahlen in einem Aufruf. Als Beispiel benutzen wir eine verbesserte Version des RANMAR- Generators, den wir im folgenden mit RANMARN bezeichnen werden. Die Initialisierungsroutine heißt entsprechend RMARINN.

Die Methode der inversen Verteilungsfunktion
Sei $p(x)$ irgendeine Dichtefunktion einer kontinuierlichen Veränderlichen $x$. Die integrale Verteilungsfunktion

$$P(\xi) = \int_{-\infty}^{\xi} dx p(x)$$ (1)

ist eine in dem Wertebereich $0 < P(\xi) \leq 1$ definierte Funktion. Wir erzeugen eine gleichverteilte Zufallszahl $u \in [0,1]$ und setzen

$$u = P(\xi).$$ (2)

Dann ergibt die Umkehrung,

$$\xi = P^{-1}(u),$$ (3)

eine gemäß der Dichtefunktion $p(x)$ verteilte Veränderliche $\xi$.

Die direkte Anwendung der inversen Verteilungsfunktion ist zwar außerordentlich beschränkt, da nur für wenige integrale Verteilungsfunktionen die Umkehrung analytisch darstellbar ist. Trotzdem bildet, wie wir im einzelnen noch sehen werden, die Methode der inversen Verteilungsfunktion die Grundlage fast aller Simulationsmethoden von zufälligen Veränderlichen.

Beispiel. Ein einfaches Beispiel ist die Exponentialfunktion

$$p(x) = \frac{1}{\tau} e^{-x/\tau}, \; \; \; x \geq 0.$$

Die integrale Verteilungsfunktion ist

$$P(\xi) = \int_{0}^{\xi} dx p(x) = 1- e^{-\xi /\tau}.$$

Die Gleichung

$$P(\xi) = u$$

läßt sich leicht nach $\xi$ auflösen und ergibt

$$\xi= P^{-1}(u) = -\tau ln (1-u) \equiv -\tau ln (u).$$

Im letzten Schritt haben wir die Tatsache benutzt, daß, wenn $u$ eine im Intervall $[0,1]$ gleichverteilte Zufallszahl ist, so ist auch $1-u$ gleichverteilt in $[0,1]$.

Beispiel. Ein zweites einfaches Beispiel ist die Dichtefunktion

$$p(x) = x e^{-x^{2}/2} , \; \; \; \; \; x \geq 0.$$

Wegen

$$P(\xi) = \int_{0}^{\xi} dx x e^{-x^{2}/2} = 1 - e^{-\xi^{2}/2}$$

und

$$P(\xi) = 1 - e^{-\xi^{2}/2} = u$$

ist

$$\xi = \sqrt{-2 ln(1-u)} \equiv \sqrt{-2 ln(u)}$$

ein Algorithmus zur Simulation der Verteilung $p(x)$.

Methode der Tabellen- Darstellung
Eine Methode, die für kontinuierliche wie für diskrete Veränderliche gleichermaßen anwendbar ist, ist die Tabellen- Darstellung der integralen Verteilungsfunktion. Wir teilen die reelle Achse in $M$ Teilintervalle mit den Stützstellen $\xi^{(0)}, \xi^{(1)},...,\xi^{(M)}$ ein und berechnen die Wahrscheinlichkeiten

$$P_{\mu} = \int_{\xi^{(\mu-1)}}^{\xi^{(\mu)}} dx p(x).$$ (4)

Die integrale Verteilungsfunktion an der Stelle $\xi^{(n)}$ ist dann

$$P(\xi^{(n)}) = \sum_{\mu=1}^{n} P_{\mu} = \int_{-\infty}^{\xi^{(n)}} dx p(x).$$ (5)

Wie bisher muß die Gleichung

$$P(\xi) = u$$

nach $\xi$ aufgelöst werden. Zu diesem Zweck erzeugen wir eine Zufallszahl $u_{1} \in [0,1]$ und bestimmen zunächst den Wert $n$ so, daß

$$P(\xi^{(n-1)}) < u_{1} \leq P(\xi^{(n)}).$$ (6)

Anschließend erzeugen wir eine zweite Zufallszahl $u_{2}\in[0,1]$ und interpolieren gemäß

$$\xi = \xi^{(n-1)} +( \xi^{(n)} - \xi^{(n-1)}) u_{2}.$$ (7)

$\xi$ ist dann eine Veränderliche der Dichtefunktion $p(x)$. Diese Methode kann offensichtlich auch für jede diskrete Veränderliche benutzt werden. Die integrale Verteilungsfunktion ist

$$P(\kappa) = \sum_{k=0}^{\kappa} P(k).$$ (8)

Mit Hilfe einer gleichverteilten Zufallszahl $u \in [0,1]$ bestimmen wir $\kappa$ gemäß

$$P(\kappa-1) < u \leq P(\kappa).$$ (9)

$\kappa$ ist dann eine Veränderliche aus der Wahrscheinlichkeitsverteilung $P(k)$.

Man beachte, daß wir mit dem gleichen Funktionssymbol $P$ die Wahrscheinlichkeitsverteilung und die integrale Verteilung dargestellt haben. Die Unterscheidung treffen wir mit Hilfe des Arguments. $P(k)$ ist die Wahrscheinlichkeitsverteilung und $P(\kappa)$ die integrale Verteilung. Diese Schreibweise werden wir auch im folgenden beibehalten. Eine sehr viel detailliertere Behandlung der Tabellen- Methode werden wir noch als Vorbemerkung zur Simulation von diskreten Veränderlichen bringen.

Mehrdimensionale Veränderliche
Ist eine Dichtefunktion mehrerer Veränderlicher vorgegeben, so lassen sich zwei Fälle unterscheiden. Falls alle Veränderlichen unabhängig voneinander sind,

$$p(x_{1},x_{2},...,x_{n}) = \prod_{\nu=1}^{n} p_{\nu}(x_{\nu}),$$ (10)

so ergibt sich keine neue Situation. Die integrale Verteilungsfunktion ist einfach

$$P(\xi_{1},\xi_{2},...,\xi_{n}) = \prod_{\nu=1}^{n} P_{\nu}(\xi_{\nu}),$$ (11)

und die Lösung ist

$$\xi_{\nu} = P^{-1}(u_{\nu}), \; \; \; \; \; \nu=1,2,...,n,$$ (12)

wobei die $u_{\nu}, (\nu=1,2,...,n)$ gleichverteilte Zufallszahlen aus $[0,1]$ sind.

Sofern andererseits die Veränderlichen $x_{\nu}$ abhängig voneinander sind, kann man auf zwei verschiedenen Wegen vorgehen. Einmal kann versucht werden, mit Hilfe einer umkehrbar eindeutigen Transformation

$$\vec{y} = \vec{\psi}(\vec{x})$$ (13)

eine Dichtefunktion $p(\vec{y})$ unabhängiger Veränderlicher $y_{\nu}$ $(\nu=1,2,...,n)$ zu erhalten. Die Simulation des Vektors $\vec{y}$ wird dann wie oben geschildert durchgeführt. Die Rücktransformation

$$\vec{x} = \vec{\psi}^{-1}(\vec{y}) \equiv \vec{\phi}(\vec{y})$$ (14)

ergibt dann den zufälligen Vektor $\vec{x}$.

Oder aber, und das ist meist der bequemere Weg, man drückt die Dichtefunktion durch bedingte Randverteilungen aus,

$$p(x_{1},x_{2},...,x_{n}) = p_{1}(x_{1}) p_{2}(x_{2}\vert x_{1}) \cdot \cdot \cdot p_{n}(x_{n}\vert x_{1},x_{2},...,x_{n-1}),$$ (15)

mit

$$p_{\nu}(x_{\nu}\vert x_{1},x_{2},...,x_{\nu-1}) = \frac{\int... ... \cdot \cdot dx_{n} p(x_{1},...,x_{\nu-1},x_{\nu},...,x_{n})}.$$ (16)

Offensichtlich führt dann der folgende Algorithmus zum Ziel: Man bilde die integralen Verteilungsfunktionen,

$$P_{1}(\xi_{1}) = \int_{-\infty}^{\xi_{1}} dx_{1} p_{1}(x_{1}),$$

$$P_{2}(\xi_{2}\vert\xi_{1}) = \int_{-\infty}^{\xi_{2}} dx_{2} p_{2}(x_{2}\vert\xi_{1}),$$

$$.$$ (17)

$$.$$

$$P_{n}(\xi_{n}\vert\xi_{1},\xi_{2},...,\xi_{n-1}) = \int_{-\infty}^{\xi_{n}} dx_{n} p_{n}(x_{n}\vert\xi_{1},\xi_{2},...,\xi_{n-1}),$$

erzeuge einen Vektor $\vec{u}=(u_{1},u_{2},...,u_{n})$ von gleichverteilten Zufallszahlen aus $[0,1]$ und löse der Reihe nach die Gleichungen

$$P_{1}(\xi_{1}) = u_{1}$$

$$P_{2}(\xi_{2}\vert\xi_{1}) = u_{2}$$

$$.$$ (18)

$$.$$

$$P_{n}(\xi_{n}\vert\xi_{1},\xi_{2},...,\xi_{n-1}) = u_{n}.$$

Dieses ist offensichtlich kein Gleichungssystem, sondern eine Rekursion. Für eine möglichst einfache Lösung dieser Rekursion ist häufig die Reihenfolge der Veränderlichen wichtig. Für zwei Veränderliche z.B. sind die beiden folgenden Darstellungen equivalent:

$$p(x_{1},x_{2})=p_{1}(x_{1})p_{2}(x_{2}\vert x_{1})=p_{2}(x_{2})p_{1}(x_{1}\vert x_{2}).$$ (19)

Beispiel. Als Beispiel zur vorherigen Bemerkung betrachten wir die Dichtefunktion

$$p(x_{1},x_{2})= \{ \begin{array}{ll} 6 x_{1} & \; \; \; \; \... ...\; x_{1},x_{2} \geq 0 \\ 0 & \; \; \; \; \; sonst. \end{array}$$

In der Darstellung

$$p(x_{1},x_{2}) = p_{1}(x_{1}) p_{2}(x_{2}\vert x_{1})$$

erhalten wir

$$p_{1}(x_{1}) = \int_{0}^{1-x_{1}} dx_{2} p(x_{1},x_{2}) = 6x_{1} (1-x_{1}),$$

$$p_{2}(x_{2}\vert x_{1}) = \frac{p(x_{1},x_{2})}{p_{1}(x_{1})} = \frac{1}{1-x_{1}},$$

und

$$P_{1}(\xi_{1}) = \int_{0}^{\xi_{1}} dx_{1} p_{1}(x_{1}) = 3\xi_{1}^{2} -2\xi_{1}^{3},$$

$$P_{2}(\xi_{2}\vert\xi_{1}) = \int_{0}^{\xi_{2}} dx_{2} p_{2}(x_{2}\vert x_{1}=\xi_{1}) = \frac{\xi_{2}}{1-\xi_{1}}.$$

Die Rekursion

$$3\xi_{1}^{2} - 2\xi_{1}^{3} = u_{1},$$

$$\frac{\xi_{2}}{1-\xi_{1}} = u_{2}$$

beinhaltet die Lösung einer kubischen Gleichung, die zwar möglich ist, aber dem Rechner viel Arithmetik abverlangt.

In der zweiten Darstellung

$$p(x_{1},x_{2}) = p_{2}(x_{2})p_{1}(x_{1}\vert x_{2})$$

erhalten wir dagegen

$$p_{2}(x_{2}) = \int_{0}^{1-x_{2}} dx_{1} p(x_{1},x_{2}) = 3(1-x_{2})^{2},$$

$$p_{1}(x_{1}\vert x_{2}) = \frac{p(x_{1},x_{2})}{p_{2}(x_{2})} = \frac{2x_{1}}{(1-x_{2})^{2}}$$

und

$$P_{2}(\xi_{2}) = \int_{0}^{\xi_{2}} dx_{2} p_{2}(x_{2}) = 1-(1-\xi_{2})^{3},$$

$$P_{1}(\xi_{1}\vert\xi_{2}) = \int_{0}^{\xi_{1}} dx_{1} p_{1}(x_{1}\vert x_{2}=\xi_{2}) = \frac{\xi_{1}^{2}}{(1-\xi_{2})^{2}}.$$

Die daraus folgende Rekursion

$$(1-\xi_{2})^{3} = u_{1},$$

$$\frac{\xi_{1}^{2}}{(1-\xi_{2})^{2}} = u_{2}$$

ist wesentlich einfacher zu lösen als die kubische Gleichung bei der ersten Darstellung, nämlich durch

$$\xi_{2} = 1 - \sqrt[3]{u_{1}},$$

$$\xi_{1} = (1-\xi_{2}) \sqrt{u_{2}}.$$

Beispiel. Das Paradebeispiel für die Transformation auf unabhängige Veränderliche ist die Simulation der $n$- dimensionalen allgemeinen Normalverteilung:

$$p(\vec{x}) = \frac{\sqrt{det A}}{(2\pi)^{n/2}} e^{-\frac{1}{2} (\vec{x}-\vec{a})^{T} A (\vec{x}-\vec{a})}.$$

Wir schreiben die Kovarianzmatrix $A^{-1}$ in der Form

$$A^{-1} = \left( \begin{array}{cccccc} \sigma_{11} & \sigma_{... ...& \sigma_{n2} & . & . & . & \sigma_{nn} \end{array} \right)$$

mit

$$\sigma_{\mu \nu} = <(x_{\nu}-a_{\nu})(x_{\mu}-a_{\mu})>.$$

Aus der numerischen Mathematik entnehmen wir die Behauptung, daß diese Matrix als Produkt zweier Dreiecksmatrizen ausgedrückt werden kann,

$$A^{-1} = C C^{T}$$

mit

$$C = \left( \begin{array}{cccccc} c_{11} & 0 & . & . & . & 0 ... ...& \\ c_{n1} & c_{n2} & . & . & . & c_{nn} \end{array} \right).$$

Die Matrix $C$ erhält man am einfachsten aus der Rekursionsrelation

$$c_{ij} = \frac{\sigma_{ij}-\sum_{k=1}^{j-1} c_{ik}c_{jk}} {\left( \sigma_{jj} - \sum_{k=1}^{j-1} c_{jk}^{2} \right)^{1/2}}.$$

Transformiert man nun die Veränderliche $\vec{x}$ gemäß

$$\vec{y} = C^{-1} (\vec{x}-\vec{a})$$

auf die Veränderliche $\vec{y}$, so erhält man die unabhängige Normalverteilung

$$p(\vec{y}) = \frac{1}{(2\pi)^{n/2}} e^{-\frac{1}{2} \vec{y}^{2}}$$

mit Mittelwerten $y_{\nu}=0$ und Kovarianzmatrix $<y_{\nu} y_{\mu}> = \delta_{\nu \mu}.$ Die Rücktransformation ist denkbar einfach, und zwar

$$\vec{x} = C \vec{y} + \vec{a}.$$

Die eigentliche Rechenarbeit besteht offensichtlich in der Berechnung der Matrix $C$.

Wir wollen dieses Verfahren noch einmal für den Fall einer zweidimensionalen Dichtefunktion mit $<x_{1}>=<x_{2}>=0$ explizit aufschreiben. Die Kovarianzmatrix lautet

$$A^{-1} = \left( \begin{array}{cc} <x_{1}x_{1}> & <x_{1}x_{2}... ...o \sigma_{1} \sigma_{2} & \sigma_{2}^{2} \end{array} \right) ,$$

wobei wir die Varianzen bzw Kovarianzen

$$\sigma_{11} = \sigma_{1}^{2} = <x_{1}x_{1}>,$$

$$\sigma_{12} = \sigma_{21} = \rho \sigma_{1} \sigma_{2},$$

$$\sigma_{22} = \sigma_{2}^{2} = <x_{2}x_{2}>$$

eingeführt haben. $\rho$ ist der hinreichend bekannte Korrelationskoeffizient. Für die Elemente der Matriz $C$ erhalten wir der Reihe nach

$$c_{11} = \sigma_{1},$$

$$c_{21} = \rho \sigma_{2},$$

$$c_{22} = \sigma_{2} \sqrt{1-\rho^{2}},$$

oder, übersichtlicher als Matriz geschrieben,

$$C = \left( \begin{array}{cc} \sigma_{1} & 0 \\ \rho \sigma_{2} & \sigma_{2} \sqrt{1-\rho^{2}} \end{array} \right) .$$

Die Dichtefunktion für $y_{1}$ und $y_{2}$ ist

$$p(y_{1},y_{2}) = \frac{1}{2\pi} e^{-\frac{1}{2}(y_{1}^{2}+y_{2}^{2})}$$

und die Rücktransformation nach $x_{1}$ und $x_{2}$ lautet

$$x_{1} = \sigma_{1} y_{1},$$

$$x_{2} = \rho \sigma_{2} y_{1} + \sigma_{2} \sqrt{1-\rho^{2}} y_{2}.$$

Dieses Ergebnis stimmt bis auf die Reihenfolge der Veränderlichen mit den Ergebnissen in Kap.1.5 überein. Wegen der Symmetrie in den Veränderlichen $x_{1}$ und $x_{2}$ kommt es auf die Reihenfolge nicht an.

Reduktionsmethoden
Besonders wichtig für die Simulation eindimensionaler Veränderlicher sind die Reduktionsmethoden. Hierbei geht man von einer mehrdimensionalen Dichtefunktion aus, für die man bereits eine Simulationsverfahren besitzten möge. Beispiele hierfür sind die $n$- dimensionale Verteilung gleichverteilter unabhängiger Zufallszahlen aus $[0,1]$ oder die $n$- dimensionale Normalverteilung. Man sucht dann geeignete Transformationen,

$$\vec{y} = \vec{\psi}(\vec{x}),$$ (20)

wobei die Randverteilung einer bestimmten Komponente $y_{j}$ dann genau der Verteilung der gesuchten Dichtefunktion entsprechen möge. Die Integration zur Bildung von Randverteilungen ist in der Simulation trivial. Man braucht alle übrigen Veränderlichen außer $y_{j}$ einfach nur zu ignorieren.

Wir beginnen mit dem einfachsten Fall, der Reduktion einer zweidimensionalen Veränderlichen auf eine eindimensionale Veränderliche. Sei also $p(x_{1},x_{2})$ eine Dichtefunktion, für die wir bereits ein Verfahren zur Simulation der Veränderlichen $x_{1}$ und $x_{2}$ besitzen mögen. Wir transformieren gemäß

$$y_{1} = \psi_{1}(x_{1},x_{2}) ,$$ (21)

$$y_{2} = \psi_{2}(x_{1},x_{2}) ,$$ (22)

und erhalten eine Dichtefunktion $p(y_{1},y_{2})$ für die Veränderlichen $y_{1}$ und $y_{2}$. Die Randverteilungen

$$p_{1}(y_{1}) = \int_{-\infty}^{\infty} dy_{2} p(y_{1},y_{2}),$$ (23)

$$p_{2}(y_{2}) = \int_{-\infty}^{\infty} dy_{1} p(y_{1},y_{2})$$ (24)

sind nun häufig gesuchte Dichtefunktionen, deren Simulation damit auf die Simulation der Dichtefunktion $p(x_{1},x_{2})$ zurückgeführt ist.

Beispiel. Als Beispiel diskutieren wir die zweidimensionale unabhängige Normalverteilung mit Mittelwerten $<x_{1}>=<x_{2}>=0$ und Varianzen $\sigma_{1}=\sigma_{2}=1$,

$$p(x_{1},x_{2}) = \frac{1}{2\pi} e^{-\frac{1}{2}(x_{1}^{2}+x_{2}^{2})}$$

Die Transformation

$$y_{1} = x_{1}+x_{2}$$

$$y_{2} = \frac{x_{1}}{x_{2}}$$

führt auf die Dichtefunktion:

$$p(y_{1},y_{2}) = \frac{1}{2\pi}\frac{\vert y_{1}\vert}{(1+y_... ...e^{-\frac{1}{2} \frac{(1+y_{2}^{2})y_{1}^{2}}{(1+y_{2})^{2}}},$$

und die Integration über $y_{1}$ liefert

$$p_{2}(y_{2}) = \frac{1}{\pi} \frac{1}{(1+y_{2}^{2})}.$$

Diese Dichtefunktion nennt man die Cauchy- Verteilung. Man erhält also eine Veränderliche der Cauchy- Verteilung, indem man die Veränderlichen $x_{1}$ und $x_{2}$ aus einer zweidimensionalen unabhängigen Normalverteilung mit $<x_{1}>=<x_{2}>=0$ und Varianzen $\sigma_{1}=\sigma_{2}=1$ durcheinander dividiert:

$$y_{2} = \frac{x_{1}}{x_{2}}.$$

Beispiele. Weitere Beispiele sind die Verteilungsfunktionen, die man aus der Reduktion von $n$ unabhängigen gleichverteilten Zufallszahlen erhält. Seien also $u_{\nu} \in [0,1], \; \nu=1,2,...,n$, und setzen wir

$$y_{1} = \sum_{\nu=1}^{n} u_{\nu} ,$$

$$y_{2} = max \{ u_{1},u_{2},...,u_{n} \} ,$$

so erhält man die Randverteilungen

$$p_{1}(y_{1}) \approx \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}(y_{1}-n/2)^{2}},$$

$$p_{2}(y_{2}) = n y_{2}^{n-1}.$$ (25)

Die erste Verteilung haben wir in vorherigen Abschnitten schon mehrmals diskutiert, die zweite Verteilung haben wir bei der Behandlung geordneter Statistiken abgeleitet.

Für $n=3$ ist insbesondere die Zufallszahl

$$y = mid \{ u_{1},u_{2},u_{3} \}$$

durch die einfache Dichtefunktion

$$p(y) = 6 y(1-y)$$

darstellbar (siehe Kap.3.1).

Beispiel. Ein anderes wichtiges Beispiel hatten wir bereits bei der Herleitung der $\chi^{2}$- Verteilung gefunden. Wenn $x_{1},x_{2},...,x_{n}$ unabhängig normalverteilt sind, mit Mittelwerten $<x_{\nu}> = 0, \; \nu=1,2,...,n,$ und Varianzen $\sigma_{\nu}=1$,

$$p(x_{1},x_{2},...,x_{n}) = \frac{1}{(2\pi)^{n/2}} e^{-\frac{1}{2} \sum_{\nu=1}^{n} x_{\nu}^{2}},$$

dann ist die Summe der Quadrate

$$y = \sum_{\nu=1}^{n} x_{\nu}^{2}$$

$\chi^{2}$- verteilt mit der Dichtefunktion

$$p(y) = \frac{1}{2^{n/2}\Gamma(\frac{n}{2})} y^{n/2-1} e^{-y/2}.$$

$n$ nennt man die Anzahl der Freiheitsgrade.

Composition- Methode
Wir kommen noch einmal zurück auf die Formel

$$p(x_{1},x_{2}) = p_{1}(x_{1}) p_{2}(x_{2}\vert x_{1}).$$ (26)

Diese Formel gestattet in vielen Fällen eine bequeme Simulation der Randverteilung

$$p_{2}(x_{2}) = \int_{-\infty}^{\infty} dx_{1} p_{1}(x_{1}) p_{2}(x_{2}\vert x_{1}),$$ (27)

sofern Algorithmen zur Simulation von $p_{1}(x_{1})$ und $p_{2}(x_{2}\vert x_{1})$ bereits existieren. Dieses Verfahren ist in der Literatur auch unter dem englischen Ausdruck Composition- Methode bekannt.

Die Composition Methode kann in einfacher Weise auf gemischte Verteilungen erweitert werden, wobei dann allerdings $x_{1} \equiv k_{1}$ eine diskrete Veränderliche wird. Wir schreiben

$$p(k_{1},x_{2}) = P_{1}(k_{1}) p_{2}(x_{2}\vert k_{1})$$ (28)

und

$$p_{2}(x_{2}) = \sum_{k_{1}=0}^{\infty} P_{1}(k_{1}) p_{2}(x_{2}\vert k_{1}) .$$ (29)

Der Algorithmus lautet also: Bilde die integralen Verteilungen

$$P_{1}(\xi_{1}) = \int_{-\infty}^{\xi_{1}} dx_{1} p_{1}(x_{1}),$$ (30)

$$P_{2}(\xi_{2}\vert\xi_{1}) = \int_{-\infty}^{\xi_{2}} dx_{2} p_{2}(x_{2}\vert\xi_{1}),$$

bzw

$$P_{1}(\kappa_{1}) = \sum_{k_{1}=0}^{\kappa_{1}} P_{1}(k_{1}),$$ (31)

$$P_{2}(\xi_{2}\vert\kappa_{1}) = \int_{-\infty}^{\xi_{2}} dx_{2} p_{2}(x_{2}\vert\kappa_{1}),$$

und löse die Bestimmungsgleichungen

$$P_{1}(\xi_{1}) = u_{1},$$ (32)

$$P_{2}(\xi_{2}\vert\xi_{1}) = u_{2},$$

bzw

$$P_{1}(\kappa_{1}-1) < u_{1} \leq P_{1}(\kappa_{1}),$$ (33)

$$P_{2}(\xi_{2}\vert\kappa_{1}) = u_{2},$$

wobei wie immer $u_{1}, u_{2} \in [0,1]$ zwei gleichverteilte Zufallszahlen sind. $\xi_{2}$ ist dann eine Veränderliche der Dichtefunktion $p_{2}(x_{2})$.

Beispiel. Das wichtigste Beispiel in diesem Rahmen ist eine Simulation der Normalverteilung. Wir beginnen mit der zweidimensionalen Normalverteilung

$$p(x_{1},x_{2}) = \frac{1}{2\pi} e^{-\frac{1}{2}(x_{1}^{2}+x_{2}^{2})}$$

und transformieren gemäß

$$x_{1} = r cos \phi$$ (34)

$$x_{2} = r sin \phi$$

mit $y_{1} \equiv r$ und $y_{2} \equiv \phi$. Dieses führt auf die Dichtefunktion der Veränderlichen $r$ und $\phi$:

$$p(r,\phi ) = \frac{1}{2\pi} r e^{-\frac{1}{2} r^{2}}.$$

Wir faktorisieren die Dichtefunktion in der Form

$$p(r,\phi ) = p_{1}(r) p_{2}(\phi \vert r),$$

und identifizieren

$$p_{1}(r) = \int_{0}^{2\pi} d\phi p(r,\phi ) = r e^{-\frac{1}{2}r^{2}},$$

$$p_{2}(\phi \vert r) = \frac{p(r,\phi )}{p_{1}(r)} = \frac{1}{2\pi}.$$

Für $p_{1}(r)$ und $p_{2}(\phi \vert r)$ haben wir aber bereits Algorithmen zur Simulation. $p_{1}(r)$ kann mit der inversen Transformationsmethode (siehe Beispiel am Anfang dieses Kapitels) gelöst werden und $p_{2}(\phi \vert r)$ ist eine Gleichverteilung in $[0,2\pi ]$. Damit erhalten wir den Algorithmus

$$r = \sqrt{-2 ln(u_{1})},$$

$$\phi = 2\pi u_{2},$$

mit $u_{1}, u_{2} \in [0,1]$. Die Rücktransformation gemäß (43) liefert zwei normalverteilte unabhängige Veränderliche.

Beispiel. Gesucht sei eine Simulation der Dichtefunktion

$$p_{2}(x_{2}) = n \int_{1}^{\infty} dx_{1} x_{1}^{-n} e^{-x_{1} x_{2}}, \; \; \; \; x_{2} > 0.$$

Wir gehen von der zweidimensionalen Dichtefunktion

$$p(x_{1},x_{2}) = n x_{1}^{-n} e^{-x_{1} x_{2}}, \; \; \; x_{1}>1, \; x_{2}>0,$$

aus und berechnen

$$p_{1}(x_{1}) = n x_{1}^{-(n+1)}$$

$$p_{2}(x_{2}\vert x_{1}) = \frac{p(x_{1},x_{2})}{p_{1}(x_{1})} = x_{1} e^{-x_{1} x_{2}}.$$

Damit ergibt sich für die integralen Verteilungsfunktionen

$$P_{1}(\xi_{1}) = 1- \xi_{1}^{-n},$$

$$P_{2}(\xi_{2}\vert\xi_{1}) = 1 - e^{-\xi_{1} \xi_{2}}.$$

Mit Hilfe von zwei gleichverteilten Zufallszahlen $u_{1}, u_{2} \in [0,1]$ erhalten wir die Lösung

$$\xi_{1} = u_{1}^{-1/n},$$

$$\xi_{2} = -\frac{1}{\xi_{1}} ln(u_{2}).$$

$\xi_{2}$ ist dann die gesuchte Veränderliche der Dichtefunktion $p_{2}(x_{2})$.

Beispiel. Die Composition Methode kann besonders vorteilhaft für die Simulation von Dichtefunktionen verwendet werden, die als Reihenentwicklung darstellbar sind,

$$p_{2}(x) = \sum_{k=0}^{\infty} (k+1) P_{k} x^{k}, \; \; \; 0 \leq x \leq 1.$$

Die Normierung verlangt, daß

$$\sum_{k=0}^{\infty} P_{k} = 1.$$

Wir setzen $P(k) \equiv P_{k}$ und

$$p(k,x) = (k+1) P(k) x^{k}, \; \; \; 0 \leq x \leq 1,\; \; k \geq 0.$$

Damit erhalten wir der Reihe nach

$$P_{1}(k) = P(k) \int_{0}^{1}dx (k+1) x^{k} = P(k) \equiv P_{k},$$

$$p_{2}(x\vert k) = \frac{p(k,x)}{P_{1}(k)} = (k+1) x^{k},$$

und

$$P_{1}(\kappa) = \sum_{k=0}^{\kappa} P_{k},$$

$$P_{2}(\xi\vert\kappa) = (\kappa+1) \int_{0}^{\xi} dx x^{k} = \xi^{\kappa+1} .$$

Wir erzeugen zwei Zufallszahlen $u_{1}, u_{2} \in [0,1]$ und erhalten $\kappa$ aus der Ungleichung

$$\sum_{k=0}^{\kappa-1} P_{k} < u_{1} \leq \sum_{k=0}^{\kappa} P_{k},$$

sowie $\xi$ aus der Gleichung

$$\xi = u_{2}^{1/(1+\kappa)}.$$

$\xi$ ist dann eine Veränderliche der Dichtefunktion $p_{2}(x)$.

Das Verfahren kann verallgemeinert werden für Summen beliebiger Dichtefunktionen $q_{k}(x)$,

$$p_{2}(x) = \sum_{k=0}^{n} P_{k} q_{k}(x),$$

mit

$$P_{0} = 0,$$

$$\sum_{k=0}^{n} P_{k} = 1,$$

$$\int_{-\infty}^{\infty} dx q_{k}(x) = 1, \; \; \; k=1,2,...,n.$$

Die zugrundeliegende zweidimensionale Dichtefunktion ist also

$$p(k,x) = P(k) q_{k}(x), \; \; \; 0 \leq k \leq n,$$

mit $P(k) \equiv P_{k}$. Die Randverteilungen sind also

$$P_{1}(k) = P(k) \equiv P_{k},$$

$$p_{2}(x\vert k) = q_{k}(x).$$

Mit Hilfe der integralen Verteilungsfunktionen

$$P_{1}(\kappa) = \sum_{k=0}^{\kappa} P_{k}, \; \; \; \kappa = 0,1,2,...,n,$$

$$P_{2}(\xi\vert\kappa) = \int_{-\infty}^{\xi} dx q_{\kappa}(x) \equiv Q_{\kappa }(\xi)$$ (35)

lautet die Formulierung des Verfahrens: Erzeuge zwei Zufallszahlen $u_{1}, u_{2} \in [0,1]$ und bestimme $\kappa$ aus

$$\sum_{k=0}^{\kappa -1} P_{k} < u_{1} \leq \sum_{k=0}^{\kappa} P_{k}.$$

Löse die Gleichung

$$P_{2}(\xi \vert \kappa ) = u_{2}$$

nach $\xi$ auf. $\xi$ ist dann eine Veränderliche der Dichtefunktion $p_{2}(x)$.

Verwerfungsmethoden Bei der Verwerfungsmethode stellen wir die gesuchte Dichtefunktion $q(x)$ in der Form

$$q(x) = C h(x) G(x)$$ (36)

dar. Dazu schlagen wir folgenden Weg ein. Wir beginnen mit einer zweidimensionalen Dichtefunktion in der üblichen Darstellung

$$p(x_{1},x_{2}) = p_{1}(x_{1}) p_{2}(x_{2}\vert x_{1}).$$

Das Integral über $x_{2}$ in den Grenzen von $-\infty$ bis $x_{1}$ kann geschrieben werden als

$$\int_{-\infty}^{x_{1}} dx_{2} p(x_{1},x_{2}) = p_{1}(x_{1}) ... ...{2}(x_{2}\vert x_{1}) = p_{1}(x_{1}) P_{2}(x_{1}\vert x_{1}).$$ (37)

Die rechts stehende Wahrscheinlichkeit $P_{2}$ lautet ausführlich geschrieben:

$$P_{2}(x_{1}\vert x_{1}) \equiv P_{2}(x_{2} \leq x_{1}\vert x_{1}).$$ (38)

Die linke Seite von (37) ist nun eine Dichtefunktion der Veränderlichen $x_{1}$, die aber im allgemeinen nicht auf 1 normiert ist. Wir setzen daher

$$q(x_{1}) = C \int_{-\infty}^{x_{1}} dx_{2} p(x_{1},x_{2})$$

und erhalten

$$q(x_{1}) = C p_{1}(x_{1}) P_{2}(x_{1}\vert x_{1}).$$ (39)

Ein Vergleich mit (36) zeigt, daß

$$h(x) \equiv p_{1}(x),$$ (40)

$$G(x) \equiv P_{2}(x\vert x).$$ (41)

Aus dieser Identifikation ergibt sich, daß für $h(x)$, $G(x)$ und $C$ die folgenden Bedingungen erfüllt sein müssen:

$$\int_{-\infty}^{\infty} dx h(x) = 1,$$

$$max \{ G(x) \} = 1,$$ (42)

$$C > 1.$$

Wir hatten die integrale bedingte Verteilungsfunktion $P_{2}(\xi_{2}\vert\xi_{1})$ definiert als Wahrscheinlichkeit dafür, daß die Veränderliche $x_{2} \leq \xi_{2}$, aber unter der Vorraussetzung, daß die Veränderliche $x_{1}$ einen festen Wert $x_{1}=\xi_{1}$ annimmt. Die Funktion $P_{2}(\xi_{1}\vert\xi_{1})$ gibt dann die Wahrscheinlichkeit an, daß die Veränderliche $x_{2}$ kleiner als der Wert $\xi_{2} = \xi_{1}$ sein soll, unter der Vorraussetzung, daß der Wert für $x_{1}=\xi_{1}$ bereits festliegt. Offensichtlich bedeutet dieses, daß $\xi_{1}$ nur dann eine Veränderliche der Dichtefunktion $p(x_{1},x_{2})$ sein kann, wenn die Bedingung

$$u_{2} \leq P_{2}(\xi_{1}\vert\xi_{1})$$ (43)

mit einer zweiten Zufallszahl $u_{2}\in[0,1]$ erfüllt ist.

Zusammengefaßt erhalten wir folgenden Algorithmus: Erzeuge zwei Zufallszahlen $u_{1}, u_{2} \in [0,1]$, löse die Gleichung

$$H(\xi) = \int_{-\infty}^{\xi} dx h(x) = u_{1}$$ (44)

nach $\xi$ auf und prüfe, ob

$$u_{2} \leq G(\xi).$$ (45)

Wenn diese Bedingung nicht erfüllt ist, verwerfen wir $\xi$ und beginnen von vorne.

Eine wichtige Größe bei dieser Methode ist die Zahl $C$. Zunächst ist klar, daß $C > 1$ ist (siehe (42)). Wie man sich leicht klarmachen kann, ist $C$ die mittlere Anzahl der Versuche, bis zu der die Ungleichung (45) erfüllt ist. Die Größe $1/C$ nennt man daher auch die Erfolgswahrscheinlichkeit oder Effizienz des Verfahrens. Je größer die Effizienz, desto schneller ist das Verwerfungsverfahren.

Der optimale Wert für $C$ wäre eigentlich $C=1$. Das ist natürlich nicht möglich. $C$ erreicht dagegen umso näher den Wert 1, je größer die Funktion $G(x)$ ist. Man muß also immer versuchen, die Funktion $G(x)$ maximal zu wählen, aber unter Einhaltung der Nebenbedingung $max \{ G(x) \} \equiv G_{max} = 1.$

Beispiel. Ein Beispiel, das auf die allgemeine Verwerfungsmethode zurückgeht, ist das folgende. Gesucht ist ein Simulationsalgorithmus der Dichtefunktion

$$q(x) = \frac{1}{\Gamma(\alpha)} x^{\alpha-1} e^{-x}, \; \; \; 0 < \alpha < 1, \; x \geq 0.$$

Wir setzen

$$h(x) = \{ \begin{array}{ll} \frac{\alpha e}{\alpha +e} x^{\... ... e}{\alpha +e} e^{-x}, \; \; \; & 1 < x < \infty, \end{array}$$

$$G(x) = \{ \begin{array}{ll} e^{-x}, \; \; \; & 0 \leq x \leq 1, \\ x^{\alpha-1}, \; \; \; & 1 < x < \infty, \end{array}$$

$$C = \frac{1}{\Gamma(\alpha)} \frac{\alpha + e}{\alpha e}.$$

Der Grund für die gewählte Aufteilung der Funktion $q(x)$ in die Faktoren $h(x)$, $G(x)$ und $C$ folgt aus unseren obigen Überlegungen. Die Wahl der Funktion $h(x)$ ergibt sich aus der Normierung

$$\int_{0}^{\infty} dx h(x) = \frac{\alpha e}{\alpha +e} \int... ... e^{-x} = \frac{e}{\alpha +e} + \frac{\alpha}{\alpha +e} = 1.$$

Für $G(x)$ ist die Bedingung

$$max \{ G(x) \} = G_{max} = 1$$

offensichtlich für $x=0$ und $x=1$ erfüllt. Der Ausdruck für $C$ folgt dann aus

$$C = \frac{q(x)}{h(x) G(x)}.$$

Die integralen Verteilungsfunktionen sind

$$H(\xi) = \{ \begin{array}{ll} \frac{e}{\alpha +e} \xi^{\alph... ...{\alpha +e} e^{1-\xi}, \; \; \; & 1 < \xi < \infty, \end{array}$$

$$G(\xi) = \{ \begin{array}{ll} e^{-\xi}, \; \; \; & 0 \leq \x... ...1, \\ \xi^{\alpha -1}, \; \; \; & 1 < \xi < \infty. \end{array}$$

Auf diesen Ergebnissen aufbauend, können wir folgenden Algorithmus entwerfen: Erzeuge zwei Zufallszahlen $u_{1}, u_{2} \in [0,1]$. Wenn

$$\xi = \left( \frac{\alpha +e}{e} u_{1} \right)^{1/\alpha} \leq 1,$$

dann prüfe, ob

$$u_{2} \leq e^{-\xi}.$$ (46)

Falls dagegen $\xi > 1$, dann ersetze $\xi$ durch

$$\xi = 1 - ln \left( \frac{\alpha +e}{\alpha} (1 -u_{1}) \right)$$

und prüfe, ob

$$u_{2} \leq \xi^{\alpha -1}.$$ (47)

Wenn die Bedingung (46) bzw (47) nicht erfüllt ist, beginne von vorne. Dieses Beispiel finden Sie in dem folgenden Applet unter der Bezeichnung Gamma.Unter dem weiteren Menupunkt BreitWigner finden Sie eine Simulation von BreitWigner- Verteilungen mit einem speziellen Verfahren, das wir im folgenden diskutieren werden.

Das von Neumann'sche Verfahren. Der Mathematiker von Neumann hat das Verwerfungsverfahren in der folgenden Form vorgeschlagen. Sei $q(x)$ eine im Intervall $a \leq x \leq b$ definierte Dichtefunktion, für die wir einen Simulationsalgorithmus suchen. Wir setzen für $h(x)$ die spezielle Gleichverteilung

$$h(x) = \frac{1}{b-a}, \; \; \; a \leq x \leq b.$$ (48)

Dann können wir $G(x)$ durch die vorgegebene Dichtefunktion $q(x)$ ausdrücken. Wegen

$$q(x) = C h(x) G(x)$$ (49)

erhalten wir

$$G(x) = \frac{b-a}{C} q(x).$$ (50)

Die Zahl $C$ muß hierbei so gewählt werden, daß $max \{ G(x) \} = 1$ für $x \in [a,b]$, d.h.

$$\frac{b-a}{C} q_{max} = 1,$$

oder

$$C = (b-a) q_{max}.$$ (51)

Hierbei ist $q_{max}$ das Maximum der Dichtefunktion $q(x)$. Wir können dann schreiben:

$$G(x) = \frac{q(x)}{q_{max}}.$$ (52)

Wir formulieren das Verfahren: Erzeuge zwei Zufallszahlen $u_{1}, u_{2} \in [0,1]$ und berechne

$$\xi = a + (b-a) u_{1}.$$ (53)

Wenn dann

$$u_{2} \leq \frac{q(\xi)}{q_{max}},$$ (54)

ist $\xi$ eine Veränderliche der Dichtefunktion $q(x)$, sonst beginne von vorne.

Das besondere an dem Verfahren von v. Neumann ist offensichtlich, daß keine inverse Funktion berechnet werden muß. Der Algorithmus ist auch für die kompliziertesten Dichtefunktionen $q(x)$ universell anwendbar. Der Nachteil des Verfahrens ist, daß für große Intervalle $[a,b]$ die Erfolgswahrscheinlichkeit $1/C$ klein wird.

Beispiel. Wir diskutieren das v. Neumannsche Verfahren an dem einfachen Beispiel

$$q(x) = \frac{2}{\pi R^{2}} \sqrt{R^{2} - x^{2}}, \; \; \; -R \leq x \leq R.$$

Das Maximum dieser Funktion wird offensichtlich für $x=0$ angenommen,

$$q_{max} = \frac{2}{\pi R}.$$

Mit $a=-R$ und $b=+R$ erhalten wir

$$h(x) = \frac{1}{2 R},$$

$$G(x) = \frac{1}{R} \sqrt{R^{2}-x^{2}}.$$

Das Verfahren lautet also: Erzeuge $u_{1}, u_{2} \in [0,1]$, berechne

$$\xi = -R + 2 R u_{1},$$

und prüfe, ob

$$u_{2} \leq \frac{1}{R} \sqrt{R^{2}-\xi^{2}}.$$

Wenn diese Bedingung erfüllt ist, ist $\xi$ eine Veränderliche der Dichtefunktion $q(x)$. Wenn nicht, verwerfe $\xi$ und beginne von vorne.

Beispiel Die Verwerfungsmethode von Neumann wird im folgenden Applet an einer pathologischen Dichtefunktionen erläutert. Wir hatten oben in diesem Kapitel bereits eine Simulation der Cauchy- Verteilung

$$p(x) = \frac{1}{\pi} \frac{1}{(1 + x^{2})}$$

diskutiert. Diese Vertelung ist im Bereich von $-\infty$ bis $+\infty$ auf eins normiert,

$$\int_{-\infty}^{+\infty} dx \; p(x) = 1.$$

Das merkwürdige ist allerdings, daß man keine höheren Momente berechnen kann, d.h.

$$\int_{-\infty}^{+\infty} dx \; x^{k} p(x) \rightarrow \infty \; \; \; \; \; fuer \; k > 0.$$

Man benötigt also einen Konvergenz erzeugenden Multiplikator. Die Chauchy- Verteilung kommt hauptsächlichst in der Teilchenphysik zur Beschreibung von Resonanzen zum Einsatz. In diesem Fall ist der Multiplikator bereits von der Physik her vorgegeben, und zwar der Phasenraum- Faktor, der für die Energie- und Impulserhaltung sorgt. Man nennt die Verteilung dann auch die Breit- Wigner- Verteilung. Mit der Transformation

$$x = \frac{M - M_{R}}{\Gamma/2}$$

erhält man

$$p(M) = \frac{1}{\pi} \; \frac{\Gamma/2}{(M - M_{R})^{2} + (\Gamma/2)^{2}}$$

Der Phasenraumfaktor kann vereinfacht geschrieben werden als

$$X(M) = \sqrt{(M- M_{0})(M_{max}-M)} \; \; \; \; \; fuer \; M_{0} < M < M_{max}$$

$$= 0 \; \; \; \; \; fuer \; M < M_{0} \; \; \; \; \; und \; M > M_{max}.$$

und die Gesamtverteilung ist

$$p(M) = \frac{1}{\pi} \; \frac{(\Gamma/2) X(M)}{(M - M_{R})^{2} + (\Gamma/2)^{2}}$$

In dem Applet studieren wir die Summe dreier Resonanzen. Sie werden bemerken, daß das Verwerfungsverfahren kein besonders effektiver Algorithmus mit einer Effizienz von nur wenigen Prozenten ist und daher auch nur im Notfall eingesetzt werden sollte.

In diesem und in allen folgenden Applets dieses Kapitels ist die Bedienung immer gleich. Im Menu Generator wählen wir zunächst einen bestimmten Generator aus. Es erscheint ein weiteres kleines Fenster, in dem Sie die Parameter der Verteilung ändern können. Auf jeden Fall müssen Sie dieses Fenster mit OK quittieren. Im Menu Plots können dann mehrere Fenster geöffnet werden, wie z.B. die Visualisierung der erzeugten Verteilung und die exakte Kurve der mathematischen Funktion der Verteilung. Bei Anklicken des Source- Buttons wird ein Fenster mit dem Source- Code des Algorithmus aufgeklappt. Bei Start beginnen wir die Simulation, bei Stop wird sie abgebrochen.

Die Methode von Forsythe
Eine etwas komplizierte Methode ist von Forsythe entwickelt worden. Wir behandeln dieses Verfahren etwas ausführlicher, weil es insbeondere gestattet, für bestimmte Lösungen des Sattelpunktverfahrens (siehe Kap.3.4) Simulationsprogramme zu entwickeln. Die Methode ist anwendbar für irgendeine Dichtefunktion in der Form

$$p(x) = C e^{-h(x)},$$ (55)

wobei $h(x)$ eine stetig anwachsende Funktion von $x$ im Bereich $[0,\infty]$ ist. Wir unterteilen die positive reelle Achse in $K$ Intervalle zwischen den Stützstellen $\xi^{(0)}=0,\xi^{(1)},\xi^{(2)},...,\xi^{(K)}$ und berechnen

$$P(\xi^{(k)}) = \int_{0}^{\xi^{(k)}} dx p(x), \; \; \; k=1,2,...,K$$ (56)

und die Intervallängen

$$d_{k} = \xi^{(k)} - \xi^{(k-1)}.$$ (57)

Wählt man die Stützstellen $\xi^{(k)}$ nun so, daß stets

$$g_{k}(x) = h(x) - h(\xi^{(k-1)}) \leq 1 \; \; \; \forall \; x \in [\xi^{(k-1)},\xi^{(k)}],$$ (58)

so führt der folgende Algorithmus auf ein Simulationsprogramm für die Dichtefunktion $p(x)$: Erzeuge zwei Zufallszahlen $u_{-1},u_{0} \in [0,1]$. Bestimme $\kappa$ so, daß

$$P(\xi^{(\kappa -1)}) < u_{-1} \leq P(\xi^{(\kappa )}),$$ (59)

und wähle einen Wert der Veränderlichen $x$ im $\kappa$-ten Intervall gemäß

$$\xi = \xi^{(\kappa -1)} + u_{0} d_{\kappa }$$ (60)

Erzeuge solange weitere Zufallszahlen $u_{\nu} \in [0,1]$, $\nu=1,2,...,$ bis die Bedingung

$$u_{n-1} < g_{\kappa }(\xi ) \leq u_{n}$$ (61)

erfüllt ist. Wenn $n$ ungerade ist, akzeptiere $\xi$ als Veränderliche der Dichtefunktion $p(x)$. Wenn dagegen $n$ gerade ist, verwerfe $\xi$ und beginne von vorne.

Die Forsythe- Methode erscheint auf den ersten Blick ziemlich umständlich und zeitraubend. Der große Vorteil liegt jedoch darin, daß die Exponentialfunktion nur in der Initialisierungsphase berechnet werden muß, ähnlich wie in der Tabellen- Darstellung. Gegenüber der Tabellen- Methode liefert die Forsythe- Methode jedoch eine exakte Simulation der zufälligen Veränderlichen

Beispiel. Als erstes Beispiel behandeln wir eine Simulation der normalisierten Exponentialfunktion

$$p(x) = e^{-x},$$

d.h.

$$h(x) = x.$$

Bei der folgenden Wahl der Stützstellen,

$$\xi^{(k)} = k, \; \; \; k=0,1,2,...$$

gilt

$$d_{k} = 1, \; \; \; k=1,2,3,....$$

$$g_{k}(x) = x-(k-1) \; \leq 1 \; \forall x \in [k-1,k]$$

und

$$P(k) = 1 - e^{-k}.$$

Ersichtlich genügen etwa 10 Stützstellen $\xi^{(k)}, k=1,2,...,10$. Das Verfahren lautet also im einzelnen: 1. Berechne die Größen

$$P(0) = 0,$$

$$P(k) = 1 - e^{-k}, \; \; \; k=1,2,...,9,$$

$$P(10) = 1$$

in einem externen Initialisierungsprogram. 2. Erzeuge zwei Zufallszahlen $u_{-1},u_{0} \in [0,1]$, bestimme $\kappa$ gemäß

$$P(\kappa-1) < u_{-1} \leq P(\kappa)$$

und $\xi$ gemäß

$$\xi = (\kappa -1) + u_{0}.$$

Erzeuge solange weitere Zufallszahlen $u_{1}, u_{2},...,u_{n} \in [0,1],$ bis die Bedingung

$$u_{n-1} < \xi - (\kappa -1) \leq u_{n}$$

erfüllt ist. Falls $n$ ungerade ist, akzeptiere $\xi$ als zufällige Veränderliche der Dichtefunktion $p(x)$, sonst beginne von vorne bei Punkt 2.

Beispiel. Als zweites Beispiel behandeln wir die Dichtefunktion

$$p(x) = \sqrt{\frac{2}{\pi}} e^{-x^{2}/2}, \; \; \; x \geq 0,$$

d.h. die normalisierte rechte Hälfte der Standard- Normalverteilung. Jetzt ist

$$h(x) = \frac{1}{2} x^{2}$$

und die folgende Wahl der Stützstellen,

$$\xi^{(0)} = 0,$$

$$\xi^{(1)} = 1,$$

$$\xi^{(k)} = \sqrt{2 k -1}, \; \; \; k \geq 2$$

erfüllt die Bedingungen der Forsythe- Methode. Wir erhalten

$$d_{1} = 1,$$

$$d_{k} = \sqrt{2k-1} - \sqrt{2k-3}, \; \; \; k \geq 2$$

und

$$g_{1}(x) = \frac{1}{2} x^{2}, \; \; \; \forall x \in [0,1],$$

$$g_{k}(x) = \frac{1}{2} [x^{2} - 2k +3], \; \; \; \forall x \in [\sqrt{2k-3}, \sqrt{2k-1} ], \; \; \; k \geq 2.$$ (62)

Die Tabellenwerte

$$P(\xi^{(k)}) = \sqrt{\frac{2}{\pi}} \int_{0}^{\sqrt{2k-1}} dx e^{-x^{2}/2}$$

müssen wiederum extern berechnet werden. Ersichtlich genügen in diesem Beispiel ungefähr 5 Stützstellen, sodaß der Algorithmus folgendermaßen formuliert werden kann:
1. Setze $P(\xi^{(0)}) = 0, P(\xi^{(5)}) = 1, \xi^{(0)} = 0, \xi^{(1)} = 1, d_{1} = 1$ und berechne in einem externen Programm

$$\xi^{(k)} = \sqrt{2k-1}, \; \; \; k \geq 2,$$

$$d_{k} = \xi^{(k)} - \xi^{(k-1)}, \; \; \; k \geq 2,$$

$$P(\xi^{(k)}) = \sqrt{\frac{2}{\pi}} \int_{0}^{\sqrt{2k-1}} dx e^{-x^{2}/2}, \; \; \; k=1,2,3,4.$$

2. Erzeuge zwei Zufallszahlen $u_{-1},u_{0} \in [0,1]$, bestimme $\kappa$ gemäß

$$P(\xi^{(\kappa-1)}) < u_{-1} \leq P(\xi^{(\kappa)})$$

und $\xi$ gemäß

$$\xi = \xi^{(\kappa-1)} + u_{0} d_{\kappa}.$$

Erzeuge solange weitere Zufallszahlen $u_{\nu} \in [0,1], \nu = 1,2,3,...$, bis die Bedingung

$$u_{n-1} < \frac{1}{2} [\xi^{2} - 2 \kappa + 3] \leq u_{n}$$

erfüllt ist. Falls $n$ ungerade ist, akzeptiere $\xi$ als Veränderliche der Dichtefunktion $p(x)$, sonst beginne von vorne bei Punkt 2.