---
title: "MCMC - Tuning und Kovergenzdiagnostik"
author: "Volker Schmid"
date: "29. Mai 2017"
output:
  beamer_presentation: default
  ioslides_presentation: default

---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```

## Beispiel MCMC bei Poisson-Lognormal
```{r, echo=TRUE}
require(coda)
data <- read.table("kaiserschnitt.raw.txt", header=TRUE)

n <-dim(data)[1]
sumx <- sum(data$n)
source("poisson-lognormal.R")
mcmc.simple<-poisson.lognormal.mcmc(sumx,n)
```

## Konvergenz
MCMC-Algorithmen erzeugen eine Markovkette aus Ziehungen aus der Posteriori-Verteilung. Probleme:

* Der Algorithmus muss konvergieren, damit er aus der stationären Verteilung zieht
* Ziehungen sind automatisch abhängig. Damit ist der Algoritmus ineffizienter als unabhängiges Ziehen
* Die Effizienz hängt stark vom genauen Algorithmus ab, z.B. Metropolis-Hastings oder Gibbs-Sampler, Wahl der Vorschlagsdichte, etc..

## MCMC mit coda {.allowframebreaks} 
```{r, echo=TRUE}
summary(mcmc.simple)
```

## MCMC mit coda
```{r, echo=TRUE}
plot(mcmc.simple)
```

## Längere Ketten
```{r, cache=TRUE, echo=TRUE}
mcmc.laenger<-poisson.lognormal.mcmc(sumx,n, I=10000)
plot(mcmc.laenger)
```

## Tuning des Random Walk Proposals 
Random Walk Proposal:
$$\theta^* = \theta^{(i-1)} + u$$
oft mit $u\sim N(0,c)$. Wie wählt man $c$?

Akzeptanzrate = Anteil akzeptierter Vorschläge
* Niedrige Akzeptanzrate: Kette bleibt oft im selben Zustand $\to$ schlecht
* Zu hohe Akzeptanzrate: Kette bewegt sich (eventuell) nur langsam $\to$ schlecht

Tuning: Finde optimalen Wert für $c$. Für Random Walk Proposal gelten Akzeptanzraten zwischen ca. 30% und 60% als optimal (?)

## Beispiel zu hohe Akzeptranzrate
```{r, cache=TRUE}
plot(as.vector(mcmc.simple[,1]),type="l",ylab="mu")
```

## Poisson-Lognormal mit Tuning
```{r, cache=TRUE, echo=TRUE}
mcmc.tuning<-poisson.lognormal.mcmc(sumx,n,do.tuning=TRUE)
```

## Poisson-Lognormal mit Tuning
```{r, cache=TRUE}
plot(as.vector(mcmc.tuning[,1]),type="l",ylab="mu")
```

## Running mean plots
```{r, cache=TRUE, echo=TRUE}
par(mfrow=c(1,2))
I=1000
plot(cumsum(mcmc.simple[,1])/1:I,type="l")
plot(cumsum(mcmc.simple[,2])/1:I,type="l")
```

## Running mean plots
```{r, cache=TRUE, echo=TRUE}
par(mfrow=c(1,2))
I=1000
plot(cumsum(mcmc.tuning[,1])/1:I,type="l",)
plot(cumsum(mcmc.tuning[,2])/1:I,type="l")
```

## Auto correlation function
```{r, cache=TRUE, echo=TRUE}
par(mfrow=c(1,2))
acf(mcmc.simple[,1])
acf(mcmc.simple[,2])
```

## Auto correlation function
```{r, cache=TRUE, echo=TRUE}
par(mfrow=c(1,2))
acf(mcmc.tuning[,1])
acf(mcmc.tuning[,2])
```

# Konvergenzdiagnostik
## Gelman-Rubin-Diagnostik

* Idee: vergleiche die Varianz von mehreren parallel gelaufenen Ketten
* Bei Konvergenz sollte die Varianz in den Ketten der Varianz zwischen den Ketten entsprechen
* Within Chain Variance (unterschätzt die wahre Varianz, wenn Ketten noch nicht konvergiert)
$$ W=\frac{1}{m}\sum_{j=1}^ms_j^2, s_j^2=\frac{1}{n-1}\sum_{i=1}^n(\theta_{ij}-\bar{\theta_j})^2 $$
* Between Chain Variance 
$$ B=\frac{n}{m-1}\sum_{j=1}^m(\bar{\theta}_j-\bar{\bar{\theta}})^2; \bar{\bar{\theta}}=\frac{1}{m}\sum_{j=1}^m\bar{\theta} = j $$

## Gelman-Rubin-Diagnostik

* Geschätzte Varianz
$$ \hat{Var}(\theta)=(1-\frac{1}{n}W+\frac{1}{n}B) $$
* Potential scale reduction factor
$$ \hat{R}=\sqrt\frac{\bar{Var}(\theta)}{W} $$

* Ist $\hat{R}$ zu groß ($>1.1$?), sollten die Ketten länger laufen.

## Gelman-Rubin-Diagnostik
```{r, echo=TRUE, cache=TRUE}
require(coda)
mc.list<-parallel::mclapply(rep(sumx, 5),
            poisson.lognormal.mcmc,n=n)
print(gelman.diag(mc.list))
```

## Gelman-Rubin-Diagnostik
```{r, echo=TRUE, cache=TRUE}
require(coda)
mc.list<-parallel::mclapply(rep(sumx, 5),
            poisson.lognormal.mcmc,n=n,I=10000)
print(gelman.diag(mc.list))
```

## Gelman-Rubin-Diagnostik

* Gelman-Rubins R muss für jeden Parameter geschätzt werden
* Ziehungen aller Ketten werden dann zusammengeworfen
* Alternativ auch Aufteilung einer Kette möglich

## Geweke-Diagnostik

* Idee: Teste, ob zwei Teile einer Kette aus der selben Verteilung stammen (Teste Differenz der Mittelwerte)
* In der Regel erste 10% und letzte 50% der Kette
```{r, echo=TRUE, cache=TRUE}
print(geweke.diag(mcmc.simple))
```

## Geweke-Diagnostik
```{r, echo=TRUE, cache=TRUE}
print(geweke.diag(mcmc.tuning))
```

## Raftery und Lewis Diagnostik

* Wir interessieren uns für ein Quantil $q$. 
* Wieviele Iterationen $N$ und welchen burn-in $M$ brauchen wir, um mit Wahrscheinlichkeit $s$ innerhalb der Toleranz $r$ das Quantil $q$ zu schätzen?
* Nach einer Pilotkette lassen wir die längere Kette entsprechend laufen.

## Raftery und Lewis Diagnostik
```{r, echo=TRUE, cache=TRUE}
print(raftery.diag(mcmc.simple, 
                   q = 0.5, r = 0.05, s = 0.9))
```

## Heidelberg und Welch Diagnostik

* Teste, ob die Kette aus einer stationären Verteilung kommt
* Zuerst: Cramer-von Mises-Test mit Niveau $\alpha$ auf ganzer Kette
* Falls Nullhypothese abgelehnt, verwirf erste 10%, 20%, ..., bis zu 50% 
* Falls Nullhypothese angenommen: Half-width-test 
* Berechne $(1-\alpha)$-Kredibilitätsintervall (KI)
* Teststatistik: Hälfte der Breite des KI durch Mittelwert

## Heidelberg und Welch Diagnostik

```{r, echo=TRUE, cache=TRUE}
print(heidel.diag(mcmc.simple))
```

## Heidelberg und Welch Diagnostik

```{r, echo=TRUE, cache=TRUE}
print(heidel.diag(mcmc.tuning))
```