In volkerschmid/bayeskurs: Bayes-Kurs

knitr::opts_chunk$set(echo = TRUE, cache = TRUE)

Beispiel MCMC bei Poisson-Lognormal

library(coda)
library(bayeskurs)
data("kaiserschnitt.raw")
data<-kaiserschnitt.raw

n <-dim(data)[1]
sumx <- sum(data$n)
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}

summary(mcmc.simple)

MCMC mit coda

plot(mcmc.simple)

Längere Ketten

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

plot(as.vector(mcmc.simple[,1]),type="l",ylab="mu")

Poisson-Lognormal mit Tuning

mcmc.tuning<-poisson.lognormal.mcmc(sumx,n,do.tuning=TRUE)

Poisson-Lognormal mit Tuning

plot(as.vector(mcmc.tuning[,1]),type="l",ylab="mu")

Running mean plots

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

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

par(mfrow=c(1,2))
acf(mcmc.simple[,1])
acf(mcmc.simple[,2])

Auto correlation function

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

require(coda)
mc.list<-parallel::mclapply(rep(sumx, 5),
            poisson.lognormal.mcmc,n=n)
print(gelman.diag(mc.list))

Gelman-Rubin-Diagnostik

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

print(geweke.diag(mcmc.simple))

Geweke-Diagnostik

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

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

print(heidel.diag(mcmc.simple))

Heidelberg und Welch Diagnostik

print(heidel.diag(mcmc.tuning))

volkerschmid/bayeskurs documentation built on Dec. 23, 2021, 4:10 p.m.