Statistische Informationsverarbeitung - Übung 8

From StatWiki
Jump to: navigation, search

Integralabschätzung mit Hilfe von MCMC

Schätzen Sie das Integral

 I = \int_{-\infty}^{\infty} \sqrt{|x|} e^{-e^{x^2}} \,dx\,

mithilfe von MCMC ab. Dazu können Sie zum Beispiel  \pi(x) \propto e^{-e^{x^2}}  \, wählen und das Integral als Erwartungswert

 I = c \cdot \operatorname{E}_{\pi} \left[ \sqrt{|X|} \right] \, mit  c = \int_{-\infty}^{\infty} e^{-e^{x^2}} \,dx \,

auffassen.

Kurze Erklärung zum Zusammenhang von Integral und Erwartungswert

 \pi(x) \propto e^{-e^{x^2}} = f(x) \, (Nur proportional weil es ja eine Dichtefunktion sein muss)

 \pi(x) = d \cdot f(x) \,

 d = \frac{1}{\int_{-\infty}^{\infty} e^{-e^{x^2}} \,dx} \,

 c = \int_{-\infty}^{\infty} e^{-e^{x^2}} \,dx \,

 \begin{align}
I 
= &  \int_{-\infty}^{\infty} \sqrt{|x|} e^{-e^{x^2}} \,dx \\
= &  \int_{-\infty}^{\infty} \sqrt{|x|} f(x) \,dx \\
= &  \frac 1 d \int_{-\infty}^{\infty} \sqrt{|x|} d \cdot f(x) \,dx \\
= &  \frac 1 d \int_{-\infty}^{\infty} \sqrt{|x|} \pi(x) \,dx \\
= &  \frac 1 d \operatorname{E}_{\pi}\left[ \sqrt{|X|} \right] \\
= &  c \cdot \operatorname{E}_{\pi}\left[ \sqrt{|X|} \right] 
\end{align} \,


(a)

Approximieren Sie erst

 \operatorname{E}_{\pi} \left[ \sqrt{|X|} \right] \,

durch das (entsprechende) arithmetische Mittel einer MH-Kette mit obigem  \pi \, und  q(x,y) \propto 1_{[x-\delta,x+\delta]}(x) \,, dh  Y_t|X_t \sim \mbox{unif}[x_t - \delta, x_t + \delta ]\, (random walk Metropolis).

Entscheiden Sie mithilfe von trace plots der benötigten fortlaufenden arithmetischen Mittel under der Gelman-Rubin Statistik, wie lange sie die MH-Kette laufen lassen.

> mh = function(n,pix,q,x0) {
+         c(x0,
+         replicate(n-1, {
+                 # q(x)
+                 yt = qm(x0)
+                 # accept?
+                 ifelse(runif(1) < (pix(yt) / pix(x0)), x0 <<- yt, x0)
+                 }))
+ }
>
> pix = function(x) exp(-exp(x^2))
> qm = function(x,delta=3) x + runif(1,-delta,delta)</pre>

 Abschätzen des Wertebereichs von pi(x)\,
REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/b0950f2138cdd906bf8d2364c64e5c05fb7b14a2_%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)
> x=seq(-1.5,1.5,length.out=5000) ; plot(x, pix(x), ylab=expression(pi(x)), type="l")
Error in xy.coords(x, y, xlabel, ylabel, log) :
  could not find function "pix"
Calls: plot -> plot.default -> xy.coords
Execution halted
in

pdf(rpdf) x=seq(-1.5,1.5,length.out=5000) ; plot(x, pix(x), ylab=expression(pi(x)), type="l")

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/30cc396b2c9ab7f6bbf197752bd66853573a026e_%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,width=24,height=12)
> library(coda)
Loading required package: lattice
>
> mcmcdata = mcmc.list(sapply(seq(-1,1,length.out=8),
+ 	function(x) {
+ 		mc = as.mcmc(mh(3000, pix, qm, pix(x)))
+ 		# fortlaufendes Mittel f. jede d. m Ketten
+ 		attr(mc, "cumMeans") = cumsum(mc)/1:length(mc)
+ 		mc
+ 	} ,simplify=F))
Error in as.mcmc(mh(3000, pix, qm, pix(x))) :
  could not find function "mh"
Calls: mcmc.list -> sapply -> lapply -> FUN -> as.mcmc
Execution halted
in

pdf(rpdf,width=24,height=12) library(coda)

mcmcdata = mcmc.list(sapply(seq(-1,1,length.out=8), function(x) { mc = as.mcmc(mh(3000, pix, qm, pix(x))) # fortlaufendes Mittel f. jede d. m Ketten attr(mc, "cumMeans") = cumsum(mc)/1:length(mc) mc } ,simplify=F))

  1. fortlaufendes Mittel ueber alle Ketten

attr(mcmcdata,"cumMeans") = rowMeans(sapply(mcmcdata, attr, "cumMeans"))

plot.chain.i<-function(i,mcdata,exact) {

  mc=mcmcdata[[i]]
  traceplot(mc,smooth=F,col=gray(.8),lwd=2 ,main=paste("Enhanced Traceplot Chain #",i,sep="")) 
  legend(0, max(mc)*.99, c("Markov Chain", "exakt","fortlaufendes Mittel", "gesamt fortlaufendes Mittel"), 

lty="solid", col=c(gray(.8), 2, 3, 5))

  points(attr(mc,"cumMeans"),t='l',col=3,lwd=2)
  points(attr(mcdata, "cumMeans"),t='l',col=5,lwd=2)
  abline(h=exact,lwd=2,col=2)

}

par(mfrow=c(2,4)) for(i in 1:length(mcmcdata)) plot.chain.i(i, mcmcdata, 0) par(mfrow=c(1,1))

In den vorherigen Plots sieht man recht schön, dass ab ca. 500 Iterationen das Burn-In startet.
 
Im folgenden Plot wird die Gelman-Rubin Statistik gezeigt. Die Anzahl der Iterationen sollte ausreichen.
REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/1ce3e33802554d6cfcfb11cf266ce8696fa97840_%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
>
> gelman.plot(mcmcdata)
Error in as.mcmc.list(x) : object 'mcmcdata' not found
Calls: gelman.plot -> as.mcmc.list
Execution halted
in

pdf(rpdf) library(coda)

gelman.plot(mcmcdata)

(b)

Schätzen Sie nun

 c = \int_{-\infty}^{\infty}  e^{-e^{x^2}} \,dx \,

zum Beispiel mithilfe von  Z_1, \ldots, Z_n \sim_{iid} N(0,1)\, ab indem Sie

 c \approx \frac 1 n \sum_{i=1}^n h(Z_i) \,

mit entsprechender Funktion h berechnen.

Bestimmung von h(x)\,

\int g(x) \,dx = \int f(x) \frac{g(x)}{f(x)} \,dx = \operatorname{E}[h(x)] 
 
 f(x) = \frac{1}{\sqrt{2 \pi}} e^{- \displaystyle \frac{x^2}{2} } \, Dichte einer  N(0,1) \, Verteilung
 
 g(x) =  e^{-e^{x^2}} \,

 h(x) = \frac{g(x)}{f(x)} = \sqrt{2 \pi} e^{\frac{x^2}{2} - e^{x^2}} \,
> g = function(x) exp(-exp(x^2))
> h = function(x) sqrt(2 * pi)*exp(x^2/2-exp(x^2))
> con=mean(h(rnorm(100000)))
> con
[1] 0.5263933
Mit numerischer Integration im Vergleich dazu:

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/507fde7fef739f53a7bdd491f8d3429dd97e0d51_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'text'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> integrate(g, -100, 100)
Error in match.fun(f) : object 'g' not found
Calls: integrate -> match.fun
Execution halted
in

integrate(g, -100, 100)

(c)

Setzen Sie Ihre Ergebnisse aus (a) und (b) zusammen, um eine Annäherung an  I \, zu bekommen.

Übrigens: Wir werden uns noch damit beschäftigen, wie man Konfidenzintervalle für solche Abschätzungen bekommen kann.

> pdf(rpdf)
> f= function(x) sqrt(abs(x)) * exp(-exp(x^2))
> x=seq(-2,2,length.out=1000) ; plot(x,f(x), type="l", main="Funktion die Integriert werden soll")
Numerische Integration:
REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/cf5a65ec01c6680c6498dab7cf8b31b743298ec1_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'text'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> integrate(f, -2, 2)
Error in match.fun(f) : object 'f' not found
Calls: integrate -> match.fun
Execution halted
in

integrate(f, -2, 2)

Die Annäherung von  I = c \cdot \operatorname{E}_{\pi} \left[ \sqrt{|X|} \right] \,:
 
(Es wird erst ab dem Burn-In gestartet)

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/61f178eee4b78bc615e92261089704e77febc230_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'text'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> con * mean(sapply(mcmcdata, function(mc) sqrt(abs(mc[-c(1:500)]))))
Error: object 'con' not found
Execution halted
in

con * mean(sapply(mcmcdata, function(mc) sqrt(abs(mc[-c(1:500)]))))