Wir setzen zuerst wieder die Hyperprameter fest:

ma = mb = 0
va = vb = 10^4

und Startwerte:

alpha = 0
beta = 0

Vektoren, in denen wir die Ziehungen abspeichern:

nriter=50000
alpha.sample = beta.sample = vector(length=nriter)

und noch einige Vorberechungen

ysum <- sum(y)
yxsum <- sum(x*y)

Wir benutzen hier als Vorschlagsverteilung (engl. proposal) die Random-Walk-Verteilung. Mit \(k\) die Iteration ist dies die zweidimensionale Normalverteilung

\[ (\alpha^*,\beta^*)\sim N_2((\alpha^{k-1},\beta^{k-1})),V) \]

also mit den “letzten” Ziehungen der Paramter als Erwartungswert und einer zu wählenden Kovarianzmatrix:

varcov <- matrix(c(1e-3, -5e-5, -5e-5, 1e-5), ncol=2)

Iteration

Eine einzelne Iteration schaut also wie folgt aus:

1. Ziehe aus der Vorschlagsverteilung (zweidimensionale Normalverteilung)

vorschlag <- LearnBayes::rmnorm(1,c(alpha,beta),varcov)
alpha.stern <- vorschlag[1]
beta.stern <- vorschlag[2]

2. Berechne die Akzeptanzwahrscheinlichkeit \(a\). Aus numerischen Gründen ist es wieder besser, zuerst die Log-Akzeptanzwahrscheinlichkeit zu berechnen. Wir berechnen dazu die Log-Posteriori-Dichte für den momentanen Zustand und für den Vorschlag, wobei wir Konstanten weglassen. Zur Erinnerung:

\[ \begin{align*} p(\alpha,\beta|y) \propto & \, \exp\left(-\frac{\alpha^2}{2v_\alpha}+\alpha\left(\frac{m_\alpha}{v_\alpha}+\sum y_t\right)\right)\\ &\cdot\exp\left(-\frac{\beta^2}{2v_\beta}+\beta\left(\frac{m_\beta}{v_\beta}+\sum x_ty_t\right)\right) \\&\cdot\exp\left(-\sum\exp(\alpha+\beta x_t)\right) \end{align*} \]

  log.dichte <- alpha*(ma/va+ysum) - 0.5*alpha^2/va +
      beta*(mb/vb+yxsum) - 0.5*beta^2/vb - sum(exp(alpha+beta*x))
  log.dichte.stern <- alpha.stern*(ma/va+ysum)  - 0.5*alpha.stern^2/va +
      beta.stern*(mb/vb+yxsum) - 0.5*beta.stern^2/vb - sum(exp(alpha.stern+beta.stern*x))
  akzeptanz <- exp(log.dichte.stern - log.dichte)

3. Ziehe \(u\) aus einer Gleichverteilung \(U[0,1]\). Ist \(u<a\) (die Akzeptanzwahrscheinlichkeit), dann akzeptiere den Vorschlag als neue Ziehung:

  if (akzeptanz>runif(1))
  {
    alpha<-alpha.stern
    beta<-beta.stern
  }

Was passiert, wenn der Vorschlag nicht akzeptiert wird? Dann wird die letzte Ziehung als neue Ziehung übernommen!

Testen Sie

Zuletzt gezogene Parameter
Vorgeschlagene Parameter
Akzeptanzwahrscheinlichkeit

Ziehen