Statistische Informationsverarbeitung - Übung 11

From StatWiki
Jump to: navigation, search

Gibbs Sampler Programmieren Sie einen Gibbs Sampler für die folgende Verteilung

\pi(x,y) \propto y \exp \{ - (x + 2xy + y) \} \quad x,y \ge 0\,

(a) Bestimmen Sie die bedingten Verteilungen X|Y\, und Y|X\,.

Partielle Integration:  \int f \cdot g = F \cdot g - \int F \cdot g' \,

 \begin{align}
\pi_X
& = \int_0^\infty \pi_{X,Y}(x,y) \, dy \\
& \propto \int_0^\infty y \cdot e^{ - (x + 2xy + y) } \, dy \\
& = \int_0^\infty y \cdot e^{ - y (2x + 1) } e^{-x} \, dy \\
& = e^{-x} \left[ y \cdot \frac{e^{ - y (2x + 1)}}{-(2x+1)} \Bigg|_{y=0}^\infty - \int_0^\infty e^{ - y (2x + 1) } \, dy \right] \\
& = e^{-x} \left[ \frac{1}{-(2x+1)} - \frac{e^{ - y (2x + 1)}}{-(2x+1)} \Bigg|_{y=0}^\infty \right] \\
& = e^{-x} \frac{1}{(2x+1)^2} \\
\end{align} 

 \begin{align}
\pi_Y
& = \int_0^\infty \pi_{X,Y}(x,y) \, dx \\
& \propto \int_0^\infty y \cdot e^{ - (x + 2xy + y) } \, dx \\
& = \int_0^\infty y \cdot e^{ - x (2y + 1) } e^{-y} \, dx \\
& = y \cdot  e^{-y} \frac{e^{ - x (2y + 1) }}{- (2y + 1)} \Bigg|_{x=0}^\infty \\
& = e^{-y} \frac{y}{2y + 1}
\end{align} 
 \begin{align}
\pi_{Y|X}(y|x)
& = \frac{\pi_{X,Y}(x,y)}{\pi_X(x)} \\
& \propto \frac{y \cdot e^{-(2xy + x + y)}}{ \frac{e^{-x}}{(2x+1)^2}} \\
& = y (2x + 1)^2 e^{-2xy - y - x + x} \\ 
& = y (2x + 1)^2 e^{-y (2x + 1)} \quad \alpha = 2, \beta=(2x+1) \\
& \propto \Gamma^1(\alpha) \beta^\alpha y^{\alpha - 1} e^{-y \beta} \\
& = \mbox{Gamma}(2, 2x + 1)
\end{align}

 \begin{align}
\pi_{X|Y}(x|y)
& = \frac{\pi_{X,Y}(x,y)}{\pi_Y(y)} \\
& \propto \frac{y \cdot e^{-(2xy + x + y)}}{ e^{-y} \frac{y}{- (2y + 1)} } \\
& = \frac{y \cdot e^{-2xy - x - y + y)} (-2y - 1)}{ y } \\
& = - (2y + 1) e^{-x (2y + 1)} \qquad \lambda = 2y + 1 \\
& = - \lambda e^{-\lambda x} \\
& = \mbox{Exponential}(\lambda)
\end{align} 
Kommentar Schneider: Die bedingten Verteilungen könnte man auch "direkt" ablesen.

(b) Programmieren Sie den Gibbs Sampler und wählen Sie eine geeignete Anzahl von Iterationen mit der Gelman-Rubin Statistik.

> pdf(rpdf,width=12)
> set.seed(12354)
> library(coda)
Loading required package: lattice
>
> gibbs = function(y0,n=5e3) {
+         d = replicate(n,
+                 {
+                         x0 <<- rexp(1,2 * y0 + 1)   # X|Y
+                         y0 <<- rgamma(1,2,2*x0 + 1) # Y|X
+                         c(x0,y0)
+                 })
+         d = t(d)
+         colnames(d) = c("x","y")
+         mcmc(d)
+ }
>
> # gelman.plot(mcmc.list(sapply(rep(1,4), function(y0) gibbs(y0,1e3), simplify=F)))
>
> mc = mcmc.list(sapply(seq(1,50,length.out=10), gibbs, simplify=F))
>
> gelman.plot(mc)
Ab ca. 2000 Iterationen wird die Gelman-Rubin Statistik schon sehr klein und daher sind 5000 Iteration ausreichend.

(c) Berechnen Sie analytisch (bis auf eine Konstante) die Dichte f_Y(y)\, der Randverteilung für Y\,.

Siehe Beispiel (a)

(d) Plotten Sie ein Histogramm für Y\, mit den simulierten Werten und zeichnen Sie im Histogramm die vorher berechnete Dichtefunktion ein (dazu können Sie verwenden das \int_0^\infty \frac{y}{1+2y} e^{-y} \, dy \approx 0.2692723\, ist)

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/6e5a1647fbdceac301067bff19679fa437469d2a_%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)
> library(coda)
Loading required package: lattice
>
> fy = function(y) exp(-y)*(y/(1+2*y))/0.2692723
>
> hist(mc[[1]][,2],freq=F,main="Histogram and Density",xlab=expression(Y[i]),breaks=50)
Error in hist(mc[[1]][, 2], freq = F, main = "Histogram and Density",  :
  object 'mc' not found
Execution halted
in

pdf(rpdf) library(coda)

fy = function(y) exp(-y)*(y/(1+2*y))/0.2692723

hist(mc[[1]][,2],freq=F,main="Histogram and Density",xlab=expression(Y[i]),breaks=50) curve(fy,0,8,col=2,add=T)