Statistische Informationsverarbeitung - Übung 2

From StatWiki
Jump to: navigation, search

Angabe: Media:StatInfMCMCue2.pdf

Zu Box-Muller

Vervollständigen Sie den Beweis zur Box-Muller Methode aus der Vorlesung. Sie dürfen dazu (als Konvention) annehmen, dass die Funktion \arctan(x) \, Werte in (0; 2 \pi) \, annimmt.

Partielle Ableitungen von g^{-1}(x_1,x_2)\,

 { \partial g_1^{-1}(x_1,x_2) \over \partial x_1} = { \partial e^{- \frac 1 2 (x_1^2 + x_2^2) } \over \partial x_1} = e^{- \frac 1 2 (x_1^2 + x_2^2) } (-x_1) 

 { \partial g_1^{-1}(x_1,x_2) \over \partial x_2} = { \partial e^{- \frac 1 2 (x_1^2 + x_2^2) } \over \partial x_2} = e^{- \frac 1 2 (x_1^2 + x_2^2) } (-x_2) 

 { \partial g_2^{-1}(x_1,x_2) \over \partial x_1} = { \partial {1 \over 2 \pi} \arctan \left( {x_2 \over x_1} \right) \over \partial x_1} = -{x_2 \over {2 \pi (x_1^2 +x_2^2)}}  

 { \partial g_2^{-1}(x_1,x_2) \over \partial x_2} = { \partial {1 \over 2 \pi} \arctan \left( {x_2 \over x_1} \right) \over \partial x_2}  = {x_1 \over {2 \pi (x_1^2 +x_2^2)}} 
 \det \left( { \partial g^{-1}(x_1,x_2) \over \partial (x_1, x_2) } \right) =  
 
 \begin{vmatrix} 
 -x_1 e^{- \frac 1 2 (x_1^2 + x_2^2) } & -x_2 e^{- \frac 1 2 (x_1^2 + x_2^2) } \\
 -{x_2 \over {2 \pi (x_1^2 +x_1^2)}} & {x_1 \over {2 \pi (x_1^2 +x_1^2)}}
 \end{vmatrix} = 


 
  =  \left(-x_1 e^{- \frac 1 2 (x_1^2 + x_2^2) } \right) \cdot \left({x_1 \over {2 \pi (x_1^2 +x_1^2)}} \right)
  - \left(- x_2 e^{- \frac 1 2 (x_1^2 + x_2^2) } \right) \cdot \left(-{x_2 \over {2 \pi (x_1^2 +x_1^2)}} \right) 




  = {-x_1^2 e^{- \frac 1 2 (x_1^2 + x_2^2) } \over {2 \pi (x_1^2 +x_1^2)}} - {x_2^2 e^{- \frac 1 2 (x_1^2 + x_2^2) } \over {2 \pi (x_1^2 +x_1^2)}}
  = - {(x_1^2 + x_2^2) e^{- \frac 1 2 (x_1^2 + x_2^2) } \over {2 \pi (x_1^2 +x_1^2)}} 
  = - \frac{1}{2 \pi} e^{- \frac 1 2 (x_1^2)} \cdot \frac{1}{2 \pi} e^{- \frac 1 2 (x_2^2)}
 

Anwendung des Transformationssatzes: 


  f_{x_1,x_2}(x_1,x_2) = f_{u_1,u_2}(g_1^{-1}(x_1,x_,2), g_2^{-1}(x_1,x_,2))  \left| \det \left( { \partial g^{-1}(x_1,x_2) \over \partial (x_1, x_2) } \right) \right| = 
  


  1 \cdot 1 \cdot \left| - \frac{1}{2 \pi} e^{- \frac 1 2 (x_1^2)} \cdot \frac{1}{2 \pi} e^{- \frac 1 2 (x_2^2)} \right| =
  \frac{1}{2 \pi} e^{- \frac 1 2 (x_1^2)} \cdot \frac{1}{2 \pi} e^{- \frac 1 2 (x_2^2)}
 

Letzters ist die Dichtefunktion 2er unabhängiger Zufallsvariablen die Normalverteilt sind.

Axiom: D(1/(2 * %pi) * atan(x_2/x_1),x_1) a:=matrix[[%e^(-1/2 * (x_1^2 + x_2^2)) * x_1, %e^(-1/2 * (x_1^2 + x_2^2)) * x_2],[-(x_2)/(2 * %pi * (x_1^2 + x_2^2)), -(x_1)/(2 * %pi * (x_1^2 + x_2^2))]] determinant a

Accept-Reject, Simulation

Simulieren Sie von

f(x) = \mathrm{1}_{[0;2]}(x) \frac 3 4 x (2 - x), \,

indem Sie gleichmäßig auf der Fläche unter dem Graphen von f(x)\, simulieren.

Umschreiben Sie dazu den Graphen mit einer geeigneten Box, wie in Teil (a) des Accept/Reject Algorithmus aus der Vorlesung.

Plotten Sie ein Histogramm Ihrer generierten Zufallszahlen und zeichnen Sie zur Überprüfung die Dichtefunktion ein.

Algorithmus:

1. generiere (y,u) \sim \mathrm{unif}([a,b] \times [0,M]) \,

> a=0
>  b=2
>  M=1
>  y=runif(10000,a,b)
>  u=runif(10000,0,M)</pre>
  
 2. wenn f(y) \ge u \,, setze x := y\,, sonst gehe 1.
  
 
REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/f68c7e51b967d80b3aef615ebdb432da92ca0fd0_%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(x) { if(x >= 0 && x <= 2) { 3/4 * x * (2-x) } else 0 }
>  x = y[f(y) >= u]
Error: object 'y' not found
Execution halted
in
pdf(rpdf)
f = function(x) { if(x >= 0 && x <= 2) { 3/4 * x * (2-x) } else 0 }
x = y[f(y) >= u]
hist(x,main="f(x)", freq=FALSE, xlim=c(-.5,2.5),ylim=c(0,.8))

x = seq(-.5,2.5,length.out=1000)
lines(x,sapply(x,f),col=2)

Accept-Reject, aus der VO

Seien f(x)\, und g(x)\, Dichtefunktionen mit f(x) \le Mg(x) \forall x \in \mathbb{R}, M > 0\,. Nehmen Sie au\erdem an, dass eine Methode bekannt ist, um von g(x)\, zu simulieren.

Geben Sie einen Algorithmus an, um mithilfe von g(x)\, von f(x)\, zu simulieren, dh vervollständigen Sie Punkt (b) aus der Vorlesung.

1. Simuliere x\, nach g(x)\,

2. Simuliere y\, als gleichverteilt im Interval [0, M g(x)]\,

2. wenn f(y) < u \,, setze x := y\,, sonst gehe 1.