Statistische Informationsverarbeitung - Übung 3

From StatWiki
Jump to: navigation, search

Angabe: Media:StatInfMCMCue3.pdf

Accept/Reject

Verwenden sie den Accept/Reject Algorithmus, um  X \sim N(0;1) \, mit  f(x) \propto e^{- \frac 1 2 x^2}\, zu simulieren.

Benutzen Sie dazu eine zweiseitige Exponentialverteilung als Kandidatsverteilung, also  g(y) \propto e^{-|y|}\,.

(" \propto \," bedeutet "proportional zu",  g(y) \propto \tilde g(y) \, heißt also, dass  g(y) = c \tilde g(y) \, für eine Konstante c \, ist.)

(a)

Um  Y \sim g(y)\, zu erzeugen, benutzen Sie die Inverse Transformationsmethode für eine "normale" Exponentialverteilung und setzen  y := -y\, mit Wahrscheinlichkeit \frac 1 2 \,.

1. generiere y \, nach Y \sim g(y) \,

 F_{\mathrm{exp}}(x) = 1 - e^{-x} \, 
 
 F_{\mathrm{exp}}^{-1}(x) \propto \log(1 - x) \, 

Eine "normale" Exponentialverteilung 

> pdf(rpdf)
>
>  G1 = function(y) { -log(1 - y) }
>  y = G1(runif(10000,0,1))
>
>  hist(y,main=expression(tilde(g)(y)), freq=FALSE, xlim=c(-4,4),ylim=c(0,1))
>  yv=seq(-4,4,length.out=1000)
>  lines(yv, dexp(yv),col=2)
Jene mit "... und setzen  y := -y\, mit Wahrscheinlichkeit \frac 1 2 \,."
REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/f134f4e020963399fe74805081a26cab9453e971_%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)
>  y = y * c(-1,1)[rbinom(10000,1,.5) + 1]
Error: object 'y' not found
Execution halted
in
pdf(rpdf)
y = y * c(-1,1)[rbinom(10000,1,.5) + 1]
hist(y,main=expression(tilde(g)(y)), freq=FALSE, xlim=c(-4,4),ylim=c(0,.45))

(b)

Generieren Sie so  n = 10000\, (oder mehr) normalverteilte Zufallszahlen, erstellen Sie ein Histogramm und plotten Sie zur Überprüfung die Dichtefunktion der Normalverteilung in Ihr Histogramm. Schreiben Sie Ihr Programm so, dass es genau die gewünschte Anzahl an Zufallszahlen liefern kann.

2. generiere u \sim \mathrm{unif}(0;1) \,
3. Falls  u \le \frac{ \tilde{f}(y) }{ \tilde M \tilde{g}(y) } \,, dann setze x := y\, sonst gehe zu 1.

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/3aa47e85e57d4f438a3d5366d95c73b8be89f5c1_%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)
>
>  g=function(x) { .5 * exp(-.5 * x) }
>  f=function(x) { 1/(2 * pi) * exp(-.5 * x^2) }
>
>  ft = function(x) { exp(-.5 * x^2) }
>  gt = function(y) { exp( - abs(y) ) }
>
>  M = ft(.5)/gt(.5)
>
>  x = 1:10000
>  yNeeded = rep(0,length(x))
>
>   for(i in x)
+   {
+     # no do { } while(); construct in R :(
+
+     repeat
+     {
+       yNeeded[i] = yNeeded[i] + 1
+       y = G1(runif(1,0,1)) * c(-1,1)[rbinom(1,1,.5) + 1]
+       u=runif(1,0,M)
+       if(ft(y) / (M * gt(y)) >= u)
+         break;
+     }
+     x[i] = y
+   }
Error: could not find function "G1"
Execution halted
in
pdf(rpdf)

g=function(x) { .5 * exp(-.5 * x) } 
f=function(x) { 1/(2 * pi) * exp(-.5 * x^2) }
ft = function(x) { exp(-.5 * x^2) }
gt = function(y) { exp( - abs(y) ) }
M = ft(.5)/gt(.5)
x = 1:10000
yNeeded = rep(0,length(x))
 for(i in x)
 {
   # no do { } while(); construct in R :(
   
   repeat
   {
     yNeeded[i] = yNeeded[i] + 1
     y = G1(runif(1,0,1)) * c(-1,1)[rbinom(1,1,.5) + 1]
     u=runif(1,0,M)    
     if(ft(y) / (M * gt(y)) >= u)
       break;
   }
   x[i] = y
 }
hist(x, main="f(x)", freq=FALSE, xlim=c(-4,4),ylim=c(0,.45))
xv = seq(-4,4,length.out=1000)
lines(xv, dnorm(xv), col=2)

(c)

Wieviele Iterationsschritte mussten Sie durchschnittlich vornehmen, um zu einem "accept-Schritt" zu gelangen? Wieviele Schritte haben Sie erwartet? Warum?

ACHTUNG das stimmt so nicht. Um m-tilde zu finden muss man e^{|y|}\, in der Ableitung verwenden. und das maximum wird bei x=1 angenommen.
REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/f5f35d6ddf9bf0c3310002764866f179a60d1a0c_%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 --->
> mean(yNeeded)
Error in mean(yNeeded) : object 'yNeeded' not found
Execution halted
in
mean(yNeeded)
Also erwartet habe ich einen Schritt, weil ich M = \sup_x \frac{f(x)}{g(x)} \, gesetzt habe.

X \sim N(0;1) \, mit f(x) = \frac{1}{2 \pi} e^{- \frac 1 2 x^2} \,

Y \sim \mathrm{exp}(1 / 2)\, mit g(y) = \lambda e^{- \lambda y}, \lambda = \frac 1 2 \,

Herleitung des Supremums ("Jener Funkionswert (nicht das x) das immer größer/gleich der Funktion ist").

\frac{\partial}{\partial x} \frac{f(x)}{g(x)} = 
\frac{\partial}{\partial x} \frac{\frac 1 \pi e^{-\frac 1 2 x^2}}{\frac 1 2 e^{-\frac 1 2 x}} = 
\frac{\partial}{\partial x} \frac 1 \pi e^{\frac{x - x^2}{2}} =
\frac 1 \pi e^{\frac{x - x^2}{2}} (1 / 2 - x)
\,

\frac 1 \pi e^{\frac{x - x^2}{2}} (1 / 2 - x) = 0 \Leftrightarrow x = \frac 1 2 \,

\frac{f(\frac 1 2)}{g(\frac 1 2)} \approx 0.3606924\,0.
> pdf(rpdf)
>
> x=seq(-5,5,length.out=1000)
> g=function(x) { .5 * exp(-.5 * x) }
> f=function(x) { 1/(2 * pi) * exp(-.5 * x^2) }
> plot(x,f(x)/g(x),type="l")
> abline(v=.5,col=2)


Ich glaub in der VO ist der Buchstabe M doppelt verwendet worden, oder ich blick grad nicht durch.

Erwartet habe ich M\, Schritte, weil T\, (die Anzahl der Iterationsschritte bis zum akzeptieren eines Y\,) geometrisch verteilt ist. 

Wie unten ersichtlich P(\mathrm{reject}) = 1 - \frac 1 M\,, daher ist P(\mathrm{accept}) = 1 - P(\mathrm{reject}) = \frac 1 M \,.
 
T \sim \mathrm{geom}(p = \frac 1 M) \, mit \operatorname{E}(T) = \frac 1 p = M \,

(c) Beweis zu ACCEPT/REJECT

Beweis das ACCEPT/REJECT ein X \sim f(x) \, von der gewünschten Verteilung mit Dichtefunktion f(x)\, generiert.

  P(X \le x) = \int_{-\infty}^x f(z) dz \,    

 
P(\mathrm{reject}) = 
\int P(\mathrm{reject} \, Y | Y=y) g(y) \, dy = 
\, nach Satz d. totalen Wahrscheinlichkeit  P(Z < a) = \sum_z (P < a|Z=z) P(Z=z) \,
 = \int P(U > \frac{f(y)}{M g(y)} | Y = y) g(y) \, dy = \,
 = \int \left(1 -  \frac{f(y)}{M g(y)} \right) g(y) \, dy = \, 
 = \int g(y) \, dy - \frac 1 M \int f(y) \, dy = 1 - \frac 1 M \, da eine Dichte auf 1 summieren muss.

 P(\mathrm{accept} \, Y, Y \le x) = 
 = \int P(\mathrm{accept} \, Y, Y \le x | Y=y) g(y) \, dy = \,
 = P(U \le \frac{f(y)}{M g(y)}, y \le x | Y=y) g(y) \, dy = \,
 = \int_{-\infty}^x \frac{f(y)}{M g(y)} g(y) \, dy = \,
 = \frac 1 M \int_{-\infty}^x f(y) \, dy 

P(X \le x) = \sum_{n=1}^\infty P(\mathrm{accept} \, Y \mathrm{beim} \, n \mathrm{-ten} \, \mathrm{Versuch}, Y \le x) = \,
 = \sum_{n=1}^\infty P(\mathrm{reject})^{n-1} \cdot \frac 1 M \int_{-\infty}^x f(y) \, dy =  (Die REJECT Wahrscheinlichkeiten dürfen miteinander multipliziert werden, da die Versuche unabhängig voneinander sind)
 = \sum_{n=1}^\infty (1 - \frac 1 M)^{n-1} \cdot \frac 1 M \int_{-\infty}^x f(y) \, dy =  (Geometrische Summenformel)
 = \int_{-\infty}^x f(y) \, dy
 \sum_{n=1}^\infty (1 - \frac 1 M)^{n-1} = \frac{1}{1 - (1 - \frac 1 M)} = M \,

(d)

Schreiben Sie nun Ihr Programm um, sodass Sie anstelle einer vorgegebenen Anzahl von  X_i^{'} s\, eine eine vorgegebene Zahl an Kandidaten (also  Y_i^{'} s\,) generieren.

Wenn Sie nun in etwa n \, (zB  n = 10000\,) normalverteilte  X_i^{'} s\, möchten, wieviele Kandidaten müssen Sie erzeugen?

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/b10ac860f4168093baf8b5b02854c7803e5d41d8_%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)
>
>   n = 10000
>
>   y = G1(runif(n,0,1)) * c(-1,1)[rbinom(n,1,.5) + 1]
Error: could not find function "G1"
Execution halted
in
 pdf(rpdf)
 n = 10000  
 
 y = G1(runif(n,0,1)) * c(-1,1)[rbinom(n,1,.5) + 1]
 u=runif(n,0,M)    
 x = y[ft(y) / (M * gt(y)) >= u]
 hist(x, main="f(x)", freq=FALSE, xlim=c(-4,4),ylim=c(0,.45))
 xv = seq(-4,4,length.out=1000)
 lines(xv, dnorm(xv), col=2)