Sei \(f(x)>0\) eine beliebige stetige Funktion mit bekanntem Wertebereich \([0,Y]\). Dann kann \(\int_a^b f(x) dx\) wie folgt approximiert werden:

• Ziehe \(n\) gleichverteilte Zufallszahlen \(x\) aus \([a,b]\)
• Ziehe unabhängig davon \(n\) gleichverteilte Zufallszahlen \(y\) aus \([0,Y]\)
• Berechne den Anteil \(h\) der Punkte \((x_i,y_i)\), die unterhalb der Funktion \(f\) liegen
• \(\int_a^b f(x) dx \approx h(b-a)Y\)

Beispiel

Als Beispiel für die Monte Carlo-Integration approximieren wir \(\int_0^1 x^2 dx\).

Analytisch ergibt sich

\[ \int_0^1 x^2 dx = \left[\frac{1}{3}x^3\right]_0^1=\frac{1}{3} \]

Hier können Sie Punkte aus dem Einheitsintervall ziehen. Der Anteil der Punkte unterhalb (grün) der Funktion \(f(x)=x^2\) schätzt die Fläche und damit das Integral.

Monte-Carlo-Schätzer

Ist \(p(x)\) eine Dichte, so können Integrale der Form

\[ E(g(x))=\int g(x)p(x) dx \]

mit einer Stichprobe \(x_1,\ldots,x_m\) aus \(p(x)\) durch den Stichprobenmittelwert

\[ \bar{g}_m=\frac{1}{m}\sum_{i=1}^mg\left(x_i\right) \]

approximiert werden. Aus dem starken Gesetz der grossen Zahlen folgt \[ \lim \frac{1}{m} \sum_{i=1}^m g(x_i)\to \int g(x)p(x)dx. \]

• Das Gesetz der großen Zahlen gibt es auch für nicht unabhängige Zufallszahlen (Ergodensatz)

Varianz des Monte-Carlo-Schätzer

Die Varianz des Monte-Carlo-Schätzers \(\bar{g}_m\) ist gegeben durch

\[ Var(\bar{g}_m)=\frac{1}{m}\int\left(g(X)-E(g(x))\right)^2 p(x)dx=\frac{1}{m}Var(g). \]

Der Approximationsfehler verringert sich mit steigendem \(m\). Es folgt sogar aus dem zentralen Grenzwertsatz

\[ \sqrt{m}\left(g_m-E(g(x))\right)\sim N(0,Var(g)). \]

Schätzer für \(Var(\bar{g}_m)\) ist

\[ \widehat{Var(\bar{g}_m)}=\frac{1}{m-1}\sum_{i=1}^m \left(g(x_i)-\bar{g}_m\right)^2 \]

Inversionsmethode

Gegeben sei die Verteilungsfunktion \(F(x)\) einer Zufallsvariablen \(X\). Sei \(u\sim U[0,1]\). Dann ist

\[ Y = F^{-1}(u) = \inf\{y:F(y)\geq u\} \sim X \]

y <- runif(10000)

lambda <- 2
x <- -log(1-y)/lambda

Acception-Rejection-Methode

Ziel: Wir wollen aus einer Dichtefunktion \(f(x)\) ziehen. Gegeben sei eine Dichtefunktion \(g(x)\), nach der wir problemlos Zufallszahlen ziehen können. Es existiere ein \(c\), so dass

\[ cg(z)\geq f(z) \]

für alle \(z\) mit \(f(z)>0\). Dann können Zufallszahlen gemäß \(f(x)\) wie folgt gezogen werden:

• ziehe \(Z\) gemäß \(g(z)\)

• akzeptiere \(Z\) mit Wahrscheinlichkeit \(\frac{f(z)}{g(z)c}\)

• Wähle \(c\) möglichst klein.

• Kann auch angewandt werden, falls die Normalisierungskonstante von \(f(x)\) nicht bekannt.

Squezed Acception-Rejection-Sampling

Gegeben sei eine untere Schranke \(s(z) \leq f(z)\).

Für \(u\) Ziehung aus \(U[0,1]\) akzeptiere \(Z\)

• wenn \(u\leq \frac{s(z)}{cg(z)}\)
• wenn \(u\leq \frac{f(z)}{cg(z)}\)

Der zweite Schritt kann ausgelassen werden, wenn bereits im ersten Schritt akzeptiert wurde.

Importance Sampling

• Idee zum bestimmen des Erwartungswerts einer Verteilung mit Dichte \(f(\theta)\)
• Ziehe aus einer Importance Dichte \(q\), die etwa der gewünschten Dichte entspricht.
• Benutze dann folgenden Zusammenhang

\[ E(g(\theta))=\int \frac{g(\theta)f(\theta)}{q(\theta)} q(\theta) d\theta = E_q\left(\frac{g(\theta)f(\theta)}{q(\theta)}\right) \]

Mit Ziehungen \(\theta_1,\ldots,\theta_m\) gemäß \(q\) gilt

\[ \hat{g}_{IS}=\frac{1}{m}\sum_{i=1}^m g(\theta_i)\frac{f(\theta_i)}{q(\theta_i)} \]

wobei \(\frac{p(\theta_i|x)}{q(\theta_i)}\) die sogenannten Importancegewichte sind.

Die Berechnung des Importanceschätzers benötigt die normierte Dichte. Die Normierungskonstante kann zuvor ebenfalls mit Importance Sampling geschätzt werden.

Die Varianz des Schätzers ist

\[ Var(\hat{g}_{IS}) = \frac{1}{m}Var_q\left(\frac{g(\theta)p(\theta|x)}{q(\theta)}\right) \]

Markoveigenschaft

Eine Kette von Zufallszahlen \(Y=\{Y_t,t\in \mathbb{N}_0\}\) mit abzählbarem Zustandsraum \(S\) heisst Markov-Kette, wenn

\[ P(Y_t=k|Y_0=j_0,Y_1=j_1,\ldots,Y_{t-1}=j_{t-1})= P(Y_t=k|Y_{t-1}=j_{t-1}) \]

für alle \(t \geq 0\) und für alle \(k,j_0,\ldots,j_{t-1} \in S\). Das heißt, wenn der nächste Schritt nur vom momentanen Zustand abhängt.

Definitionen

• \(P(Y_t=k|Y_{t-1}=j)\) heisst Übergangswahrscheinlichkeit
• die Markovkette ist homogen, wenn \(P(Y_t=k|Y_{t-1}=j)\) nicht von \(t\) abhängt. Dann ist

\[ p_{jk}=P(Y_t=k|Y_{t-1}=j)=P(Y_1=k|Y_0=j) \]

• Die Matrix \(\mathbf{P}=(p_{jk})\) heisst Übergangsmatrix
• Eine Markovkette heißt irreduzibel, falls für alle \(j,k\in S\) eine positive Zahl \(1\leq t \leq \infty\) existiert, so dass

\[ P(Y_t=k|Y_0=j)>0, \]

also jeder Zustand \(k\) von jedem Zustand \(j\) in endlicher Zeit erreicht werden kann * Die Periode eines Zustands \(k\) ist der größte gemeinsame Teiler der Zeitpunkte \(n\), zu denen eine Rückkehr möglich ist. Falls alle Zustände einer Markov-Kette die Periode 1 haben, ist die Kette aperiodisch.

Beispiel:

Gegeben sei die Übergangsmatrix

\[ P=\left(\begin{array}{ccc} 0.1 & 0.7 & 0.2 \\ 0.1 & 0.5 & 0.4 \\ 0.1 & 0.45 & 0.45 \end{array} \right) \]

• Das heißt, ist man im Zustand 1, springt man mit Wahrscheinlichkeit 0.7 in den Zustand 2, mit Wahrscheinlichkeit 0.2 in den Zustand 3 und mit Wahrscheinlichkeit 0.1 bleibt man im Zustand 1.
• Man kann leicht nachrechnen, dass mit \(\pi=(0.1, 0.5, 0.4)\) gilt: \(\pi=\pi P\).
• Erzeugen wir nun eine Markovkette mit diesen Übergangswahrscheinlichkeiten.
• Wir starten im Zustand 1 (für viele Iterationen ist der Startzustand jedoch irrelevant).
• Ist die Kette lang genug, dann sind die relativen Häufigkeiten der Zustände gleich der invarianten Verteilung \(\pi\).

Invariante Verteilung

• Statt eines Startzustandes kann man sich auch eine Startverteilung \(\pi_0\) vorgeben (man startet zufällig in einem Zustand)

• Eine Markovkette heisst ergodisch, wenn die Zustandsverteilung \(\pi_t\) von \(Y_t\) für jede beliebige Startverteilung \(\pi_0\) gegen die selbe Wahrscheinlichkeitsverteilung \(\pi\) konvergiert:

  \[ \lim_{t \to \infty} \pi_t = \lim_{t \to \infty}\pi_0 \mathbf{P}^t = \pi \]

  • Die Grenzverteilung einer ergodischen Markovkette ist die invariante Verteilung \(\pi\)
  • Eine homogene Markovkette mit Übergangsmatrix \(\mathbf{P}\) ist ergodisch, wenn sie irreduzibel und aperiodisch ist
  • Die Zustandsverteilung einer irreduziblen und aperiodischen, homogenen Markovkette kovergiert daher gegen die stationäre Verteilung \(\pi\)
  • Eine ergodische Markovkette wird asymptotisch stationär, d.h., der Einfluss der Startverteilung geht verloren

• Hier haben wir nur diskrete Zustände, die Definitionen lassen sich auf stetige Zustände übertragen (statt von der Übergangsmatrix spricht man dann vom Übergangskern)
• Beim MCMC-Algorithmus müssen wir also eine homogene, irreduzible und aperiodische Markovkette erzeugen.
• Dies gewährleisten die bekannten MCMC-Methoden.

Vorschlagsdichten

• Independence Proposal: Vorschlagsverteilung ist unabhängig vom aktuellen Wert

• Symetrisches Proposal: \(q(\theta^*|\theta^{(k-1)})=q(\theta^{(k-1)}|\theta^*)\), die Vorschlagsdichte kürzt sich aus der Akzeptanzwahrscheinlichkeit (Metropolis-Algorithmus):

  \[ \alpha=\frac{p(\theta^*|x)}{p(\theta^{(k-1)}|x)} \]

  Jeder Vorschlag mit \(p(\theta^*|x)>p(\theta^{(k-1)}|x)\) wird angenommen!

  • Random Walk Proposal: Vorschlagsverteilung ist ein Random Walk

\[ \theta^* = \theta^{(k-1)}+\epsilon, \epsilon\sim f \]

also \(q(\theta^*|\theta^{(k-1)})=f(\theta^*-\theta^{(k-1)})\).

Random Walk

Random Walk wird in der Regel mit Normalverteilung konstruiert:

\[ \theta^* \sim N(\theta^{(k-1)},C) \]

mit vorgegebener Kovarianzmatrix \(C\).

• Eine zu kleine Varianz führt zu hohen Akzeptanzraten, aber ineffizienten, da stark autokorrelierten Ziehungen. Im Extremfall \(C\to 0\) führt zu \(\alpha = 1\), \(\tau \to \infty\).
• Eine zu große Varianz führt zu zu großen Schritten, Vorschläge liegen in den Enden der Posterioriverteilung, sehr kleine Akzeptanzraten.
• Tuning der Kovarianzmatrix notwendig

Multivariate Verteilung

• Metropolis-Hastings-Algorithmus kann für \(\theta\)-Vektoren durchgeführt werden

• Akzeptanzraten jedoch i.d.R. geringer mit höherer Dimension

• Alternative ist der komponentenweise Metropolis-Hastings: Jede Komponente des Parameters wird einzeln (skalar oder blockweise) aufdatiert. Sei \(\theta=(\theta_1,\theta_2)\):

  \[ \alpha=\min\left(1,\frac{p(\theta_1^*|x,\theta_2^{(k-1)})q(\theta^{(k-1)}_1|\theta_1^*)}{p(\theta^{(k-1)}_1|x,\theta_2^{(k-1)})q(\theta^*_1|\theta_1^{(k-1)})}\right) \]

  • Updates können in fester oder zufälliger Weise erfolgen

Gibbs-Block-Algorithmus

Allgemeiner können die Komponenten des Parametervektors in Blöcke zerlegt: \(\mathbf\theta=(\mathbf\theta_1,\ldots,\mathbf\theta_p)\). Jeder Block wird dann aus der vollständig bedingten Dichte \(p(\mathbf\theta_j|\mathbf{x},\mathbf\theta_{-j})\) gezogen.