Wir kennnen nun zwei MCMC-Algorithmen:

• Gibbs-Sampler
• Metropolis-Algorithmus

MCMC steht für Markov Chain Monte Carlo.

• Monte Carlo bezeichnet alle numerischen Algorithmen, die auf dem Ziehen von Zufallsziehen beruhen (nach dem bekannten Spielkasino in Monaco).
• Markov Chain, zu Deutsch Markov-Kette bezieht sich aus die Markov-Eigenschaft der Ziehungen:

Jede Ziehung \(k\) ist von der vorherigen Ziehung \(k-1\) abhängig. Gegeben der vorherigen Ziehung \(k-1\), ist die Ziehung \(k\) aber unabhängig von der Ziehung \(k-2\) und allen Ziehungen davor (\(k-3, k-4, \ldots\)).

Man spricht für die Menge der Ziehungen daher auch von der Kette bzw. chain

• Beim Gibbs-Sampler ziehen wir jeden Parameter einzeln abwechselnd aus den full conditionals. Das geht natürlich nur, wenn wir daraus ziehen können (z.B. wenn diese Standardverteilungen sind)
• Beim Metropolis-Algorithmus ziehen wir einen Vorschlag für alle Parameter gemeinsam und akzeptieren diesen dann mit einer Akzeptanzwahrscheinlichkeit. Bei der Berechnung der Akzeptanzwahrscheinlichkeit reicht es, die Posteriori-Dichte bis auf Konstanten zu wissen. Nachteil des Metropolis-Algorithmus ist es, dass wir eine Vorschlagsdichte wählen müssen, deren Wahl die Effizienz des Algorithmus beeinflußt.

Eine Kombination der beiden Algorithmen ist der Metropolis- within-Gibbs-Algorithmus:

• Ziehe wie beim Gibbs-Sampler aus \(\theta_i|\theta_{-i}\)
• Falls die Verteilung von \(\theta_i|\theta_{-i}\) nur bis auf eine Proportionalitätskonstante bekannt ist, benutze den Metropolis-Algorithmus nur für diesen Schritt

Als Beispiel erweiteren wir unser Poisson-Regressionsmodell. Wir hatten:

\[ y_t \sim Po(\lambda_t e_t) \]

Bei der Poisson-Verteilung ist der Erwartungswert gleich der Varianz. Das ist bei realen Daten oft nicht realistisch. Oft ist die Varianz der Daten größer als erwartet. Um das zu berücksichtigen, kann man einen sogenannten Überdispersionsparameter zum Modell hinzu, der zusätzliche Varianz berücksichtigt:

\[ \log(\lambda_t) = \alpha+\beta x_t + \epsilon_t \]

Vergleichen Sie dazu auch das Modell der linearen Regression, wo \(\epsilon_t\) der Fehlerterm ist.

Wie bei der linearen Regression geben wir eine Normalverteilungspriori vor:

\[ \epsilon_t \sim N(0,\sigma^2) \]

\(\sigma^2=0.001\) geben wir uns hier vor.

Regularisierungspriori

• Warum geben wir uns \(\sigma^2=0.001\), also eine relativ kleine Varianz vor?

• Wir gehen hier erstmals wieder vom Ansatz weg, möglichst wenig Information in die Priori zu stecken. Diese sehr kleine Priori-Varianz sagt, dass (a priori) die \(\epsilon_t\) sehr nahe am Erwartungswert 0 sind.

• Das ist hier beabsichtigt! Die Überdispersionsparameter sollen möglichst 0 sein. Sie sollen nur zusätzliche Varianz erklären, die über unser erstes Modell hinausgeht. Deswegen: Wir “drücken \(\epsilon_t\) Richtung 0”.

• Man spricht dabei von einer Regularisierungspriori.

• Die Regularisierung ist hier auch deswegen notwendig, weil wir mehr Parameter (T+2) als Daten (T) haben. Würden wir hier mit nicht-informativen Prioris arbeiten, wäre einfach

\[ \epsilon_t = \log(\lambda_t), \]

wir würden die Daten also einfach deterministisch festlegen.

Unbekannte Parameter sind also \(\alpha, \beta\) und \(\epsilon_t\). Alternativ können wir \(\epsilon_t\) aber auch aus \(\eta_t=\log(\lambda_t), \alpha\) und \(\beta\) berechnen. Wir betrachten daher die Posteriori-Verteilung von \(\alpha, \beta, (\eta_1,\ldots,\eta_T)|y\).

Dafür müssen wir die Priori-Dichte von \(\epsilon_t\) umformen:

\[ p(\epsilon_t)\propto\exp\left(-\frac{1}{2\sigma^2}\sum_{t=1}^T\epsilon_t^2\right) = \exp\left(-\frac{1}{2\sigma^2}\sum_{t=1}^T\left(\eta_t-\alpha-\beta x_t\right)^2\right) \]

Posteriori

Wir setzten nun in der Datendichte \(\lambda_t=\exp(\eta_t)\) und erhalten die Posteriori-Dichte

\[ \begin{align*} p(\alpha,\beta,\eta|y) \propto & \, \exp\left(\sum_{t=1}^T\eta_t y_t\right)\exp\left(-\sum_{t=1}^T\exp(\eta_t)\right)\\ &\cdot\exp\left(-\frac{1}{2v_\alpha}(\alpha-m_\alpha)^2\right) \cdot\exp\left(-\frac{1}{2v_\beta}(\beta-m_\beta)^2\right)\\ &\cdot\exp\left(-\frac{1}{2\sigma^2}\sum_{t=1}^T\left(\eta_t-\alpha-\beta x_t\right)^2\right) \end{align*} \]

Für den Metropolis-within-Gibbs-Algorithmus müssen wir wieder die full conditionals berechnen.

Full conditional

In der Posteriori sind jetzt nur noch zwei Terme, die von \(\alpha\) abhängen: ein Term aus der Priori von \(\alpha\), ein Term aus der Priori von \(\epsilon\):

\[ p(\alpha|\beta,\lambda;y) \propto \exp\left(-\frac{1}{2v_\alpha}(\alpha-m_\alpha)^2\right) \cdot\exp\left(-\frac{1}{2\sigma^2}\sum_{t=1}^T\left(\eta_t-\alpha-\beta x_t\right)^2\right) \]

Dabei haben wir ist in der Exponentialfunktion die Summe von zwei quadratischen Formen. Wir landen also bei der Normalverteilung:

• \(\alpha\sim N\left(\tilde{v_\alpha}\tilde{m_\alpha},\tilde{v_\alpha}\right)\)
• \(\tilde{v_\alpha}=(1/v_\alpha+T/\sigma^2)^{-1}\)
• \(\tilde{m_\alpha}=\left(m_\alpha/v_\alpha+\sum_{t=1}^T(\eta_t-\beta x_t)/\sigma^2\right)\)

Analog erhalten wir auch für \(\beta|\alpha,\lambda;y\) eine Normalverteilung:

• \(\beta\sim N\left(\tilde{v_\beta}\tilde{m_\beta},\tilde{v_\beta}\right)\)
• \(\tilde{v_\beta}=\left(1/v_\beta+\sum_{t=1}^Tx_t^2/\sigma^2\right)^{-1}\)
• \(\tilde{m_\beta}=\left(m_\beta/v_\beta+\sum_{t=1}^T(x_t\eta_t-\alpha x_t)/\sigma^2\right)\)

Full conditional

Jetzt brauchen wir noch die Full conditional der \(\eta_t\):

\[ \begin{align*} p(\eta_t|\alpha,\beta,\eta_{-t};y) \propto \,& \exp(\eta_t y_t)\exp\left(-\exp(\eta_t)\right)\\ &\cdot\exp\left(-\frac{1}{2\sigma^2}\sum_{t=1}^T\left(\eta_t-\alpha-\beta x_t\right)^2\right) \end{align*} \]

Das ist nicht der Kern einer uns bekannten Verteilung. Hier können wir mit einem Metropolis-Schritt ansetzen.

Eine Idee, eine sinnvolle Vorschlagsverteilung zu konstruieren ist, die full conditional möglichst gut zu approximieren. Ist \(q\left(\theta^{*}|\theta^{(k-1)}\right)\approx f\left(\theta^{*}\right)\), dann wird die Akzeptanzwahrscheinlichkeit fast 1, wir sind also “fast” beim Gibbs-Sampler.

Die Implementierung des Metropolis-Hastings-Samplers in R finden Sie hier

Wir verwerfen wieder 5000 Iterationen als burn-in.

alpha.sample<-alpha.sample[-(1:5000)]
beta.sample<-beta.sample[-(1:5000)]
epsilon.sample<-epsilon.sample[-(1:5000),]

Wir haben nun 24 Parameter! \(\alpha\), \(\beta\) und \(T=22\) \(\epsilon_t\)-Parameter. Allmählich wird es schwierig, sich alle trace plots anzuschauen.

Im Folgenden können Sie sich die Trace Plot der einzelnen \(\epsilon_t\) anzeigen:

Die Schätzer selbst sind ähnlich zu den Schätzern aus dem Modell ohne Überdispersion:

	Mittelwert	Median	95%-Quantil
\(\alpha\)	9.8253	9.8251	[9.77, 9.8809]
\(\beta\)	0.0807	0.0808	[0.0526, 0.1084]

Den Überdispersionsparameter \(\epsilon\) können wir uns graphisch anschauen. Die folgende Graphik zeigt jeweils Mittelwert und (als Linien) 95%-Kredibilitätsintervalle.

Es zeigt sich ein Trend in \(\epsilon\). Hier könnte entweder eine weitere, nicht berücksichtige Kovariable versteckt sein, oder unsere Kovariable ist selbst nicht ganz passend (z.B. ist der Zeitverzug von 20 Jahren eine fragliche Annahme).