Statistische Informationsverarbeitung - Übung 3

From StatWiki

Jump to: navigation, search

Angabe: Media:StatInfMCMCue3.pdf

Contents

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)
 
#
> sys.load.image("/var/www/localhost/htdocs/StatWiki/Rfiles/true", TRUE)
Error in loadNamespace(name) : there is no package called 'distrEx'
Calls: sys.load.image ... tryCatch -> tryCatchList -> tryCatchOne -> 
Execution halted



Jene mit "... und setzen  y := -y\, mit Wahrscheinlichkeit \frac 1 2 \,."

 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))
 
#
> sys.load.image("/var/www/localhost/htdocs/StatWiki/Rfiles/true", TRUE)
Error in loadNamespace(name) : there is no package called 'distrEx'
Calls: sys.load.image ... tryCatch -> tryCatchList -> tryCatchOne -> 
Execution halted

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


 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)
 
#
> sys.load.image("/var/www/localhost/htdocs/StatWiki/Rfiles/true", TRUE)
Error in loadNamespace(name) : there is no package called 'distrEx'
Calls: sys.load.image ... tryCatch -> tryCatchList -> tryCatchOne -> 
Execution halted

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

 mean(yNeeded)
 
#


Error in loadNamespace(name) : there is no package called 'distrEx'
Calls: sys.load.image ... tryCatch -> tryCatchList -> tryCatchOne -> 
Execution halted



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? 

  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)

 
#
> sys.load.image("/var/www/localhost/htdocs/StatWiki/Rfiles/true", TRUE)
Error in loadNamespace(name) : there is no package called 'distrEx'
Calls: sys.load.image ... tryCatch -> tryCatchList -> tryCatchOne -> 
Execution halted
Retrieved from "http://www.eisber.net/StatWiki/index.php/Statistische_Informationsverarbeitung_-_%C3%9Cbung_3"
Personal tools