Kontinuierliche Verteilungen II
Übersicht
Nachdem wir im vorherigen Kapitel die Normalverteilung, die Exponentialverteilung und die Gammaverteilung diskutiert haben, setzen wir im vorliegenden Kapitel die Behandlung der kontinuierlichen Verteilungen fort und beginnen mit der Betaverteilung. Zunächst jedoch müssen Sie das Applet auf Ihren Bildschirm laden.

Betaverteilung
Die Betaverteilung $B[\alpha, \beta]$ ist definiert als

\begin{displaymath}
p(x) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\...
... \; \; \; \; \; \alpha > 0, \; \beta > 0, \;
0 \leq x \leq 1,
\end{displaymath} (1)

und ist eine der wenigen Verteilungen, bei denen die erzeugenden Funktionen nicht analytisch darstellbar sind. Die Momente dagegen sind leicht ausrechenbar und zwar

\begin{displaymath}
m_{j} = <x^{j}> = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alph...
...amma(\beta)}
\int_{0}^{1} dx x^{j+\alpha -1} (1-x)^{\beta -1}.
\end{displaymath}

Einer Integrationstabelle entnehmen wir die Formel

\begin{displaymath}
\int_{0}^{1} dx x^{a} (1-x)^{b} = \frac{\Gamma(a+1) \Gamma(b+1)}{\Gamma(a+b+2)}.
\end{displaymath}

Damit erhalten wir
\begin{displaymath}
m_{j} = \frac{\Gamma(\alpha + \beta) \Gamma(j+\alpha)}{\Gamma(\alpha)
\Gamma(j+\alpha + \beta)}.
\end{displaymath} (2)

Die erzeugende Funktion kann dann natürlich als unendliche Reihe dargestellt werden:

\begin{displaymath}
M_{x}(v) = \sum_{j=0}^{\infty} \frac{m_{j}}{j!} v^{j} =
\su...
...a)}
{\Gamma(j+1) \Gamma(\alpha) \Gamma(j+\alpha+\beta)} v^{j},
\end{displaymath}

wobei wir noch $\Gamma(j+1) = j!$ für ganzzahliges, positives $j$ gesetzt haben. Für Erwartungswert und Varianz erhalten wir
$\displaystyle m_{1}$ $\textstyle =$ $\displaystyle \frac{\alpha}{\alpha+\beta},$ (3)
$\displaystyle \sigma_{x}^{2} = m_{2} - m_{1}^{2}$ $\textstyle =$ $\displaystyle \frac{\alpha \beta}
{(\alpha + \beta)^{2} (\alpha + \beta +1)}.$ (4)

Simulation der speziellen Betaverteilung. Für den speziellen Fall $\alpha = m$ und $\beta = n$ mit ganzzahligem, poitivem $m$ und $n$ existiert ein sehr einfacher Simulations- Algorithmus aus der geordneten Statistik (siehe auch Kap.3). Seien $u_{1},u_{2},...u_{m+n-1} \in U[0,1]$ gleichverteilte Zufallszahlen und $u_{(1)},u_{(2)},...,u_{(m+n-1)}$ die gleichen Zahlen in geordneter, aufsteigender Reihenfolge, d.h.

\begin{displaymath}
u_{(1)} \leq u_{(2)} \leq ...... \leq u_{(m+n-1)}.
\end{displaymath}

Dann ist die Zahl $\xi = u_{(m)}$ eine Veränderliche der speziellen Betaverteilung $B[m,n]$,

\begin{displaymath}
p(x) = \frac{\Gamma(m+n)}{\Gamma(m)\Gamma(n)} x^{m-1} (1-x)^{n-1}.
\end{displaymath}

Die entsprechende Programmroutine ist JBeta1. Die SORT- Methode hatten wir bereits bei der Simulation der Gammaverteilung eingeführt.

Ein Verwerfungsverfahren. Aufbauend auf dieser Routine können wir jetzt ein Verwerfungsverfahren für die allgemeine Betaverteilung $B[\alpha, \beta]$ entwickeln. Mit $\alpha = m + \delta_{1}, \beta = n + \delta_{2},
0 < \delta_{1} < 1, 0 < \delta_{2} < 1$, erfüllt

$\displaystyle h(x)$ $\textstyle =$ $\displaystyle \frac{\Gamma(m+n)}{\Gamma(m) \Gamma(n)} x^{m-1} (1-x)^{n-1},$  
$\displaystyle G(x)$ $\textstyle =$ $\displaystyle \frac{(\delta_{1}+\delta_{2})^{\delta_{1}+\delta_{2}}}
{\delta_{1}^{\delta_{1}} \delta_{2}^{\delta_{2}}} x^{\delta_{1}}
(1-x)^{\delta_{2}},$  
$\displaystyle C$ $\textstyle =$ $\displaystyle \frac{\delta_{1}^{\delta_{1}} \delta_{2}^{\delta_{2}}}
{(\delta_{...
...a(m) \Gamma(n) \Gamma(\alpha+\beta)}{\Gamma(m+n) \Gamma(\alpha)
\Gamma(\beta)},$  

die Voraussetzungen des Verwerfungsverfahrens. Das Maximum der Funktion $G(x)$ liegt an der Stelle $x_{max} = \delta_{1}/(\delta_{1}+\delta_{2})$. Einsetzen in die obige Formel ergibt tasächlich $G(x_{max}) = 1.$ Die Erfolgswahrscheinlichkeit dieser Prozedur genügt der Ungleichung

\begin{displaymath}
\frac{(\delta_{1}+\delta_{2})^{\delta_{1}+\delta_{2}}}{(m+n+...
...{\delta_{1}} \delta_{2}^{\delta_{2}}} \leq \frac{1}{C} \leq 1.
\end{displaymath}

Der Algorithmus lautet also:
1. Berechne $m, n, \delta_{1}, \delta_{2}$ und $\epsilon$ gemäß

$\displaystyle \alpha$ $\textstyle =$ $\displaystyle m + \delta_{1}, \; \; \; \; \; 0 < \delta_{1} < 1,$  
$\displaystyle \beta$ $\textstyle =$ $\displaystyle n + \delta_{2}, \; \; \; \; \; 0 < \delta_{2} < 1,$  
$\displaystyle \epsilon$ $\textstyle =$ $\displaystyle \frac{(\delta_{1}+\delta_{2})^{\delta_{1}+\delta_{2}}}
{\delta_{1}^{\delta_{1}} \delta_{2}^{\delta_{2}}}.$  

2. Erzeuge $\xi \in B[m,n]$ und $u \in U[0,1]$. Falls

\begin{displaymath}
u \leq \epsilon \xi^{\delta_{1}} (1-\xi)^{\delta_{2}},
\end{displaymath}

akzeptiere $\xi$ als Veränderliche aus $B[\alpha, \beta]$, sonst beginne von vorne bei Punkt 2. JBeta2 zeigt den Algorithmus.

Reduktionsmethode aus der Gammaverteilung. Seien $y_{1}$ und $y_{2}$ zwei Veränderliche aus Gammaverteilungen $G[\alpha, 1]$ und $G[\beta, 1]$,

$\displaystyle p(y_{1})$ $\textstyle =$ $\displaystyle \frac{1}{\Gamma(\alpha)} y_{1}^{\alpha-1} e^{-y_{1}},$  
$\displaystyle p(y_{2})$ $\textstyle =$ $\displaystyle \frac{1}{\Gamma(\beta)} y_{2}^{\beta-1} e^{-y_{2}}.$  

Da $y_{1}$ und $y_{2}$ unabhängig voneinander sind, ist ihre gemeinsame Dichtefunktion:

\begin{displaymath}
p(y_{1},y_{2}) = \frac{1}{\Gamma(\alpha) \Gamma(\beta)} y_{1...
...ta-1} e^{-(y_{1}+y_{2})}, \; \; \; y_{1} \geq 0, y_{2} \geq 0.
\end{displaymath}

Mit Hilfe der Transformation
$\displaystyle x_{1}$ $\textstyle =$ $\displaystyle \frac{y_{1}}{y_{1}+y_{2}},$  
$\displaystyle x_{2}$ $\textstyle =$ $\displaystyle y_{1} + y_{2},$  

bzw deren Umkehrung
$\displaystyle y_{1}$ $\textstyle =$ $\displaystyle x_{1} x_{2},$  
$\displaystyle y_{2}$ $\textstyle =$ $\displaystyle x_{2} - x_{1} x_{2},$  

erhalten wir die Dichtefunktion für $x_{1}$ und $x_{2}$:

\begin{displaymath}
p(x_{1},x_{2}) = \frac{1}{\Gamma(\alpha) \Gamma(\beta)}
x_{...
... e^{-x_{2}},
\; \; \; \; \; 0 \leq x_{1} \leq 1, 0 \leq x_{2}.
\end{displaymath}

Die Randverteilung in $x_{1}$ ist dann tatsächlich eine Betaverteilung $B[\alpha, \beta]$,

\begin{displaymath}
p_{1}(x_{1}) = \int_{0}^{\infty} dx_{2} p(x_{1},x_{2}) =
\f...
...a(\alpha) \Gamma(\beta)} x_{1}^{\alpha-1}
(1-x_{1})^{\beta-1}.
\end{displaymath}

Dieses Ergebnis ist nicht selbsverständlich. Die Transformation von $(y_{1},y_{2})$ auf $(x_{1},x_{2})$ ist nicht umkehrbar eindeutig. Daher darf der von uns angewandte Transformationsalgorithmus eigentlich nicht verwendet werden. Wir erinnern uns aber der Aussage von Kap.3, wonach die Transformation trotzdem das richtige Ergebnis liefert, sofern das Resultat der Transformation eine auf 1 normierte Dichtefunktion ergibt. Dieses ist im vorliegenden Beispiel offensichtlich der Fall.

Der Algorithmus lautet: Erzeuge $y_{1} \in G[\alpha,1], y_{2} \in G[\beta,1]$. Dann ist

\begin{displaymath}
\xi = \frac{y_{1}}{y_{1}+y_{2}}
\end{displaymath}

eine Veränderliche der Betaverteilung $B[\alpha, \beta]$. Das Programm JBeta3 zeigt den Sourcecode.

Eine Reduktionsmethode mit Verwerfungsverfahren. Die folgende Methode beruht auf einem ähnlichen mathematischen Formalismus wie das soeben diskutierte Verfahren. Seien $y_{1},y_{2}$ zwei Veränderliche mit den Dichtefunktionen

$\displaystyle p_{1}(y_{1})$ $\textstyle =$ $\displaystyle \alpha y_{1}^{\alpha-1}, \; \; \; \; \; 0 \leq y_{1} \leq 1.$  
$\displaystyle p_{2}(y_{2})$ $\textstyle =$ $\displaystyle \beta y_{2}^{\beta-1}, \; \; \; \; \;
0 \leq y_{2} \leq 1.$  

Die gemeinsame Dichtefunktion ist

\begin{displaymath}
p(y_{1},y_{2}) = \alpha \beta y_{1}^{\alpha-1} y_{2}^{\beta-1}, \; \; \;
0 \leq y_{1} \leq 1, 0 \leq y_{2} \leq 1.
\end{displaymath}

Wir studieren wiederum die Veränderlichen
$\displaystyle x_{1}$ $\textstyle =$ $\displaystyle \frac{y_{1}}{y_{1}+y_{2}},$  
$\displaystyle x_{2}$ $\textstyle =$ $\displaystyle y_{1} + y_{2}.$  

Die transformierte Funktion in $x_{1}, x_{2}$,

\begin{displaymath}
p^{\ast}(x_{1},x_{2}) = \alpha \beta x_{1}^{\alpha-1} (1-x_{...
...a+\beta-1}, \; \; \; 0 \leq x_{1} \leq 1, 0 \leq x_{2} \leq 2,
\end{displaymath}

ist nun leider keine auf 1 normierte Dichtefunktion. Wie wir schon im vorherigen Beispiel diskutiert haben, ist die Transformation von $(y_{1},y_{2})$ auf $(x_{1},x_{2})$ nicht umkehrbar eindeutig. Eine genauere Analyse zeigt jedoch, daß die bedingte Wahrscheinlichkeitsdichte in $x_{1}$, unter der Voraussetzung, daß $0 \leq x_{2} \leq 1$, durch

\begin{displaymath}
p(x_{1}\vert 0 \leq x_{2} \leq 1) = \frac{\Gamma(\alpha+\bet...
...a-1} (1-x_{1})^{\beta-1},
\; \; \; \; \; 0 \leq x_{1} \leq 1,
\end{displaymath}

gegeben ist. Damit lautet der Algorithmus:
1. Erzeuge $u_{1},u_{2} \in U[0,1]$ und setze

$\displaystyle y_{1}$ $\textstyle =$ $\displaystyle u_{1}^{1/\alpha},$  
$\displaystyle y_{2}$ $\textstyle =$ $\displaystyle u_{2}^{1/\beta}.$  

Wenn $y_{1}+y_{2} > 1$, gehe zurück zu Punkt 1.
2. Die Größe

\begin{displaymath}
\xi = \frac{y_{1}}{y_{1}+y_{2}}
\end{displaymath}

ist eine Veränderliche aus $B[\alpha, \beta]$. JBeta4 zeigt das Programmlisting.

Von Neumannsches Verwerfungsverfahren. Die Betaverteilung ist auf das endliche Intervall von $x=0$ bis $x=1$ beschränkt. Daher kann man das von Neumannsche Verwerfungsverfahren anwenden. Das Maximum der Betaverteilung liegt an der Stelle

\begin{displaymath}
x_{max} = \frac{\alpha-1}{\alpha+\beta-2},
\end{displaymath}

und hat den Wert

\begin{displaymath}
p_{max} = p(x_{max}) =
\frac{\Gamma(\alpha+\beta)}{\Gamma(\...
...ha-1} \left(
\frac{\beta-1}{\alpha+\beta-2} \right)^{\beta-1}.
\end{displaymath}

Die Verwerfungsabfrage des von Neumannschen Verfahrens kann dann geschrieben werden als

\begin{displaymath}
u_{2} \leq \frac{p(\xi)}{p_{max}} = \left( \frac{\alpha+\bet...
...}{\beta-1} \right)^{\beta-1}
\xi^{\alpha-1} (1-\xi)^{\beta-1}.
\end{displaymath}

Zusammengefaßt lautet das Verfahren:
1. Berechne in einem externen Programm

\begin{displaymath}
\epsilon = \left( \frac{\alpha+\beta-2}{\alpha-1} \right)^{\alpha-1}
\left( \frac{\alpha+\beta-2}{\beta-1} \right)^{\beta-1}.
\end{displaymath}

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

\begin{displaymath}
u_{2} \leq \epsilon u_{1}^{\alpha-1} (1-u_{1})^{\beta-1},
\end{displaymath}

akzeptiere $\xi = u_{1}$ als Veränderliche von $B[\alpha, \beta]$, sonst beginne von vorne bei Punkt 2. JBeta5 zeigt das Listing.

Zwei weitere Verwerfungsmethoden. Da die beiden Funktionen $h_{\alpha}(x) = \alpha x^{\alpha-1}$ und $h_{\beta}(x) = \beta (1-x)^{\beta-1}$ für $0 \leq x \leq 1$ jeweils auf 1 normierte Dichtefunktionen sind, drängen sich zwei weitere Verwerfungsmethoden auf:
1. $\alpha$- Verfahren:

$\displaystyle h_{\alpha}(x)$ $\textstyle =$ $\displaystyle \alpha x^{\alpha-1},$  
$\displaystyle G_{\alpha}(x)$ $\textstyle =$ $\displaystyle (1-x)^{\beta-1},$  
$\displaystyle C_{\alpha}$ $\textstyle =$ $\displaystyle \frac{1}{\alpha} \frac{\Gamma(\alpha+\beta)}
{\Gamma(\alpha)\Gamma(\beta)}.$  

2. $\beta$- Verfahren:
$\displaystyle h_{\beta}(x)$ $\textstyle =$ $\displaystyle \beta(1-x)^{\beta-1},$  
$\displaystyle G_{\beta}(x)$ $\textstyle =$ $\displaystyle x^{\alpha-1},$  
$\displaystyle C_{\beta}$ $\textstyle =$ $\displaystyle \frac{1}{\beta} \frac{\Gamma(\alpha+\beta)}
{\Gamma(\alpha)\Gamma(\beta)}.$  

Vergleich der Erfolgswahrscheinlichkeiten $C_{\alpha}$ und $C_{\beta}$ zeigt, daß das $\alpha$- Verfahren für $\alpha > \beta$ und das $\beta$- Verfahren für $\beta > \alpha$ gewählt werden sollte. Beide Verfahren sind in JBeta6 vereinigt.

Cauchy-Verteilung
Die Cauchy- Verteilung $C[\alpha,\beta]$,

\begin{displaymath}
p(y) = \frac{1}{\pi \beta [1 + (\frac{y-\alpha}{\beta})^{2}]...
... \alpha < \infty, \; \beta > 0, \; -\infty \leq y \leq \infty,
\end{displaymath} (5)

ist eine auf der gesamten reellen Achse definierte Dichtefunktion. In physikalischen Anwendungen heißt diese Verteilung auch Breit- Wigner Verteilung. Die physikalische Bedeutung liegt im folgenden. Die Heisenbergsche Unschärferelation

\begin{displaymath}
\Delta E \cdot \Delta t \geq \hbar
\end{displaymath}

besagt, daß man die Energie $E$ und die Zeit $t$ nicht gleichzeitig genau messen kann. Setzt man für $\Delta t$ die Lebensdauer $\tau = 1/\Gamma$ eines sich in Ruhe befindliches quantenmechanischen Systems und für $E$ die Ruhemasse, so muß man folgern, daß die Masse eines instabilen Systems keinen scharf definierten Wert hat. Eine genauere Analyse zeigt, daß die Masse eine Veränderliche der Breit-Wigner Verteilung ist:
\begin{displaymath}
p(M) = \frac{\Gamma/2}{\pi [(M-M_{0})^{2} + (\Gamma/2)^{2}]}.
\end{displaymath} (6)

Dieses ist offensichtlich eine Cauchy- Verteilung mit $y = M, \alpha = M_{0}$ und $\beta = \Gamma/2$.

Die erzeugenden Funktionen sind analytisch nicht darstellbar. Das erste Moment $m_{1}$ kann man wegen der Symmetrie der Verteilung um $y = \alpha$ sofort hinschreiben,

\begin{displaymath}
m_{1} = <y> = \alpha.
\end{displaymath} (7)

Wie der Leser selbst feststellen wird, bereiten die höheren Momente erhebliche Schwierigkeiten. So ist z.B. das zweite zentrale Moment

\begin{displaymath}
\mu_{2} = <(y-\alpha)^{2}> = \lim_{y_{1} \to -\infty} \lim_{...
...\frac{\beta (y-\alpha)^{2}}{\pi [\beta^{2} +
(y-\alpha)^{2}]}
\end{displaymath}

nicht endlich integrierbar ( $\mu_{2} \to \infty$ für $y_{1} \to -\infty$ und $y_{2} \to \infty$). Wir verzichten daher auf eine weitere Diskussion.

Aus der Standard- Cauchy- Verteilung $C[0,1]$,

\begin{displaymath}
p(x) = \frac{1}{\pi [ 1 + x^{2}]}, \; \; \; \; \; -\infty \leq x \leq \infty,
\end{displaymath} (8)

kann die allgemeine Cauchy- Verteilung mit Hilfe der Transformation $y = \alpha + \beta x$ gewonnen werden. Es genügt daher, Simulationsalgorithmen für die Standard- Cauchy- Verteilung zu diskutieren.

Inverse Transformationsmethode. Die integrale Verteilungsfunktion der Standard- Cauchy- Verteilung $C[0,1]$,

\begin{displaymath}
P(\xi) = \int_{-\infty}^{\xi} dx p(x) = \frac{1}{2} + \frac{1}{\pi}
arctg(x),
\end{displaymath}

kann leicht invertiert werden. Mit $u \in U[0,1]$ ist dann

\begin{displaymath}
\xi = tg[ \pi (u - \frac{1}{2})] = \frac{1}{tg(\pi u)}
\end{displaymath}

eine Veränderliche aus $C[0,1]$. Das Programmlisting ist in JCauchy1 gezeigt.

Eine Reduktionsmethode. Ein zweites Verfahren benutzt eine Reduktion zweier normalverteilter Veränderlicher. Die zugehörige Mathematik hatten wir bereits diskutiert. Wir formulieren den Algorithmus: Erzeuge $v_{1},v_{2} \in N[0,1]$ und setze

\begin{displaymath}
\xi = \frac{v_{1}}{v_{2}}.
\end{displaymath}

$\xi$ ist dann eine Veränderliche aus $C[0,1]$. Das Programm ist in JCauchy2 gezeigt. Man beachte, daß in Routine JNormal1 zwei unabhängige normalverteilte Veränderliche erzeugt werden.

Ein Verwerfungsverfahren. Ein drittes Verfahren benutzt eine Verwerfung und lautet: Erzeuge $u_{1},u_{2} \in U[-\frac{1}{2}, +\frac{1}{2}]$ und setze

$\displaystyle x_{1}$ $\textstyle =$ $\displaystyle \frac{u_{1}}{u_{2}},$  
$\displaystyle x_{2}$ $\textstyle =$ $\displaystyle u_{1}^{2} + u_{2}^{2}.$  

Wenn dann $x_{2} \leq \frac{1}{4}$, ist $\xi = x_{1}$ eine Veränderliche aus $C[0,1]$. Sonst beginne von vorne. Wir überlassen es dem Leser, die Mathematik dieses Verfahrens zu überprüfen. Das Programmlisting zeigt JCauchy3.

Weibul-Verteilung
Die allgemeine Weibull- Verteilung $W[\alpha,\beta,\gamma]$,

\begin{displaymath}
p(y) = \frac{\alpha}{\beta^{\alpha}} (y-\gamma)^{\alpha-1}
e...
...^{\alpha}}, \; \; \; \alpha > 0, \; \beta > 0, \;
y > \gamma,
\end{displaymath} (9)

können wir mit Hilfe der Transformation $y = \gamma + \beta x$ aus der Standard- Weibul- Verteilung $W[\alpha, 1, 0]$,
\begin{displaymath}
p(x) = \alpha x^{\alpha-1} e^{-x^{\alpha}}, \; \; \; \; \; 0 \leq x,
\end{displaymath} (10)

erhalten. Die Weibul- Verteilung kann zur Darstellung von nichtresonanten effektiven Massen in der Teilchenphysik benutzt werden. Die Momente der Standard- Verteilung sind einfach zu berechnen und ergeben
\begin{displaymath}
<x^{j}> = \Gamma(1 + \frac{j}{\alpha}).
\end{displaymath} (11)

Daraus erhalten wir Erwartungswert und Varianz der allgemeinen Verteilung zu
$\displaystyle <y>$ $\textstyle =$ $\displaystyle \gamma + \beta \Gamma(1 + \frac{1}{\alpha}),$ (12)
$\displaystyle \sigma_{y}^{2} = <(y-<y>)^{2}>$ $\textstyle =$ $\displaystyle \beta^{2} \left( \Gamma(1+\frac{2}{\alpha})
- \Gamma^{2}(1+\frac{1}{\alpha}) \right).$ (13)

Ähnliche Ausdrücke können für die höheren Momente abgeleitet werden.

Inverse Transformationsmethode. Die Simulation erfolgt fast ausschließlich mit der inversen Transformationsmethode. Die integrale Verteilungsfunktion der Standard- Weibul- Verteilung,

\begin{displaymath}
P(\xi) = 1 - e^{-x^{\alpha}}
\end{displaymath}

kann leicht invertiert werden, sodaß mit $u \in U[0,1]$ die Größe

\begin{displaymath}
\xi = [-ln(u)]^{1/\alpha}
\end{displaymath}

eine Veränderliche aus $W[\alpha, 1, 0]$ ist. Das Programm zeigt JWeibul.

Verwerfungsverfahren. Natürlich gibt es auch hier eine Reihe von Verwerfungsverfahren. Man kann z.B. die beiden Faktoren $\alpha x^{\alpha-1}$ bzw $e^{-x^{\alpha}}$ als Dichtefunktion $h(x)$ und den jeweiligen anderen Faktor als Testfunktion $C \cdot G(x)$ ansetzen. Approximationen mit Hilfe der Gammafunktionen sind ebenfalls möglich. Alle diese Verwerfungsverfahren können aber nicht mit der in diesem Fall sehr einfachen inversen Transformationsmethode mithalten. Wir verzichten daher auf eine weitere Diskussion.

$\chi^{2}$-Verteilung
Die $\chi^{2}$- Verteilung $Chi[n]$ hatten wir bereits mehrfach diskutiert. Die Definition ist:

\begin{displaymath}
p(x) = \frac{1}{2^{n/2}\Gamma(\frac{n}{2})} x^{n/2 -1} e^{-x/2}.
\end{displaymath} (14)

Dieses ist eine spezielle Gammaverteilung $G[\alpha, \beta]$ mit $\alpha=n/2$ und $\beta=2$. Wir können also alle Algorithmen zur Simulation der Gammaverteilung auch zur Simulation der $\chi^{2}$- Verteilung benutzen. Insbesondere ist für geradzahliges $n$ die $\chi^{2}$- Verteilung identisch mit der Erlang- Verteilung $Er[\frac{n}{2}, 2]$. Für ungeradzahliges $n$ setzen wir

\begin{displaymath}
\frac{n}{2} = m + \frac{1}{2}
\end{displaymath}

mit geradzahligem $m$ und

\begin{displaymath}
x = y + z,
\end{displaymath}

wobei $y \in Er[m,2]$ und $z$ eine Veränderliche mit der Dichtefunktion

\begin{displaymath}
p(z) = \frac{1}{\Gamma(\frac{1}{2})} z^{-1/2} e^{-z/2}, \; \; \; \; \; z \geq 0,
\end{displaymath}

ist. Die Variablentransformation

\begin{displaymath}
w = \sqrt{z}, \; \; \; w \geq 0,
\end{displaymath}

führt auf die Dichtefunktion

\begin{displaymath}
p(w) = \sqrt{\frac{2}{\pi}} e^{-w^{2}/2}, \; \; \; w \geq 0,.
\end{displaymath}

Dieses ist nichts anderes als die rechte positive Hälfte der Standard- Normalverteilung. Zusammengefaßt erhalten wir den Algorithmus: Berechne $m$ und $\delta$ aus

\begin{displaymath}
\frac{n}{2} = m + \delta, \; \; \; \delta = 0 \; oder \; \frac{1}{2},
\end{displaymath}

und erzeuge $y \in Er[m]$. Falls $\delta = \frac{1}{2}$, erzeuge $w \in N[0,1]$. Die Größe

\begin{displaymath}
\xi = \{ \begin{array}{lll} y & f''ur & \delta = 0, \\
y + w^{2} & f''ur & \delta = \frac{1}{2}, \end{array}\end{displaymath}

ist dann eine Veränderliche der $\chi^{2}$- Verteilung $Chi[n]$. Das Verfahren zeigt JChisqr1.

Bei der $\chi^{2}$- Verteilung hat man es häufig mit großen Werten von $n$ zu tun. Hier empfiehlt sich die Approximation der Erlang- Verteilung durch eine Normalverteilung. In der Programmroutine JChisr1 wird in diesem Fall der Aufruf von JErlang1 durch JErlang2 ersetzt. Die Simulation der Standard- Normalverteilung JNormal3 kann natürlich durch jeden anderen Algorithmus ersetzt werden.

Faltung von $n$ Normalverteilungen. Ein zweites Verfahren zur Simulation der $\chi^{2}$- Verteilung folgt aus der schon mehrfach diskutierten Tatsache, daß die quadratische Summe von $n$ normalverteilten Veränderlichen $x_{\nu}, \nu=1,2,...,n$ einer $\chi^{2}$- Vertteilung gehorcht. Der Algorithmus lautet: Erzeuge $x_{1},x_{2},...,x_{n} \in N[0,1]$. Dann ist

\begin{displaymath}
\xi = \sum_{\nu=1}^{n} x_{\nu}^{2}
\end{displaymath}

eine Veränderliche aus $Chi[n]$. Das Programm ist in JChisqr2 gelistet.

Statistisches Wortmodell. Zu erwähnen ist ferner die Simulation der $\chi^{2}$- Verteilung, die wir bereits bei der Herleitung des $\chi^{2}$- Tests in Kap.1 benutzt hatten. Der Leser wird allerdings sofort bemerken, daß dieser Algorithmus mit den im vorliegenden Kapitel abgeleiteten Verfahren nicht mithalten kann.

Physikalische Bedeutung der Chisqr- Verteilung. Die Chisqr- Verteilung hat nicht nur in der Statistik eine überragende Bedeutung, sondern auch in physikalischen Anwendungen. Setzt man $n=3$ und führt die Transformation

\begin{displaymath}
x = \frac{m v^{2}}{kT}
\end{displaymath}

aus, so wird man auf die Maxwellsche Formel für die Geschwindigkeitsverteilung von Molekülen in einem idealen Gas geführt:

\begin{displaymath}
p(v) = 4 \pi \left( \frac{m}{2 \pi k T} \right)^{3/2} v^{2} e^{-m v^{2} /2kT}
\end{displaymath}

Diese Verteilung ist also eine Chisqr- Verteilung mit 3 Freiheitsgraden.

$t$-Verteilung
Wir betrachten zwei unabhängige Veränderliche $y_{1} \in N[0,1]$ und $y_{2} \in Chi[n]$, d.h.

$\displaystyle p_{1}(y_{1})$ $\textstyle =$ $\displaystyle \frac{1}{\sqrt{2 \pi}} e^{-y_{1}^{2}/2}, \; \; \; \; \;
-\infty < y_{1} < \infty,$  
$\displaystyle p_{2}(y_{2})$ $\textstyle =$ $\displaystyle \frac{1}{2^{n/2} \Gamma(n/2)} y_{2}^{n/2 -1}
e^{-y_{2}/2}, \; \; \; \; \; 0 \leq y_{2}.$  

Die gemeinsame Dichtefunktion ist:

\begin{displaymath}
p(y_{1},y_{2}) = \frac{1}{\sqrt{2\pi} 2^{n/2} \Gamma(n/2)}
y...
... y_{2}/2}, \; \; \; -\infty < y_{1} < \infty, \;
0 \leq y_{2}.
\end{displaymath}

Die Transformation
$\displaystyle x_{1}$ $\textstyle =$ $\displaystyle \frac{y_{1}}{\sqrt{y_{2}/n}},$  
$\displaystyle x_{2}$ $\textstyle =$ $\displaystyle \sqrt{y_{2}/n},$  

führt auf die Dichtefunktion in $x_{1}, x_{2}$:

\begin{displaymath}
p(x_{1},x_{2}) = \frac{n^{n/2}}{\sqrt{2\pi} 2^{n/2 -1} \Gamm...
..._{2}^{n} e^{-\frac{n}{2} (1 + \frac{x_{1}^{2}}{n}) x_{2}^{2}}.
\end{displaymath}

Integration über $x_{2}$ liefert die Randverteilung (mit $x \equiv x_{1}$):
\begin{displaymath}
p(x) = \frac{1}{\sqrt{n\pi}} \frac{\Gamma(\frac{n+1}{2})}{\G...
...}{(1+x^{2}/n)^{(n+1)/2}}, \; \; \; \; \; -\infty < x < \infty.
\end{displaymath} (15)

Diese Verteilung heißt t- Verteilung $T[n]$. Ein Simualtions- Algorithmus folgt sofort aus unserer Herleitung und ist als Unterroutine JTdist gelistet.

$F$-Verteilung
Als letztes diskutieren wir noch die $F$-Verteilung. Wenn $y_{1}$ eine Veränderliche einer $\chi^{2}$- Verteilung mit $n$ Freiheitsgraden und $y_{2}$ eine Veränderliche einer $\chi^{2}$- Verteilung mit $m$ Freiheitsgraden ist, dann ist die Variable

\begin{displaymath}
x = \frac{y_{1}/n}{y_{2}/m}
\end{displaymath}

eine Veränderliche einer $F$- Verteilung. Mit einer ähnlichen Rechnung wie vorher bei der $t$- Verteilung erhalten wir die Dichtefunktion $F[n,m]$.

\begin{displaymath}
p(x) = \left( \frac{n}{m} \right)^{n/2} \frac{\Gamma((n+m)/2...
...a(m/2)}
x^{n/2-1} \left( 1 + \frac{n}{m} x \right)^{-(n+m)/2}
\end{displaymath}

Der Code von diesem Algorithmus ist in JFdist gelistet.

Ausblick. Es gibt noch eine ganze Reihe weiterer kontinuierlicher Verteilungen, die in den Natur-, Ingenieur- und Sozialwissenschaften eine gewisse Rolle spielen. Falls der Leser bis hierhin den Ausführungen gefolgt ist, sollte es ihm keine Schwierigkeiten bereiten, auch hierfür einen geeigneten Simulationsalgorithmus zu entwickeln. Falls einem überhaupt nichts einfällt, kann man natürlich immer auf das von Neumann'sche Verwerfungsverfahren zurückgreifen, sofern man den Wertebereich der Variablen auf ein endliches Intervall einschränken kann. Dieses Verfahren habem wir in vorliegenden Abschnitt und im vorherigen Abschnitt bei der Diskussion der kontinuierlichen Verteilungen total ignoriert.




Harm Fesefeldt
2006-07-11