Statistische Informationsverarbeitung - Übung 12

From StatWiki

Jump to: navigation, search

Contents


Beispiel 18

Beispiel aus der VO zu Hybrid MCMC

\pi_{X_1,X_2}(x_1,x_2) = \frac 3 2 e^3 \exp\{-x_1 - |x_2| e^{2x_1} \}, \quad x_1 >1, x_2 \in \mathbb{R} \,

Programmieren Sie den sampler und plotten Sie Ihre simulierten Punkte in der Ebene mit Konturlinien der gemeinsamen Dichtefunktion. 

Als Kandidatsverteilung wird 

q_1(y) =1_{(1,\infty)} (y) e^{-y} \,

gewählt.

Die passende Konstante für eine Dichte ist
 
 \begin{align}
c \cdot \int_1^\infty e^{-y} 
& = c \cdot  -e^{-y} \Big|_1^\infty \\
& = c \cdot e^{-1}  \\
& = 1 \Rightarrow c = e
\end{align}


Durch Integration kommt man auf die Verteilungsfunktion

\begin{align}
F(x) 
& = \int_{y=1}^x e^{1-y} \, dy \\
& = -e^{1-x} \Big|_1^x \\
& = -e^{1-x} + e^0 \\
& = 1 - e^{1-x}
\end{align}

Das ist aber die Verteilungsfunktion einer Expontentialverteilten Zufallsvariable und daher kann y mit 'rexp' simuliert werden.

pdf(rpdf)
f = function(x1,x2) 3/2 * exp(3) * exp(-x1-abs(x2)*exp(2*x1))

x1 = seq(1.00001, 2.7, length.out=50)
x2 = seq(-.4, .4, length.out=50)

# library(distr)
library(coda)

f = function(x1,x2) 3/2 * exp(3 - x1 - abs(x2) * exp(2*x1))

gibbs = function(x1,n=4000) {
 
	# sehr langsam
	# D = DExp(rate = 1)

	d = matrix(0, nrow=n, ncol=2)
	colnames(d) = c("x1","x2")

	qx1 = 1 + rexp(n,1)

	ru = runif(n)
	rb = (rbinom(n,1,.5)*2)-1

	for(i in 1:n) {
		# sehr langsam
		# X2|X1
		# rate(D) = exp(2 * x1)
		# d[i,2] = x2 = r(D)(1)

		d[i,2] = x2 = rb[i] * rexp(1, exp(2*x1))

		y = qx1[i]
		d[i,1] = x1 = ifelse(ru[i] < f(y, x2) / f(x1, x2), y, x1)
	}

        mcmc(d)
}

# mc = mcmc.list(sapply(rep(1,2), gibbs, simplify=F))
mc = mcmc.list(gibbs(1,10000))

plot(c(1,2.5),c(-.4,.4),type="n",xlab="x1",ylab="x2")
points(unlist(mc[,1]), unlist(mc[,2]), pch=".", cex=1.5)
contour(x1, x2, z, add=T)

#

Hier die Dichtefunktion noch in 3D. 

pdf(rpdf)
persp(x1,x2,z,theta=60,col = "lightblue")

#

Beispiel 19

Zeigen Sie, dass der Gibbs Sampler ein Spezialfall der allgemeinen Version des MH- Algorithmus ist. 

\vec X = (X_1, X_2, \ldots, X_n) \, ist Zufallsvektor.

\vec{X}_{-1}^{(t)} = \left(X_1^{(t+1)}, \ldots, X_{i-1}^{(t+1)}, X_{i+1}^{(t)}, \ldots, X_{p}^{(t)} \right)\, ist Zufallsvektor.

Bezeichne \pi_i(x_i| \vec{x}_{-i}^{(t)})\, die bedingte Dichten.

Als Vorschlagsdichte: q_i(y|x_i, \vec{x}_{-1}^{(t)}) = \pi_i(y_i| \vec{x}_{-i}^{(t)}) \,

Die Annahmewahrscheinlichkeit ist 

\begin{align}
\alpha_i( (x_i^{(t)}, \vec{x}_{-1}^{(t)}), (y_i, \vec{x}_{-1}^{(t)}))
& = \min\left( 1, \frac{ \pi_i(y_i| \vec{x}_{-i}^{(t)}) }{ \pi_i(x_i| \vec{x}_{-i}^{(t)}) } \frac{ q_i(x_i|y_i,\vec{x}_{-i}^{(t)})} { q_i(y_i|x_i,\vec{x}_{-i}^{(t)}) } \right) \\
& = \min\left( 1, \frac{ \pi_i(y_i| \vec{x}_{-i}^{(t)}) }{ \pi_i(x_i| \vec{x}_{-i}^{(t)}) } \frac{ \pi_i(x_i| \vec{x}_{-i}^{(t)}) } { \pi_i(y_i| \vec{x}_{-i}^{(t)}) } \right) \\
& = 1
\end{align}
Retrieved from "http://www.eisber.net/StatWiki/index.php/Statistische_Informationsverarbeitung_-_%C3%9Cbung_12"