Statistische Informationsverarbeitung - Übung 7

From StatWiki
Jump to: navigation, search

13. Random Walk Metropolis

Berechnen Sie den Erwartungswert von X \sim \Gamma(\alpha = 2.1, \beta = 1.5) \, mithilfe von Metropolis-Hastings, also durch

\operatorname{E}(X) \approx \frac 1 n \sum_{t=1}^n X_t \,

wobei X_t\, die entsprechende MH-Kette für \pi(x) \sim \Gamma(\alpha = 2.1, \beta = 1.5) \, mit einem Startpunkt X_0 = x_0\, Ihrer Wahl bezeichnet.

Wählen Sie n nach Ihrem Ermessen. Beachten Sie,dass Sie hier - im Gegensatz zum Irrfahrt-Beispiel - nur eine Markovkette laufen lassen müssen.

Allerdings können Sie zum Schluss auch über mehrere Ergebnisse mitteln.

Verwenden Sie die proposal distribution q(y|x) \propto e^{-\frac 1 2 (y-x)^2}\, also Y_t|X_t = x \sim N(x,1)\,, bzw. Y_t = X_t + \epsilon_t\, mit \epsilon_t \sim N(0,1)\,.

Dabei müssen Sie beachten, dass Y_t\, auch negativ sein kann, und \pi\, dann an diesem Wert 0 ist.

(Für die Implementierung bedeutet das, dass Sie ein Y_t < 0 \, immer verwerfen.)

Übrigens: Eine Zufallsvariable X \sim \Gamma(\alpha, \beta)\, hat Dichtefunktion

f_X(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}\, mit x > 0 \,

Für uns gilt also \pi(x) \propto x^{1.1} e^{-1.5x}\,.

Dieses Beispiel soll v.a. dazu dienen, den Metropolis-Hastings Algorithmus einmal anhand eines einfachen Beispiels zu implementieren (für das auch sehr grosse n möglich sind).

Wir können für die \Gamma\,-Verteilung den Erwartungswert natürlich auch explizit berechnen, es gilt \operatorname{E}(X) = \frac \alpha \beta \,. Sie können damit Ihr Ergebnis überprüfen.

> pdf(rpdf)
> pix = function(x) { x^1.1 * exp(-1.5 * x) }
> prop = function(x,y) { exp(-.5 * (y-x)^2) }
>
> mh = function(n) {
+         # path = rep(0,n)
+         # xt = path[1] = pix(1)
+         xt = s = pix(1)
+         # alot faster
+         rn = rnorm(n)
+         un = runif(n)
+         for(i in 2:n) {
+                 # yt = xt + rnorm(1)
+                 yt = xt + rn[i]
+                 # if(yt >= 0 && runif(1) < ((pix(yt) * prop(xt, yt))/(pix(xt) * prop(yt, xt)))) {
+                 # if(yt >= 0 && un[i] < ((pix(yt) * prop(xt, yt))/(pix(xt) * prop(yt, xt)))) {
+                 # Weil prop(x,y) == prop(y,x) ist, muss das nicht berechnet werden
+                 if(yt >= 0 && un[i] < (pix(yt) / pix(xt))) {
+                         xt = yt
+                 }
+                 # path[i] = xt
+                 s = s + xt
+         }
+         # path
+         s/n
+ }
>
> # mean(replicate(10,mh(1000)))
> # mean(mh(10000))
> # boxplot(replicate(100,mh(100000)))
>
> boxplot(replicate(10,mh(1000)))
> abline(h=2.1/1.5,col=2)