Statistische Informationsverarbeitung - Übung 2

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.