Statistische Informationsverarbeitung - Übung 12

From StatWiki
Jump to: navigation, search

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.
REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/6de4d3b0f7b6165d96cb4584ca2125de7df325a0_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'display'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> 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)
Loading required package: lattice
>
> 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)
Error in contour.default(x1, x2, z, add = T) : object 'z' not found
Calls: contour -> contour.default
Execution halted
in

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)

  1. 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)

}

  1. 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.

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/206aed0324cb237f6a8e40f3723da97ea485f66b_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'display'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> pdf(rpdf)
> persp(x1,x2,z,theta=60,col = "lightblue")
Error in persp(x1, x2, z, theta = 60, col = "lightblue") :
  object 'x1' not found
Execution halted
in

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}