Kontinuierliche Verteilungen I

Übersicht
In den nächsten drei Abschnitten des vorliegenden Kapitels werden wir noch einmal eine Zusammenstellung der wichtigsten Algorithmen zur Simulation einfacher zufälliger Veränderlicher geben. Die hier getroffene Auswahl wird sicherlich nicht erschöpfend sein. Sie soll vielmehr dem Leser einen Eindruck davon vermitteln, wie vielfältig die Ideen und Methoden sind, und damit den Leser auch zu eigener Entwicklungsarbeit anregen. Den Source Code zu allen Algorithmen finden Sie in unsreren Java- Applets.

Zur klareren Bezeichnung der Verteilungen führen wir noch einige mnemotechnische Abkürzungen ein. So werden wir z.B. eine gleichverteilte Veränderliche aus dem Zahlenintervall von $a$ bis $b$ mit $u \in U[a,b]$ kennzeichnen. Eine normalverteilte Veränderliche entsprechend mit $x \in N[a,\sigma ]$. Die Dichtefunktion bzw Wahrscheinlichkeitsverteilung wird dabei mit einem mnemotechnisch einfachen Symbol ($U$ = uniform, $N$ = normal) gekennzeichnet, die Parameter der Verteilung in eckigen Klammern angefügt. Eine Zusammenstellung der wichtigsten kontinuierlichen Verteilungen und ihrer Symbole ist in Tabelle 1 gelistet. Für die folgenden Erörterungen verweisen wir auf das erste Applet dieses Kapitels.

Tabelle: Zusammenstellung der wichtigsten Verteilungen kontnuierlicher Veränderlicher und ihrer mnemotechnischen Symbole.

Dichtefunktion		Symbol	Bezeichnung der Verteilung
$\frac{1}{b-a}$	$a \leq x \leq b$	$U[a,b]$	Gleichverteilung
$\frac{1}{\sqrt{2\pi}\sigma} e^{-(x-a)^{2}/2\sigma^{2}}$		$N[a,\sigma ]$	Normalverteilung
$\frac{1}{\sqrt{2\pi} \sigma x} e^{-(ln x -a)^{2}/2\sigma^{2}}$	$0 \leq x$	$LN[a,\sigma ]$	Lognormalverteilung
$\frac{1}{\tau} e^{-x/\tau}$	$\tau > 0, x \geq 0$	$E[\tau]$	Exponentialverteilung
$\frac{x^{\alpha -1} e^{-x/\beta}}{\beta^{\alpha} \Gamma(\alpha)}$	$\alpha > 0, \beta > 0, 0 \leq x$	$G[\alpha, \beta]$	Gammaverteilung
$\frac{x^{m-1} e^{-x/\beta}}{\beta^{m} \Gamma(m)}$	$m > 0, \beta > 0, 0 \leq x$	$Er[m,\beta]$	Erlang- Verteilung

Um die HTML- Seite nicht zu lang werden zu lassen, haben wir die kontnuierlichen Verteilungen in zwei Kapitel unterteilt. Im vorliegenden Kapitel werden wir die Normalverteilung, die Exponentialverteilung und die Gammaverteilung diskutieren, im nächsten Kapitel werden dann alle übrigen Verteilungen behandelt.

Normalverteilung
Wir beginnen mit der Normalverteilung $N[a,\sigma ]$ in der allgemeinen Form

$$p(y) = \frac{1}{\sqrt{2\pi}\sigma } e^{-(y-a)^{2}/2\sigma^{2}}.$$ (1)

Es genügt ersichtlich ein Algorithmus für die sogenannte Standard- Normalverteilung $N[0,1]$,

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

da die Transformation $y = a + \sigma x$ eine zufällige Veränderliche der allgemeinen Normalverteilung ergibt. Die erzeugende Funktion der Normalverteilung ist

$$M_{y}(v) = e^{av+\frac{1}{2}\sigma^{2} v^{2}},$$ (3)

sowie die logarithmisch erzeugende Funktion ist

$$H_{y}(v) = av + \frac{1}{2} \sigma^{2} v^{2},$$ (4)

mit den zwei Momenten

$$<y> = a,$$

$$<(y-a)^{2}> = \sigma^{2}.$$

Wir hatten die Normalverteilung und ihre Eigenschaften im Verlaufe dieses Kursus schon häufig diskutiert. Der breite Anwendungsbereich wird uns auch noch in den folgenden Kapiteln begegnen.

Die wichtigste Eigenschaft der Normalverteilung folgt aus dem schon mehrfach erwähnten zentralen Grenzwertsatz. Aus diesem Satz folgt, daß, wenn $x_{1}, x_{2},...,x_{n}$ irgendwelche unabhängigen zufälligen Veränderlichen sind, die Summe $x=x_{1}+x_{2}+ .... + x_{n}$ im Grenzfall $n \to \infty$ gegen eine Normalverteilung konvergiert. Auch für endliche Werte von $n$ ist diese Behauptung schon sehr gut erfüllt. Aus dieser Aussage folgt insbesondere, daß Meßfehler, die sich ja im allgemeinen aus einer Vielzahl von unabhängigen statistischen und systematischen Einzelfehlern additiv zusammenstzen, normalverteilt sind. Aus diesem Grunde gibt man als Maß für einen Meßfehler die einfache Varianz $\sigma$ oder den zweifachen Wert $2 \sigma$ der Varianz an.

Methode der Randverteilungen.
Den Algorithmus mit Hilfe der Methode der Randverteilungen hatten wir bereits ausführlich im vorigen Abschnitt diskutiert. Es genügt daher die folgende Zusammenfassung: Erzeuge $u_{1},u_{2} \in U[0,1]$, dann sind

$$\xi_{1} = \sqrt{-2 ln(u_{1})} cos(2\pi u_{2}),$$ (5)

$$\xi_{2} = \sqrt{-2 ln(u_{1})} sin(2\pi u_{2}),$$ (6)

zwei unabhängige normalverteilte Veränderliche aus $N[0,1]$. Die entsprechende Methode heißt haben wir in dem Applet mit JNormal1 bezeichnet.

Verwerfungsmethode. Ein sehr einfacher Algorithmus folgt aus der Verwerfungsmethode. Sei

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

die auf 1 normierte rechte Hälfte einer Standard- Normalverteilung. Da die Normalverteilung $N[0,1]$ symmetrisch um Null ist, genügt es, eine gemäß (7) erzeugte Veränderliche mit einem zufälligen Vorzeichen zu versehen. Wir setzen

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

$$G(x) = e^{-(x-1)^{2}/2},$$ (8)

$$C = \sqrt{\frac{2 e}{\pi}}.$$

Die Funktion $G(x)$ erreicht ihr Maximum $G_{max}=1$ an der Stelle $x=1$, die Dichtefunktion $h(x)$ ist für $x \geq 0$ auf 1 normiert und $C$ ist $\approx 1.32$, also größer als 1. Einsetzen in den Ausdruck

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

ergibt die Dichtefunktion (7). Alle Bedingungen der Verwerfungsmethode sind also erfüllt. Die Erfolgswahrscheinlichkeit ist mit $1/C \approx 0.76$ verhältnismäßig gut. Der Algorithmus lautet dann wie folgt: Erzeuge $u_{1},u_{2},u_{3} \in U[0,1]$ und setze

$$\xi = - ln(u_{1}).$$ (9)

Falls

$$u_{2} \leq e^{-(\xi -1)^{2}/2} ,$$ (10)

akzeptiere $\xi$ als Veränderliche, sonst beginne von vorne. Wenn

$$u_{3} \leq 0.5,$$ (11)

ersetze $\xi$ durch $-\xi$. $\xi$ ist dann eine Veränderliche aus $N[0,1]$.

Die Bedingung (10) kann ersetzt werden durch

$$- ln(u_{2}) \geq \frac{1}{2} (\xi -1)^{2}.$$

Mit dieser Feststellung kann man das Verfahren auch in der folgenden Form formulieren: Erzeuge $v_{1},v_{2} \in E[1]$ und $u \in U[0,1]$ und setze

$$\xi = v_{1}.$$ (12)

Wenn dann

$$v_{2} \geq \frac{1}{2} (\xi -1)^{2},$$ (13)

akzeptiere $\xi$ als Veränderliche, sonst beginne von vorne. Falls

$$u \leq 0.5,$$ (14)

ersetze $\xi$ durch $-\xi$. $\xi$ ist dann eine Veränderliche aus $N[0,1]$. JNormal2 zeigt das Verfahren in der zweiten Formulierung an Hand einer Unteroutine.

Approximation mit Hilfe des zentralen Grenzwertsatzes. Wir hatten früher bereits festgestellt, daß die Summe

$$x = u_{1} + u_{2} + ... + u_{n}$$

von $n$ unabhängigen gleichverteilten Zufallszahlen aus dem Intervall $[-\frac{1}{2},\frac{1}{2}]$ im Grenzwert $n \to \infty$ gegen eine Normalverteilung mit Mittelwert $<x> = a = 0$ und Varianz $\sigma = \sqrt{n/12}$ konvergiert. Wir hatten darüberhinaus festgestellt, daß man bereits für $n=12$ eine außerordentlich gute Approximation der Normalverteilung mit $\sigma = 1$ erhält. JNormal3 zeigt das Programm für diesen speziellen Fall. Ersichtlich glänzt dieses Programm durch seine Schlichtheit hinsichtlich mathematischer Operationen.

Eine weitere Approximation. Durch numerisches Rechnen kann man sich leicht davon überzeugen, daß für $x > 0$ die Approximation

$$e^{-x^{2}/2} \approx \frac{4 e^{-kx}}{[1+e^{-kx}]^{2}}$$

sehr gut erfüllt ist, sofern $k = \sqrt{8/\pi}$. Wir setzen also

$$\sqrt{\frac{2}{\pi}} e^{-x^{2}/2} \approx \frac{2 k e^{-kx}... ...{-kx}]^{2}} = q(x), \; \; \; x > 0, \; k=\sqrt{\frac{8}{\pi}}.$$ (15)

Die integrale Verteilungsfunktion für diese Approximation erhält man am einfachsten mit der Substitution $y = e^{-kx}$. Dieses führt auf

$$Q(\xi) = \int_{0}^{\xi} dx q(x) = \frac{2}{1+e^{-k\xi}} -1.$$ (16)

Mit Hilfe der inversen Transformations- Methode muß die Gleichung

$$Q(\xi) = \frac{2}{1+e^{-k\xi}} - 1 = u$$ (17)

mit $u \in U[0,1]$ nach $\xi$ aufgelöst werden. Dieses ergibt

$$\xi = \frac{1}{k} ln \left( \frac{1+u}{1-u} \right) = 0.56399 \; ln \left( \frac{1+u}{1-u} \right) .$$ (18)

Anschließend wird $\xi$ mit einem zufälligen Vorzeichen versehen. Diesen Algorithmus zeigt JNormal4.

Die $n$-dimensionale Normalverteilung. Die allgemeine $n$- dimensionale Normalverteilung.

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

kann mit Hilfe der Transformation $\vec{y} = C \vec{x} + \vec{a}$ auf die $n$- dimensionale unabhängige Standard- Normalverteilung

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

transformiert werden. Es genügt also, $n$ Veränderliche $x_{1},x_{2},...,x_{n} \in N[0,1]$ zu erzeugen und zum Vektor $\vec{x} = (x_{1},x_{2},...,x_{n})$ zusammenzufassen. Eine Programmroutine ist in JNormalN1 gelistet. In dieser Routine verwenden wir den Algorithmus JNormal1 zur Simulation einer einfachen normalverteilten Veränderlichen. Dieser Teil kann natürlich durch jeden anderen besprochenen Algorithmus ersetzt werden.

Die inverse $\chi^{2}$- Methode. Bei dem soeben beschriebenen Verfahren muß $n$- mal eine eindimensionale normalverteilte Veränderliche erzeugt werden. Wir diskutieren daher noch ein weiteres Verfahren, das zwar vom logischen Aufbau komplizierter ist, aber mit wesentlich weniger mathematischen Operationen auskommt. Wir erinnern in diesem Zusammenhang an die Definition der $\chi^{2}$- Verteilung: Wenn $x_{1},x_{2},...,x_{n} \in N[0,1]$, dann ist die Größe

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

$\chi^{2}$- verteilt mit $n$ Freiheitsgraden,

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

Diese Vorschrift kann umgekehrt werden. Dazu erzeugen wir zunächst eine $\chi^{2}$- verteilte Größe $y \in Ch[n]$. Einen hierzu geeigneten Algorithmus JChi2 werden wir später diskutieren. Die Größe $y$ ist dann nichts anderes als das Quadrat des Radius $r$ des $n$- dimensionalen Vektors $\vec{x} = (x_{1},x_{2},...,x_{n})$,

$$r \equiv \sqrt{y} = \sqrt{\sum_{\nu=1}^{n} x_{\nu}^{2}}.$$

Als nächstes erzeugen wir eine Gleichverteilung von Zufallszahlen auf der Oberfläche der $n$- dimensionalen Einheitskugel. Diese bezeichnen wir mit $\vec{z}=(z_{1},z_{2},...,z_{n})$, sodaß also

$$\sum_{\nu=1}^{n} z_{\nu}^{2} = 1.$$

Den Vektor $\vec{x}$ der Standard- Normalverteilung erhält man dann aus

$$\vec{x} = r \vec{z}.$$

Das Verfahren zur Erzeugung einer Gleichverteilung von Zufallszahlen auf der Oberfläche der Einheitskugel hatten wir im Prinzip bereits behandelt, zumindest für den Fall $n=2$ und $n=3$. Wir erzeugen zunächst eine Gleichverteilung $u_{1},u_{2},...,u_{n} \in U[-1,+1]$ im Volumen des $n$- dimensionalen Quaders und verwerfen alle Punkte mit

$$t = \sqrt{u_{1}^{2}+u_{2}^{2}+...+u_{n}^{2}} \geq 1.$$

Dieses liefert offensichtlich eine Gleichverteilung im Volumen der $n$- dimensionalen Kugel. Die Normierung gemäß

$$z_{\nu} = \frac{u_{\nu}}{t}, \; \; \; \nu=1,2,...,n,$$

führt dann auf eine Gleichverteilung auf der Oberfläche der Kugel.

Das vollständige Verfahren lautet: Erzeuge $u_{1},u_{2},...,u_{n} \in U[-1,+1]$ und berechne

$$t^{2} = \sum_{\nu=1}^{n} u_{\nu}^{2}.$$ (21)

Wenn

$$t^{2} \leq 1,$$ (22)

akzeptiere $t^{2}$, sonst beginne von vorne. Berechne $t=\sqrt{t^{2}}$ und setze

$$z_{\nu} = \frac{u_{\nu}}{t}, \; \; \; \nu=1,2,...,n.$$ (23)

Erzeuge $y \in Ch[n]$, berechne $r = \sqrt{y}$ und setze

$$x_{\nu} = r z_{\nu}, \; \; \; \nu=1,2,...,n.$$ (24)

Der zufällige Vektor $\vec{x} = (x_{1},x_{2},...,x_{n})$ ist dann eine Veränderliche der $n$- dimensionalen unabhängigen Standard- Normalverteilung.

Das Programm zu diesem Verfahren ist in JNormalN2 gelistet. Die wesentlichen mathematischen Operationen sind zwei Wurzelausdrücke und die Simulation der $\chi^{2}$- verteilten Größe $y$. Das Verfahren wäre außerordentlich schnell, wenn nicht die Verwerfungsbedingung (22) enthalten wäre. Die Erfolgswahrscheinlichkeit ist beschränkt auf

$$\frac{1}{C} = \frac{Volumen \; der \; Einheitskugel}{Volumen... ... = \frac{1}{n 2^{n-1}} \frac{\pi^{n/2}}{\Gamma (\frac{n}{2})}.$$ (25)

Ersichtlich gilt für große $n$:

$$\lim_{n \to \infty} \frac{1}{C} = 0.$$ (26)

Zum Beispiel ist für $n=2$ die Erfolgswahrscheinlichkeit $1/C \approx 0.79$ und für $n=10$ nur $\approx 0.0025$. Die Routine JNormalN2 ist daher nur für Werte von $n$ zwischen 2 und ungefähr 6 schneller als die Routine JNormalN1.

Lognormalverteilung.
Eine aus der allgemeinen Normalverteilung $N[a,\sigma ]$ durch die Transformation

$$y = e^{x}$$

hergeleitete Verteilung ist die Lognormalverteilung $LN[a,\sigma ]$:

$$p(y) = \frac{1}{\sqrt{2\pi} \sigma y} e^{-(ln y -a)^{2}/2\sigma^{2}}, \; \; \; 0 \leq y.$$ (27)

Der übliche Algorithmus ist daher: Erzeuge $x \in N[0,1]$ und setze

$$\eta = e^{a + \sigma x}.$$

$\eta$ ist dann eine zufällige Veränderliche aus der Lognormalverteilung $LN[a,\sigma ]$ (JLogNormal1).

Die Lognormalverteilung hat vielfache Anwendung in der Ökonomie bei der Darstellung von Einnahmen und Ausgaben. Erwartungswert und Varianz sind durch

$$<y> = e^{a + \sigma^{2}/2},$$

$$\sigma_{y} = \sqrt{<(y-<y>)^{2}>} = e^{2a + \sigma^{2}} (e^{\sigma^{2}}-1)$$

gegeben.

Exponentialverteilung
Zerfall quantenmechanischer Systeme. Während die Normalverteilung eine zentrale Rolle im Rahmen der mathematischen Statistik spielt, ist die Exponentialverteilung sicherlich die wichtigste Verteilung der mikroskopischen Physik, in der Physik also, die mit Hilfe der Quantenmechanik beschrieben werden muß. Ein in Ruhe befindliches instabiles quantenmechanisches System (Atom, Kern, Teilchen) wird mit Hilfe der Wellenfunktion

$$\psi(t) = e^{-i M t - \Gamma t /2} \psi(0)$$ (28)

beschrieben, wobei $M$ die Masse, $\tau= 1/\Gamma$ die mittlere Lebensdauer und $\psi(0)$ die Wellenfunktion am Zeitpunkt $t=0$ ist. Die Größe $\psi(0)$ beschreibt alle anderen, nicht von der Zeit abhängenden Größen des Systems. Die Wahrscheinlichkeit dafür, daß wir das System am Zeitpunkt $t$ noch unverändert antreffen, ist nach den Grundaxiomen der Quantenmechanik durch

$$P(t) = \vert \psi(t) \vert^{2} = \psi(t) \psi^{\ast}(t)$$ (29)

gegeben. Einsetzen von $\psi(t)$ aus (28) ergibt

$$P(t) = e^{-\Gamma t} \vert\psi(0)\vert^{2}.$$ (30)

Wir sehen also, daß die Zerfallswahrscheinlichkeit eines instabilen quantenmechanisches Systems einer Exponentialverteilung gehorcht.

Die normierte allgemeine Exponentialverteilung

$$p(y) = \frac{1}{\tau} e^{-y/\tau}$$ (31)

erhält man aus der Standard- Normalverteilung

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

mit Hilfe der Transformation $y = \tau x$. Wir bezeichnen die Exponentialverteilung mit dem mnemotechnischen Symbol $E[\tau]$. Die erzeugenden und logarithmisch erzeugenden Funktionen sind

$$M_{y}(v) = \frac{1}{1-v \tau}, \; \; \; v \tau < 1,$$ (33)

$$H_{y}(v) = - ln(1 - v \tau ), \; \; \; v \tau < 1.$$ (34)

Aus der Reihenentwicklung der logarithmisch erzeugenden Funktion,

$$H_{y}(v) = \sum_{\nu=1}^{\infty} \frac{\tau^{\nu}}{\nu} v^{\nu} = \tau v + \frac{1}{2} \tau^{2} v^{2} + .....$$

folgen sofort die Ausdrücke für Erwartungswert und Varianz:

$$<y> = \tau,$$ (35)

$$\sigma_{y} = \sqrt{<(y-<y>)^{2}>} = \tau.$$ (36)

Inverse Transformationsmethode. Zur Simulation einer Veränderlichen der Standard- Exponentialverteilung $E[1]$ verwendet man fast immer die inverse Transformationsmethode. Das Verfahren wurde bereits ausgiebig diskutiert. Wir können daher den Algorithmus sofort anschreiben: Erzeuge $u \in U[0,1]$ und setze

$$\xi = - ln(u).$$

$\xi$ ist dann eine Veränderliche aus $E[1]$. Die Programmroutine JExponential ist trivial und braucht nicht weiter erläutert zu werden.

Obwohl dieser Algorithmus bestechend einfach ist, sind vielfältige andere Algorithmen entwickelt worden, u.a. mit Hilfe der Composition- Methode, der Verwerfungsmethode und der Forsythe- Methode. Allen diesen Algorithmen gemeinsam ist, daß sie die Berechnung der Logarithmus- Funktion vermeiden, dafür aber logische Operationen einführen müssen. Diese Verfahren sind auf den heute zur Verfügung stehenden Rechnern im Vergleich zu JExponential nicht empfehlenswert, wir haben daher auf eine weitere Diskussion verzichtet. Die ersten $\mu$- Prozessoren dagegen, wie zum Beispiel der legendäre INTEL 8080, waren nur für logische Operationen und einfache arithmetische Algebra ausgelegt. Hier waren andere Verfahren durchaus von Wert.

n-dimensionale unabhängige Exponentialverteilung. Die Simulation der allgemeinen n- dimensionalen Standard- Exponentialverteilung

$$p(x_{1},x_{2},...,x_{n}) = \prod_{\nu=1}^{n} e^{-x_{\nu}} = e^{-\sum_{\nu=1}^{n} x_{\nu}}$$ (37)

erfordert die Berechnung von $n$ Logarithmus- Funktionen,

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

mit $u_{1},u_{2},...,u_{n} \in U[0,1]$. Für gewisse Werte von $n$ kann hierfür ein schnellerer Algorithmus angewendet werden, den wir im folgenden herleiten. Das Verfahren beruht auf einer Reduktionsmethode mit nachfolgender Variablen- Transformation.

Wir erzeugen zunächst $n$ gleichverteilte Zufallszahlen $u_{1},u_{2},...,u_{n}$ und bilden

$$z = - ln \prod_{\nu=1}^{n} u_{\nu}.$$

Wir behaupten jetzt, daß $z$ durch die sogenannte Erlang- Verteilung $Er[n,1]$,

$$p(z) = \frac{z^{n-1} e^{-z}}{(n-1)!}, \; \; \; z > 0,$$ (38)

beschrieben wird. Zum Beweis bemerken wir, daß die logarithmisch erzeugende Funktion der Summe von $n$ Veränderlichen

$$z = \sum_{\nu=1}^{n} z_{\nu}$$

der Exponentialfunktion

$$p(z_{\nu}) = e^{-z_{\nu}}$$

durch

$$H_{z}(v) = \sum_{\nu=1}^{n} H_{z_{\nu}}(v) = - n ln(1 - v)$$

gegeben ist. Auf Grund der Ergebnisse in Kap.3 ist dieses die logarithmisch erzeugende Funktion von

$$p(z) = \frac{z^{n-1} e^{-z}}{(n-1)!}.$$

Da die $z_{\nu} = -ln(u_{\nu})$ exponentiell verteilt sind, falls $u_{\nu} \in U[0,1]$, folgt unsere Behauptung, daß

$$z = \sum_{\nu=1}^{n} z_{\nu} = - ln \prod_{\nu=1}^{n} u_{\nu}.$$

Als nächstes erzeugen wir weitere $n-1$ Zufallszahlen $u_{n+1},u_{n+2},...,u_{2n-1} \in U[0,1]$ und sortieren diese in aufsteigender Reihenfolge. Diese letzteren Zahlen bezeichnen wir mit $u_{(1)},u_{(2)},...,u_{(n-1)}$. Weiterhin setzen wir $u_{(0)}=0$ und $u_{(n)}=1$. Dann gilt also:

$$u_{(0)} \leq u_{(1)} \leq u_{(2)} \leq ..... \leq u_{(n-1)} \leq u_{(n)}.$$

Wir setzen

$$y_{k} = \lbrace \begin{array}{lll} u_{(k)} - u_{(k-1)} & , & k=1,2,...,n-1 \\ z & , & k=n. \end{array}$$

Wie man zeigen kann (siehe geordntete Statistik in Kap.3), ist der $(n-1)$- dimensionale Vektor $(y_{1},y_{2},...,y_{n-1})$ gemäß der Dichtefunktion

$$p(y_{1},y_{2},...,y_{n-1}) = (n-1)!, \; \; \; \; \; \sum_{k=1}^{n-1} y_{k} \leq 1, \; y_{k} \geq 0, k=1,2,...,n-1,$$

verteilt. Da die Zufallszahlen $u_{1},u_{2},...,u_{n}$ unabhängig von den $(n-1)$ Zufallszahlen $u_{(1)},u_{(2)},...,u_{(n-1)}$ sind, so ist auch $y_{n} = z$ unabhängig von den Veränderlichen $y_{1},y_{2},...,y_{n-1}$. Daher ist die gemeinsame Dichtefunktion des Vektors $\vec{y} = (y_{1},y_{2},...,y_{n})$ als Produkt darstellbar:

$$p(y_{1},y_{2},...,y_{n}) = y_{n}^{n-1} e^{-y_{n}}, \; \; \; ... ... \sum_{k=1}^{n-1} y_{k} \leq 1, \; y_{k} \geq 0, k=1,2,...,n.$$

Mit Hilfe der Transformation

$$x_{k} = \lbrace \begin{array}{lll} y_{k} y_{n} & , & k=1,2,...,n-1 \\ (1-y_{1}-y_{2}-...-y_{n-1}) y_{n} & , & k=n \end{array}$$

erhalten wir schließlich

$$p(x_{1},x_{2},...,x_{n}) = \prod_{k=1}^{n} e^{-x_{k}}, \; \; \; \; \; x_{k} \geq 0, \; k=1,2,...,n,$$

also die einfache $n$- dimensionale unabhängige Standard- Exponentialverteilung.

Zusammengefaßt erhalten wir folgenden Algorithmus: Erzeuge $u_{1},u_{2},...,u_{n} \in U[0,1]$ und bilde

$$z = - ln \prod_{\nu=1}^{n} u_{\nu}.$$

Erzeuge weitere $(n-1)$ Zufallszahlen $u_{n+1},u_{n+2},...,u_{2n-1} \in U[0,1]$, sortiere diese in aufsteigender Reihenfolge und bezeichne diese mit $u_{(1)},u_{(2)},...,u_{(n-1)}$. Setze außerdem $u_{(0)}=0$ und $u_{(n)}=1$. Bilde

$$\xi_{k} = (u_{(k)} - u_{(k-1)}) z, \; \; \; k=1,2,...,n.$$

$\vec{\xi} = (\xi_{1},\xi_{2},...,\xi_{n})$ ist dann ein Vektor der unabhängigen $n$- dimensionalen Standard- Exponentialverteilung $E[1,1,...,1]$.

Ersichtlich muß in diesem Verfahren nur einmal der Logarithmus berechnet werden, im Gegensatz zu $n$ Logarithmus- Berechnungen bei der inversen Transformationsmethode. Demgegenüber stehen doppelt so viele Aufrufe des Zufallszahlen- Generators und eine zusätzliche Sortierung von $(n-1)$ Zahlen. Der letztere Algorithmus erweist sich daher nur für Werte von $n$ zwischen etwa 3 und 6 der inversen Transformationsmethode überlegen. Programmroutinen sind in JExponentialN1 und JExponentialN2 gelistet.

Gammaverteilung
Vorbemerkung. Die Verallgemeinerung der Exponentialverteilung und der bereits im vorherigen Abschnitt behandelten Erlang- Verteilung ist die Gammaverteilung $G[\alpha, \beta]$,

$$p(y) = \frac{y^{\alpha -1} e^{-y/\beta}}{\beta^{\alpha} \Gam... ...lpha)}, \; \; \; \; \; \alpha > 0, \; \beta > 0, \; 0 \leq y.$$ (39)

Ersichtlich genügt es, Algorithmen zur Simulation der Standard- Gammaverteilung $G[\alpha, 1]$,

$$p(x) = \frac{x^{\alpha -1} e^{-x}}{\Gamma(\alpha)}, \; \; \; \alpha > 0, \; 0 \leq x,$$ (40)

zu diskutieren, da die Transformation $y = \beta x$ eine Veränderliche der allgemeinen Gammaverteilung liefert. Wie bereits diskutiert, erhalten wir für $\alpha =1$ die Exponentialverteilung $E[\beta]$ und für ganzzahliges $\alpha = m$ die Erlang- Verteilung $Er[m,\beta]$. Die erzeugende Funktion der Gammverteilung ist

$$M_{y}(v) = \frac{1}{(1-v \beta)^{\alpha}}$$ (41)

und die logarithmisch erzeugende Funktion ist

$$H_{y}(v) = - \alpha ln(1 - v \beta).$$ (42)

Erwartungswert und Varianz erhält man am einfachsten aus der Reihenentwicklung der logarithmisch erzeugenden Funktion

$$H_{y}(v) = \alpha \sum_{\nu=1}^{\infty} \frac{\beta^{\nu}}{\... ... \beta v + \frac{1}{2} (\sqrt{\alpha} \beta)^{2} v^{2} + .....$$ (43)

zu

$$<y> = \alpha \beta,$$ (44)

$$\sigma_{y} = \sqrt{<(y-<y>)^{2}>} = \sqrt{\alpha} \beta .$$ (45)

Reproduktivität der Gammaverteilung. Die wichtigste Eigenschaft der Gammaverteilung ist ihre Reproduktivität: Seien $x_{1}$ und $x_{2}$ zwei Veränderliche aus Gammaverteilungen mit Parametern $\alpha_{1}$ und $\alpha_{2}$, d.h. $x_{1} \in G[\alpha_{1},\beta]$ und $x_{2} \in G[\alpha_{2},\beta]$. Dann ist $x=x_{1}+x_{2}$ eine Veränderliche aus der Gammaverteilung $G[\alpha_{1}+\alpha_{2}, \beta]$. Der Beweis ist einfach und erfolgt über die logarithmisch erzeugende Funktion.

$$H_{x_{1}}(v) = -\alpha_{1} ln(1-\beta v),$$

$$H_{x_{2}}(v) = -\alpha_{2} ln(1-\beta v)$$

sind die logarithmisch erzeugenden Funktionen von

$$p(x_{1}) = \frac{x_{1}^{\alpha_{1}-1} e^{-x_{1}/\beta}} {\beta^{\alpha_{1}} \Gamma(\alpha_{1})},$$

$$p(x_{2}) = \frac{x_{2}^{\alpha_{2}-1} e^{-x_{2}/\beta}} {\beta^{\alpha_{2}} \Gamma(\alpha_{2})}.$$

Nach Kap.3 ist dann, da $x_{1}$ und $x_{2}$ unabhängig voneinander sind,

$$H_{x}(v) = H_{x_{1}}(v) + H_{x_{2}}(v) = -(\alpha_{1} + \alpha_{2}) ln(1 - \beta v)$$

die logarithmisch erzeugende Funktion der Summe $x=x_{1}+x_{2}$ der zwei Veränderlichen $x_{1}$ und $x_{2}$. Hieraus folgt sofort unsere Behauptung. Wichtig ist offensichtlich, daß der Parameter $\beta$ in beiden Verteilungen identisch ist.

Die Verallgemeinerung auf $n$ Veränderliche $x_{1}, x_{2},...,x_{n}$ aus den Gammaverteilungen $G[\alpha_{1},\beta]$, $G[\alpha_{2}, \beta]$, ...., $G[\alpha_{n}, \beta]$ ist trivial und ergibt, daß die Veränderliche $x=x_{1}+x_{2}+...+x_{n}$ aus $G[\alpha_{1}+\alpha_{2}+...+ \alpha_{n}, \beta]$ ist.

Reproduktivitäts- Verfahren. Ein auf der Reproduktivitätseigenschaft basierendes Verfahren soll im folgenden vorgestellt werden. Wir setzen

$$\alpha = \sum_{\nu=1}^{m} \alpha_{\nu} + \delta,$$

wobei

$$\alpha_{\nu} = 1, \; \; \; \nu = 1,2,...,m,$$

$$0 < \delta < 1,$$

und erzeugen $m$ Zufallszahlen $x_{1},x_{2},...,x_{m}$ aus der Exponentialverteilung $E[1]$,

$$p(x_{\nu}) = e^{-x_{\nu}}, \; \; \; \nu=1,2,...,m,$$

sowie eine Zufallszahl $z$ aus der speziellen Gammaverteilung $G[\delta,1]$,

$$p(z) = \frac{z^{\delta -1} e^{-z}}{\Gamma(\delta)}, \; \; \; 0 \leq z,$$

mit

$$0 < \delta < 1 .$$

Auf Grund der Reproduktivität ist dann

$$x = \sum_{\nu=1}^{m} x_{\nu} + z$$

eine Veränderliche der Standard- Gammaverteilung $G[\alpha, 1]$

$$p(x) = \frac{x^{\alpha -1} e^{-x}}{\Gamma(\alpha)}, \; \; \; 0 \leq x.$$

Die Simulation der speziellen Gammaverteilung $G[\delta,1]$ mit $0 < \delta < 1$ kann wie folgt durchgeführt werden.

Hierzu führen wir noch die Betaverteilung $B[\alpha, \beta]$ ein. Diese ist allgemein definiert durch

$$p(y) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\bet... ...ta-1}, \; \; \; \alpha > 0, \; \beta > 0, \; 0 \leq y \leq 1.$$

Wir benötigen im folgenden die spezielle Betaverteilung $B[\delta, 1-\delta]$,

$$p(w) = \frac{1}{\Gamma(\delta)\Gamma(1-\delta)} w^{\delta -1} (1-w)^{-\delta}.$$

Wegen $\delta > 0$ und $1-\delta > 0$ ist diese Verteilung ersichtlich nur für $0 < \delta < 1$ definiert. Seien also $w \in B[\delta,1-\delta]$ und $v \in E[1]$, dann erhalten wir mit Hilfe einer Variablen- Transformation für die Veränderlichen

$$u = v,$$

$$z = v w,$$

die Dichtefunktion

$$p(u,z) = \frac{1}{\Gamma(\delta)\Gamma(1-\delta)} z^{\delta-1} (u-z)^{-\delta} e^{-u}, \; \; \; 0 \leq z, \; 0 \leq u.$$

Die Randverteilung in $z$ ist

$$p(z) = \frac{1}{\Gamma(\delta)} z^{\delta-1} e^{-z}, \; \; \; 0 \leq z.$$

Dieses ist offensichtlich die gesuchte Gammaverteilung $G[\delta,1]$.

Wir fassen das Verfahren zusammen: Setze $\alpha = m + \delta$ mit $0 < \delta < 1$. Erzeuge $u_{1},u_{2},...,u_{m} \in U[0,1]$ und bilde

$$y = - \sum_{\nu=1}^{m} ln u_{\nu} = - ln \prod_{\nu=1}^{m} u_{\nu}.$$

Erzeuge $w \in B[\delta,1-\delta]$ und $v \in E[1]$. Die Veränderliche

$$\xi = y + v w$$

ist dann eine Veränderliche aus $G[\alpha, 1]$ mit $\alpha = m + \delta$, $0 < \delta < 1$. Das Programmlisting JGamma2 im zweiten Applet dieses Kapitels zeigt diesen Algorithmus. JGamma1 zeigt den einfachen Programmcode für ganzzahliges $\alpha = m$. Diese Routine haben wir auch unter dem Namen JErlang1 hinzugefügt..

Das soeben diskutierte Verfahren erfordert im wesentlichen die Erzeugung einer Veränderlichen aus der Betaverteilung. Wie wir im nächsten Abschnitt sehen werden, kann dieses zu erheblichem Rechenzeit- Bedarf führen. Günstiger sind daher im allgemeinen Methoden, die auf dem Verwerfungsverfahren beruhen.

Approximation für große $\alpha$. Für große Werte von $\alpha$ bzw $m$ kann die Erlang- Verteilung $Er[m,1]$ durch eine Normalverteilung ersetzt werden. Dieses folgt aus der Tatsache, daß eine Veränderliche der Erlang- Verteilung als Summe von $m$ Veränderlichen der Exponentialverteilung $E[1]$ geschrieben werden kann. Erwartungswert und Varianz der Erlang- Verteilung sind für große Werte von $m$ durch

$$<x> = m,$$

$$\sigma_{x} = \sqrt{m}$$

gegeben. Die Approximation kann also geschrieben werden als

$$p(x) \approx \frac{1}{\sqrt{2 \pi m}} e^{-(x-m)^{2}/2m}.$$

Der Algorithmus lautet: Erzeuge $v \in N[0,1]$, dann ist

$$\xi = m + v \sqrt{m}$$

eine Veränderliche der Erlang- Verteilung $Er[m,1]$ für große $m$ ($m > 10$). Diese Routine ist unter dem Namen JErlang2 gelistet.

1. Verwerfungsverfahren. Als erstes betrachten wir ein Verfahren, in dem die Simulation der speziellen Gammaverteilung

$$p(z) = \frac{1}{\Gamma(\delta)} z^{\delta -1} e^{-z}, \; \; \; \; \; 0 < \delta < 1, \; 0 \leq z,$$

mit Hilfe eines Verwerfungsverfahrens durchgeführt wird. Der Algorithmus hierzu wurde bereits ausführlich diskutiert. Wir fassen das Verfahren noch einmal in einer mathematisch etwas abgewandelten Formulierung zusammen:
1. Berechne $m$ und $\delta$ gemäß

$$\alpha = m + \delta, \; \; \; \; \; 0 < \delta < 1.$$

2. Erzeuge $u_{1},u_{2},...,u_{m} \in U[0,1]$ und setze

$$y = - ln \prod_{\nu=1}^{m} u_{\nu}.$$

Erzeuge $u_{m+1},u_{m+2} \in U[0,1]$ und berechne

$$\kappa = ln(\frac{\delta + e}{e}),$$

$$ln(z) = \frac{1}{\delta} (\kappa + ln(u_{m+1})).$$

Wenn $ln(z) > 0$, fahre fort bei Punkt 3. Ansonsten berechne $z = e^{ln(z)}$. Falls

$$ln(u_{m+2}) \leq -z,$$

akzeptiere $\xi = y + z$ als Veränderliche der Gammaverteilung $G[\alpha, 1]$, sonst wiederhole von Punkt 2.
3. Setze

$$z = 1 - \kappa - ln(1-u_{m+1}).$$

Wenn

$$ln(u_{m+2}) \leq (1-\delta) ln(z),$$

akzeptiere $\xi = y + z$ als Veränderliche der Gammaverteilung $G[\alpha, 1]$, sonst wiederhole bei Punkt 2. JGamma3 zeigt den Programm- Code.

2. Verwerfungsverfahren. Das zweite Verwerfungsverfahren beruht auf dem gleichen mathematischen Formalismus wie die Routine JGam2, nur daß jetzt die Simulation der Betaverteilung durch ein anderes Verwerfungsverfahren ersetzt wird. Die Funktionen

$$h(x) = (1-\delta)\frac{x^{m-1}e^{-x}}{(m-1)!} + \delta \frac{x^{m}e^{-x}} {m!},$$

$$G(x) = \frac{(\frac{x}{m})^{\delta}}{1+\delta (\frac{x}{m}-1)},$$

$$C = \frac{(m-1)! m^{\delta}}{\Gamma(m+\delta)},$$

wobei, wie bisher, $\alpha = m + \delta$, $0 < \delta < 1$, erfüllen ersichtlich die Bedingungen der Verwerfungsmethode. Die Simulation der Dichtefunktion $h(x)$ ist einfach und beruht auf der alternativen Simulation zweier Erlang- Verteilungen $Er[m-1,1]$ und $Er[m,1]$ mit der Wahrscheinlichkeit $(1-\delta)$ bzw $\delta$. Die Erfolgswahrscheinlickeit $1/C$ konvergiert für große Werte von $\alpha$ bzw $m$ gegen 1:

$$\lim_{m \to \infty} \frac{1}{C} = \lim_{m \to \infty} \frac{\Gamma(m+\delta)} {(m-1)! m^{\delta}} = 1.$$

Wegen $\Gamma(m) = (m-1)!$ für ganzzahliges positives $m$ ist für $\delta =0$ die Erfolgswahrscheinlichkeit immer gleich 1. Der Algorithmus lautet:
1. Berechne $m$ und $\delta$ gemäß

$$\alpha = m + \delta, \; \; \; \; \; 0 < \delta < 1.$$

2. Erzeuge $u_{1},u_{2},...,u_{m},u_{m+1},u_{m+2},u_{m+3} \in U[0,1]$ und setze

$$\xi = \lbrace \begin{array}{lll} - ln \prod_{\nu=1}^{m} u_{\... ...{m+1} u_{\nu} & \; f''ur \; & u_{m+2} \leq \delta. \end{array}$$

Falls

$$u_{m+3} \leq \frac{(\frac{\xi}{m})^{\delta}}{1 + \delta (\frac{\xi}{m} -1)},$$

akzeptiere $\xi$ als Veränderliche der Gammaverteilung $G[\alpha, 1]$, sonst beginne bei Punkt 2. JGamma4 zeigt den Programmcode.

3. Verwerfungsverfahren. Ein drittes Verwerfungsverfahren beruht auf den Festsetzungen

$$h(x) = \frac{1}{\alpha} e^{-x/\alpha},$$

$$G(x) = x^{\alpha -1} \frac{e^{-x(1-1/\alpha)}} {\alpha^{\alpha -1} e^{1-\alpha}},$$

$$C = \frac{\alpha^{\alpha} e^{1 -\alpha}}{\Gamma(\alpha)},$$

und ist nur anwendbar für $\alpha \geq 1$. Die Erfolgswahrscheinlichkeit $1/C$ kann für große $\alpha$ durch

$$\frac{1}{C} \approx \frac{1}{e} \sqrt{\frac{2 \pi}{\alpha}}$$

approximiert werden. Der Algorithmus lautet: Erzeuge $u_{1},u_{2} \in U[0,1]$ und berechne

$$\xi = - \alpha ln(u_{1}).$$

Wenn

$$ln(u_{2}) \leq (\alpha-1) (ln(\xi) - \frac{\xi}{\alpha} - ln(\alpha) + 1),$$

akzeptiere $\xi$ als Veränderliche von $G[\alpha, 1]$, sonst beginne von vorne. Hierbei haben wir, wie auch schon in vielen vorherigen Beispielen, die komplizierte Funktion $x^{y}$ dadurch vermieden, daß auf beiden Seiten der Verwerfungsabfrage der Logarithmus gebildet wurde. Auf den meisten Rechnern ist der Logarithmus $ln(x)$ und die Exponentialfunktion $e^{x}$ eine Standard- Routine, während ein Ausdruck der Form $x^{y}$ mit Hilfe der Umformung $x^{y} = e^{y ln(x)}$ berechnet wird, also aus zwei Operationen besteht. Die Routine JGamma5 zeigt den Algorithmus.

4. Verwerfungsverfahren. Die nächste Prozedur faktorisiert $p(x)$ in der Form

$$h(x) = \frac{\lambda \mu x^{\lambda -1}}{(\mu+x^{\lambda})^{2}}, \; \; \; \; \; 0 \leq x,$$

$$G(x) = x^{\alpha -\lambda} (\mu+x^{\lambda})^{2} \frac{e^{\alpha -x}} {4 \alpha^{\alpha + \lambda}},$$

$$C = \frac{4 \alpha^{\alpha}}{\Gamma(\alpha) e^{\alpha} \lambda},$$

mit $\mu=\alpha^{\lambda}$ und $\lambda = \sqrt{(2\alpha -1)}$. Die Erfolgswahrscheinlichkeit ist mit

$$\frac{1}{C} = \lbrace \begin{array}{lll} 0.68 & \; f''ur \; ... ... \\ . & & \\ 0.88 & \; f''ur \; & \alpha \to \infty \end{array}$$

für alle Werte von $\alpha$ relativ gut. Die Simulation der Dichtefunktion $h(x)$ beruht auf der inversen Transformationsmethode. Wir transformieren gemäß

$$y = x^{\lambda}$$

und erhalten die Dichtefunktion der Veränderlichen $y$ zu

$$h(y) = \frac{\mu}{(\mu+u)^{2}}.$$

Die integrale Verteilungsfunktion ist

$$H(\eta) = 1 - \frac{\mu}{\mu+\eta}.$$

Mit $u_{1},u_{2} \in U[0,1]$ ergeben sich die Gleichungen

$$\eta = \frac{\mu u_{1}}{1-u_{1}},$$

$$\xi = \eta^{1/\lambda},$$

sowie die Ungleichung

$$u_{2} \leq \xi^{\alpha -\lambda} (\mu + \xi^{\lambda})^{2} \frac{e^{\alpha - \xi}}{4 \alpha^{\alpha + \lambda}}.$$

Die letzte Ungleichung zur Verwerfungsabfrage kann vereinfacht werden zu

$$ln(u_{1}^{2} u_{2}) \leq (\alpha + \lambda) ln(\xi) - \xi + ln(\frac{e^{\alpha}}{4 \alpha^{\alpha} \mu})$$

mit

$$ln(\xi) = \frac{1}{\lambda} ln(\frac{\mu u_{1}}{1-u_{1}}),$$

$$\xi = e^{ln (\xi)}.$$

Damit ergibt sich der Algorithmus: 1. Berechne in einer Initialisierungsroutine die Größen

$$\lambda = \sqrt{2 \alpha -1},$$

$$\mu = \alpha^{\lambda},$$

$$\epsilon = \alpha + \lambda,$$

$$\kappa = ln\left(\frac{e^{\alpha}}{4 \alpha^{\alpha} \mu}\right).$$

2. Erzeuge $u_{1},u_{2} \in U[0,1]$ und setze

$$v = \frac{1}{\lambda} ln(\frac{\mu u_{1}}{1-u_{1}}).$$

$$\xi = e^{v}.$$

Wenn

$$ln(u_{1}^{2} u_{2}) \leq \epsilon v - \xi + \kappa,$$

akzeptiere $\xi$ als Veränderliche von $G[\alpha, 1]$, sonst beginne von vorne bei Punkt 2. JGamma6 zeigt das Programmlisting, das aus einem Initialisierungsteil und dem Generator selbst besteht.

5. Verwerfungsverfahren. Ein weiteres Verwerfungsverfahren arbeitet mit Hilfe der Cauchy- Verteilung $C[\gamma, \lambda]$,

$$h_{C}(x) = \frac{1}{\pi \lambda [ 1 +(\frac{x-\gamma}{\lambda})^{2}]}, \; \; \; \; \; -\infty < x < \infty.$$

Diese Verteilung ist allgemein auf der gesamten reellen Achse definiert. Bei Einschränkung des Definitionsbereiches auf $x \geq 0$ muß die Normierung geändert werden. Man erhält

$$h(x) = \frac{h_{C}(x)}{1 - H_{C}(x)}, \; \; \; \; \; 0 \leq x,$$

wobei

$$H_{C}(x) = \int_{-\infty}^{\xi} dx h_{C}(x) = \frac{1}{2} + \frac{1}{\pi} arctg(\frac{\xi-\gamma}{\lambda})$$

die integrale Verteilungsfunktion der Cauchy- Verteilung ist. Mit Hilfe der inversen Transformationsmethode ergibt sich, daß

$$\xi = \gamma + \lambda \; tg\{\pi [u_{1} (1-H_{C}(0)) + H_{C}(0) -\frac{1}{2}] \}$$

eine Veränderliche von $h(x)$ ist, wobei, wie immer, $u_{1} \in U[0,1]$. Das Verwerfungsverfahren für die Gammaverteilung $G[\alpha, 1]$ erhält man mit den weiteren Festsetzungen

$$G(x) = [1 + (\frac{x-\gamma}{\lambda})^{2}] (\frac{x}{\gamma})^{\gamma} e^{-x + \gamma},$$

$$C = \pi \lambda \gamma^{\gamma} \frac{[1-H_{C}(0)]}{\Gamma(\alpha) e^{\gamma}},$$

mit $\gamma = \alpha -1$ und $\lambda = \sqrt{2 \alpha -1}$. Der Algorithmus lautet: 1. Berechne in einem externen Programm

$$\gamma = \alpha -1,$$

$$\lambda = \sqrt{2 \alpha -1},$$

$$\epsilon = \frac{\pi}{2} + arctg(\frac{\gamma}{\lambda}),$$

$$\rho = \frac{\pi}{2} - \epsilon.$$

2. Erzeuge $u_{1},u_{2} \in U[0,1]$ und setze

$$\xi = \gamma + \lambda \; tg\{ \epsilon u_{1} + \rho \}.$$

Wenn

$$ln(u_{2}) \leq ln\{ 1 + (\frac{\xi-\gamma}{\lambda})^{2} \} + \gamma ln(\frac{\xi}{\gamma}) - \xi + \gamma,$$

akzeptiere $\xi$ als Veränderliche von $G[\alpha, 1]$, sonst beginne von vorne bei Punkt 2. Die entsprechende Routine ist JGamma7.

6. Verwerfungsverfahren. Ein Verwerfungsverfahren für die Gammaverteilung $G[\alpha, \beta]$ mit $\beta \not= 1$ ist das folgende:

$$h(x) = \frac{x^{m-1} e^{-x/\beta}}{\beta^{m} (m-1)!}, \; \; \; \; \; \beta > 0, \alpha > 0, x \geq 0.$$

$$G(x) = \frac{x^{\delta} e^{-x(1-1/\beta)}}{[\delta \beta/(\beta-1)]^{\delta} e^{-\delta}}, \; \; \; \; \; \beta \not= 1,$$

$$C = \frac{\beta^{m} (m-1)!}{\Gamma(\alpha)} \left( \frac{\delta \beta}{\beta -1} \right)^{2} e^{-\delta}, \; \; \; \; \; \beta \not= 1,$$

mit $\alpha = m + \delta, 0 < \delta < 1$. Ersichtlich kann die Faktorisierung für $G(x)$ und $C$ nur für $\beta \not= 1$ durchgeführt werden. Das Verfahren ist ähnlich dem 3. Verwerfungsverfahren, nur daß jetzt die Gammaverteilung $G[\alpha, \beta]$ durch die allgemeine Erlang- Verteilung $Er[m,\beta]$ approximiert wird, anstatt durch die Exponentialverteilung $E[\alpha]$ im 3. Verwerfungsverfahren. Die Mathematik dieses Verfahrens ist sehr einfach, sodaß wir den Algorithmus sofort hinschreiben können: 1. Berechne $m$ und $\delta$ aus

$$\alpha = m + \delta, \; \; \; 0 < \delta < 1,$$

und setze

$$\gamma = 1 - \frac{1}{\beta},$$

$$\kappa = \delta \left( 1 - ln( \frac{\delta \beta}{\beta -1}) \right).$$

2. Erzeuge $u_{1},u_{2},...,u_{m} \in U[0,1]$ und setze

$$\xi = - \beta \; ln \prod_{\nu=1}^{m} u_{\nu}.$$

Wenn

$$ln(u_{2}) \leq \delta ln(\xi) - \gamma \xi + \kappa,$$

akzeptiere $\xi$ als Veränderliche von $G[\alpha, \beta]$, sonst beginne von vorne bei Punkt 2. Das Programm ist in JGamma8 gelistet.

Wir haben die Gamma- Verteilung benutzt, um die vielfältigen Methoden aufzuzeigen. Die Rechenzeiten für die einzelnen Methoden sind zwar nicht dramatisch verschieden, aber natürlich schon etwas unterschiedlich. Im nächsten Abschnitt werden wir die Simulation weiterer kontinuierlicher Verteilungen besprechen.