Simulation diskreter Veränderlicher

Vorbemerkung. Wir hatten die Simulation diskreter Veränderlicher schon im Zusammenhang mit der Simualtion kontinuierlicher Veränderlicher kurz angesprochen, und tatsächlich, die meisten Verfahren lassen sich ohne weiteres auf diskrete Veränderliche übertragen. Hinzu kommen aber einige Besonderheiten, die wir im folgenden Unterkapitel diskutieren wollen. Im Mittelpunkt unserer Betrachtungen wird die inverse Transformationsmethode stehen, die bei diskreten Veränderlichen immer anwendbar ist. Verwerfungsmethoden spielen zwar eine weniger wichtige Rolle, jedoch gibt es auch bei diskreten Veränderlichen einige schöne Anwendungsbeispiele. Hinzu kommen dagegen Methoden der Diskretisierung kontinuierlicher Veränderlicher und insbesondere Approximationen durch kontinuierliche Funktionen. Auch die Simulation des verbalen Wortmodells führt bei diskreten Veränderlichen zu außerordentlich schönen und schnellen Algorithmen.

Inverse Transformationsmethode. Die Wahrscheinlichkeitsverteilung $P(x_{k})$ einer diskreten Veränderlichen $x_{k}, -\infty < k < \infty$, mit ganzzahligem $k$, kann in praktischen Anwendungen immer auf eine Wahrscheinlichkeitsverteilung $P(k)$ mit einer positiven ganzzahligen Veränderlichen $k$ transformiert werden. Wir legen daher im folgenden immer eine Verteilung $P(k)$ mit $k \geq 0$ zugrunde. Die integrale Verteilung ist

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

und die inverse Transformationsmethode lautet: Erzeuge $u \in U[0,1]$ und bestimme $\kappa$ so, daß

$$\sum_{k=0}^{\kappa -1} P(k) < u \leq \sum_{k=0}^{\kappa} P(k).$$ (2)

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

Wir haben in unseren bisherigen Betrachtungen jedoch nicht über optimale Algorithmen zur Lösung der Ungleichung (2) gesprochen. Das einfachste Verfahren basiert natürlich auf den Rekursionsrelationen

$$P(k+1) = a'_{k+1} P(k)$$

$$P(k-1) = a''_{k-1} P(k)$$

wobei $a'_{k+1}$ und $a''_{k-1}$ im allgemeinen einfache Funktionen von $k$ sind. Mit Hilfe dieser Funktion $a'_{k+1}$ kann man das Verfahren folgendermaßen formulieren:
1. Erzeuge $u \in U[0,1]$, berechne $P_{0} = P(k=0)$ und setze

$$c = P_{0},$$

$$b = c,$$

$$k = 0.$$

2. Wenn $u \leq b$, dann ist $k$ eine Veränderliche der Wahrscheinlichkeitsverteilung $P(k)$. Wenn dagegen $u < b$, dann berechne $a'_{k+1}$ und setze

$$c \gets a'_{k+1} c,$$

$$b \gets b + c,$$

$$k \gets k + 1,$$

und gehe zurück zu Punkt 2.

Dieser Algorithmus ist natürlich nur ein Spezialfall des allgemeineren Verfahrens, bei dem man bei irdendeinem Wert $k = m$ mit den Vergleichsoperationen beginnt. Das Verfahren lautet dann entsprechend:
1. Erzeuge $u \in U[0,1]$ und setze

$$d = \sum_{k=0}^{m} P(k),$$

$$e = P(m),$$

$$k = m.$$

Falls $u > d$, fahre fort bei Punkt 3, sonst bei Punkt 2.
2. Setze

$$d \gets d - e.$$

Falls $u > d$, so ist $k$ eine Veränderliche der Wahrscheinlichkeitsverteilung $P(k)$. Wenn dagegen $u \leq d$, dann berechne $a''_{k-1}$ und setze

$$e \gets a''_{k-1} e,$$

$$k \gets k - 1,$$

und gehe zurück zu Punkt 2.
3. Berechne $a'_{k+1}$ und setze

$$e \gets a'_{k+1} e,$$

$$d \gets d + e,$$

$$k \gets k + 1.$$

Falls $u \leq d$, so ist $k$ eine Veränderliche von $P(k)$, sonst gehe zurück zu Punkt 3.

Entscheidend für die Schnelligkeit dieses allgemeineren Verfahrens ist natürlich einmal die Berechnung des Startwertes für $d = \sum_{k=0}^{m} P(k)$. Hierfür kann man im allgemeinen keinen analytischen Ausdruck angeben, sodaß die Summe tatsächlich aus $m$ arithmetischen Operationen berechnet werden muß. Dieses kann die Vorteile, die im folgenden diskutiert werden, durchaus wieder zunichte machen. Zweitens hängt die Schnelligkeit des Verfahrens von der Anzahl der Vergleichsoperationen ab. Diese Anzahl ist offensichtlich gegeben durch

$$\eta = \{ \begin{array}{lll} 2 + (m-k) & f''ur & k \leq m, \\ 1 + (k-m) & f''ur & k > m. \end{array}$$

gegeben. Der Erwartungswert von $\eta$ ist

$$C = <\eta> = \sum_{k=0}^{m} [2+(m-k)] P(k) - \sum_{k=m+1}^{\infty} [1+(k-m)] P(k) = 1 + <k> - \gamma(m),$$

mit

$$\gamma(m) = m + \sum_{k=0}^{m} (2k - 2m - 1) P(k).$$

Ähnlich wie bei den Verwerfungsverfahren nennen wir die Größe $1/C$ die Erfolgswahrscheinlichkeit:

$$\frac{1}{C} = \frac{1}{1+<k>-\gamma(m)}.$$ (3)

Wie man sich leicht klarmachen kann, kann $\gamma(m)$ positive wie negative Werte annehmen. Den optimalen Wert für $m$ erhält man aus dem folgenden Satz:
Sei $g(m)$ und der Median $m_{med}$ definiert durch

$$g(m) = P(m+1) + 2 \sum_{k=0}^{m} P(k),$$

$$m_{med} = max \{ m ; \sum_{k=0}^{m} p(k) \leq \frac{1}{2} \}.$$

Dann ist der optimale Startwert $m_{0}$ durch

$$m_{0} = \{ \begin{array}{lll} m_{med} & f''ur & g(m_{0}) \leq 1 \\ m_{med} + 1 & f''ur & g(m_{0}) > 1. \end{array}$$

Vorraussetzung für beide Fälle ist jedoch, daß

$$P_{0} = P(k=0) \leq \sum_{k=1}^{m_{0}+1} (2k-1) P(k).$$

In allen Programmbeispielen werden wir mit $m_{med}$ als Parameter beginnen, dann ist der Startwert $d = 1/2$. Es muß aber betont werden, daß dieses Verfahren nicht notwendig das optimale Verfahren ist.

Verwerfungsverfahren. Die Verwerfungsverfahren werden ganz analog wie bei kontinuierlichen Veränderlichen eingeführt. Wir stellen die Wahrscheinlichkeitsverteilung in der Form

$$P(k) = C H(k) G(k)$$ (4)

mit

$$\sum_{k=0}^{\infty} H(k) = 1,$$

$$max G(k) = 1,$$

$$C \geq 1,$$

dar. Hierbei sollte $H(k)$ eine möglichst gute Approximation von $P(k)$ sein und wir sollten für $H(k)$ bereits einen Simulations- Algorithmus besitzen. Das Verfahren lautet also: Erzeuge $u \in U[0,1]$ und eine Veränderliche $k$ aus der Wahrscheinlichkeitsverteilung $H(k)$. Wenn dann

$$u \leq G(k),$$

akzeptiere $k$ als Veränderliche der Wahrscheinlichkeitsverteilung $P(k)$, sonst beginne von vorne.

Besonders einfach ist bei diskreten Veränderlichen das von Neumannsche Verfahren, sofern die Veränderliche $k$ auf endlich viele Werte beschränkt ist, d.h. $k=0,1,2,...,n$. Für $H(k)$ setzen wir die diskrete Gleichverteilung

$$H(k) = \frac{1}{n+1}, \; \; \; \; \; k=0,1,2,...,n,$$ (5)

und für $G(k)$ die Funktion

$$G(k) = \frac{P(k)}{P_{max}}.$$ (6)

Das Verfahren lautet: Erzeuge $u_{1},u_{2} \in U[0,1]$ und bestimme $k$ gemäß

$$(n+1) u_{1} = k + \delta, \; \; \; \; \; 0 < \delta \leq 1.$$

Wenn dann

$$u_{2} \leq \frac{P(k)}{P_{max}},$$

akzeptiere $k$ als Veränderliche der Wahrscheinlichkeitsverteilung $P(k)$, sonst beginne von vorne.

Diskretisierung kontinuierlicher Funktionen. Jede Wahrscheinlichkeitsverteilung kann als Integral über eine kontinuierliche Dichtefunktion dargestellt werden:

$$P(k) = \int_{k}^{k+1} dx p(x),$$ (7)

mit

$$\sum_{k=0}^{\infty} P(k) = \int_{0}^{\infty} dx \; p(x) = 1.$$

Für große Werte von $k$ kann das Inegral approximiert werden durch (man beachte: $\Delta k = 1$)

$$P(k) \approx p(x=k+\frac{1}{2}) \Delta k = p(x=k+\frac{1}{2}).$$ (8)

Das Verfahren lautet: Erzeuge eine Veränderliche $\xi$ aus der Dichtefunktion $p(x)$ und bestimme $k$ gemäß

$$\xi = k + \delta, \; \; \; \; \; 0 < \delta < 1.$$

$k$ ist dann eine Veränderliche aus $P(k)$.

Simulation des verbalen Wortmodells. Bei diskreten Veränderlichen liefert die Übertragung des verbalen Wortmodells in ein Rechnerprogramm häufig sehr schnelle Algorithmen. Ein Beispiel möge dieses belegen. Eine Veränderliche der Binominal- oder Bernoulli- Verteilung hatten wir definiert als Anzahl der positiven Versuche in einer Reihe von insgesamt $n$ Versuchen, wobei der Einzelversuch mit einer Wahrscheinlichkeit $p$ zu einem positiven Ergebnis führt. Wir simulieren den Einzelversuch durch den einmaligen Aufruf des Zufallszahlen- Generators. Der Versuch wird als positiv bewertet, sofern die Zufallszahl $u \leq p$. Diesen Einzelversuch wiederholen wir $n$-mal und zählen die Anzahl $k$ der positiven Ereignisse. Die Zahl $k$ ist dann offenbar eine Veränderliche der Binominalverteiling $B[p,n]$. Der Rechenaufwand in diesem Algorithmus besteht im wesentlichen aus der Erzeugung von $n$ Zufallszahlen aus $U[0,1]$. Ähnliche Wortmodelle können für fast alle diskreten Veränderlichen formuliert werden.

Bevor wir fortfahren, sollten Sie das Applet über die diskreten Verteilungen auf Ihren Bildschirm laden.

Binomial-Verteilung
Wir bezeichnen die Binomialverteilung

$$P(k) = \left( \begin{array}{c} n \\ k \end{array} \right) p^{k} (1-p)^{n-k}$$

mit $B(n,p)$. Die Methode JBinomial1 beschreibt das Wortmodell. Wir erzeugen uns $n$ Zufallszahlen $u_{i}$ mit dem Gleichverteilungsgenerator Ranmar und fragen nach der Anzahl der Zahlen $k$ mit $u_{i} < p$. In der Methode JBinomial2 haben wir die inverse Transformationsmethode realisiert. Hierfür wurden die folgenden Ausdrücke berechnet:

$$m_{med} = [(n+1)p]$$

$$a_{k+1}' = \frac{(n-k)}{(k+1)} \frac{p}{(1-p)}$$

$$a_{k-1}'' = \frac{k}{(n-k+1)} \frac{(1-p)}{p}$$

Die eckigen Klammer $[...]$ bedeuten hier und im folgenden die Abrundung auf eine ganze Zahl. Die Anzahl der Vergleichsoperationen ist im Mittel durch

$$C = 1.5 + \sqrt{\frac{2}{\pi}} \sqrt{np(1-p)}$$

gegeben. Die Verteilung dieser Zahlen zeigen wir im Applet im Plot RandomCalls.

Für große Werte von $n$ kann die Binomial- Verteilung durch eine Normalverteilung approximiert werden (siehe auch Kap.1 dieses Tutorials). Die Größe

$$z = \frac{k - np + 0.5}{\sqrt{np(1-p)}}$$ (9)

ist näherungsweise normal verteilt, mit einer Standardabweichung von $\sigma = 1$. In der Methode JBinomial3 erzeugen wir uns daher eine normal verteilte Variable $z$ aus $N[0,1]$ und lösen die Gleichung (9) nach $k$ auf:

$$k = [-0.5 + np + z \sqrt{np(1-p)}]$$

Poisson- Verteilung
Wir schreiben die Poisson- Verteilung in der Form

$$P(k) = \frac{1}{k!} \lambda^{k} e^{-\lambda}, \; \; \; \; \; k \geq 0, \; \; \; \lambda > 0.$$

und bezeichnen diese mit $P[\lambda]$. Unser erstes Modell kommt aus dem radioaktiven Zerfall. Wenn das Zeitintervall $T$ zwischen zwei Zerfällen aus einer Exponentialverteilung $E[1/\lambda]$ ist, dann ist die Anzahl der Zerfälle in der Zeiteinheit durch eine Poissonverteilung $P[\lambda]$ gegeben. Mathematisch können wir diesen Sachverhalt formulieren als

$$\sum_{i=0}^{k} T_{i} \leq 1 \leq \sum_{i=0}^{k+1} T_{i}$$

mit $T_{i} = -(1/\lambda) ln(u_{i})$ aus $E[1/\lambda]$ und $u_{i}$ aus $U[0,1]$. Daher gilt auch

$$- \sum_{i=0}^{k} ln(u_{i}) \leq \lambda \leq - \sum_{i=0}^{k+1} ln(u_{i})$$

oder

$$\prod_{i=0}^{k} u_{i} \geq e^{-\lambda} \geq \prod_{i=0}^{k+1} u_{i}$$

Der in JPoisson1 entwickelte Algorithmus basiert auf dieser letzten Ungleichung.

In JPoisson2 wird wiederum die inverse Transformationsmethode angewendet. Die entsprechenden Formeln sind

$$m_{med} = [\lambda]$$

$$a_{k+1}' = \frac{\lambda}{k+1}$$

$$a_{k-1}'' = \frac{k}{\lambda}$$

Die Anzahl der notwendigen Vergleichsoperationen ist ähnlich wie bei der Binomial- Verteilung durch

$$C = 1.5 + \sqrt{\frac{2}{\pi}} \sqrt{\lambda}$$

gegeben. Auch die Poisson- Verteilung kann für große Werte von $\lambda$ durch eine Normalverteilung approximiert werden. Die Variable

$$z = \frac{k - \lambda + 0.5}{\sqrt{\lambda}}$$

ist dann aus $N[0,1]$. Wir lösen wieder nach $k$ auf und erhalten

$$k = [\lambda + z \sqrt{\lambda} - 0.5].$$

JPoisson3 zeigt diesen Algorithmus.

Geometrische Verteilung
Die geometrische Verteilung

$$P(k) = p (1-p)^{k}, \; \; \; \; \; k=0,1,2,.... \; \; \; \; \; 0<p<1,$$

wird mit $Ge[p]$ abgekürzt und beschreibt die Anzahl der Versuche in einem Bernoulli- Experiment bis zum ersten erfogreichen Versuch. Dieses Wortmodell ist in JGeometric1 in einen Simulations- Algorithmus umgesetzt worden. In der inversen Transformationsmethode in JGeometric2 starten wir die Vergleichsoperationen ausnahmsweise bei $P_{0} = P(0) = p$ und benutzen $a'_{k+1} = 1-p$. Die mittlere Zahl der Vergleichsoperationen ist dann

$$C = 1 + \frac{1-p}{p}.$$

Ein sehr schöner Algorithmus folgt aus der Verbindung der geometrischen Verteilung und der Exponentialverteilung. Wenn $y$ eine Variable aus $E[\beta]$ ist, dann gilt

$$Pr(k < y < k+1) = \frac{1}{\beta} \int_{k}^{k+1} dy \; e^{-y/\beta} = e^{-y/\beta} ( 1 - e^{-y/\beta})$$

Der rechts stehende Ausdruck ist aber nichts anderes als $Ge[p = 1-e^{-y/\beta}]$. Wir erzeugen also eine Variable $u$ aus einer Gleichverteilung und berechnen

$$k = \left[ \frac{ln(u)}{ln(1-p)} \right].$$

$k$ ist dann eine Varible aus der geometrischen Verteilung. Der Algorithmus ist in JGeometric3 gezeigt.

Negative Binomialverteilung
Als letztes diskutieren wir noch die negative Binomialverteilung $NB[r,p]$:

$$P(k) = \left( \begin{array}{c} k + r -1 \\ k \end{array} \right) p^{r} (1-p)^{k}.$$

Hier ist $k$ die Anzahl der positiven Versuche in einem Bernoulli- Experiment bis zum Auftreten von $r$ negativen Versuchen. Dieses ist exakt das Wortmodell in JNegBinomial1. Für die inverse Transformationsmethode in JNegBinomial2 haben wir die folgenden Ausdrücke berechnet:

$$m_{med} = \left[ \frac{(r-p)(1-p)}{p} \right]$$

$$a'_{k+1} = \frac{(r+k)(1-p)}{k+1}$$

$$a''_{k-1} = \frac{k}{(r+k-1)(1-p)}$$

Die mittlere Anzahl der Vergleichsoperationen ist

$$C = 1 + \frac{r(1-p)}{p}.$$

Der dritte Algorithmus ist hergeleitet aus der Verbindung der negativen Binomialverteilung mit der Gammaverteilung und der Poissonverteilung. Wir nehmen an, daß der Parameter $\lambda$ in der Poissonverteilung

$$P(k\vert\lambda) = \frac{1}{k!} e^{-\lambda} \lambda^{k}$$

mit einer Gammaverteilung variiert:

$$f(\lambda) = \frac{1}{\beta^{\alpha} \Gamma(\alpha)} \lambda^{\alpha-1} e^{-\lambda/\beta}$$

Dann erhalten wir die gefaltete Verteilung

$$P(k) = \int_{0}^{\infty} P(k\vert\lambda) f(\lambda) \; d\lambda.$$

Die Berechnung führt auf

$$P(k) = \left( \begin{array}{c} \alpha + k -1 \\ k \end{array... ...beta+1} \right)^{k} \left( \frac{1}{1+\beta} \right)^{\alpha}.$$

Also ist $k$ eine Variable aus $NB[\alpha,1/(1+\beta)]$. Unser Algorithmus lautet also: Erzeuge $\lambda$ aus Gammaverteilung $G[r,p/(1-p)]$ und $k$ aus der Poissonverteilung $P[\lambda]$. Der Source- Code wird in JNegBinomial3 gezeigt. Der letzte Algorithmus in JNegBinomial4 benutzt die Verbindung zwischen der negativen Binomialverteilung und der geometrischen Verteilung. Der Algorithmus lautet: Erzeuge $r$ Variable $y_{1}, y_{2}, ...., y_{r}$ aus der geometrischen Verteilung $Ge[p]$ und setze $k = \sum_{i=1}^{r} y_{i}$.

Schlußbemerkung
Hiermit wollen wir unser Tutorial über Statistik mit Schwerpunkt auf Simulationsalgorithmen zufälliger Veränderlicher beschliessen. Der aufmerksame Leser sollte hiermit in der Lage sein, auch für andere Veränderliche effektive Algorithmen zu entwickeln. Die Simulation ist heute eine wichtige Methode in allen Gebieten der Wissenschaften geworden.