Eine Idee, eine sinnvolle Vorschlagsverteilung zu konstruieren ist, die full conditional möglichst gut zu approximieren. Ist \(q\left(\theta^{*}|\theta^{(k-1)}\right)\approx f\left(\theta^{*}\right)\), dann wird die Akzeptanzwahrscheinlichkeit fast 1, wir sind also “fast” beim Gibbs-Sampler.

Gleichzeitig müssen wir gut aus der Vorschlagsverteilung ziehen können.

Im vorliegenden Fall ist die full conditional

\[ p(\eta_t|\alpha,\beta,\eta_{-t},y) \propto \exp(\eta_t y_t)\exp\left(-\exp(\eta_t)\right)\exp\left(-\frac{1}{2\sigma^2}\left(\eta_t-\alpha-\beta x_t\right)^2\right) \]

Der erste und der dritte Term passen zu einer Normalverteilung. Nur der Term \(\exp\left(-\exp(\eta_t)\right)\) passt nicht dazu. Benutzen wir die Taylorreihe

\[ f(x; a) = \sum_{n=0}^\infty \frac{f^{(n)}(a)}{n!} (x-a)^n \]

wobei \(f^{(n)}\) die \(n\)-te Ableitung ist und \(a\) ein Punkt, um den die Reihe entwickelt wird.

Setzen wir \(f(x)=\exp(\eta_t)\), dann bekommen wir

\[ \exp(\eta_t)\approx \frac{\exp(a)}{1}+ \frac{\exp(a)}{1}(\eta_t-a) +\frac{\exp(a)}{2}(\eta_t-a)^2 \]

und brechen nach \(n=2\) ab.

Setzen wir das in die full conditional von \(\eta_t\) ein erhalten (mit \(\underset{\sim}{\propto}\) approxmitativ proportional zu)

\[ \begin{align*} p(\eta_t|\alpha,\beta,y) \underset{\sim}{\propto}\, & \exp\left(\eta_t y_t-\eta_t\exp(a)-\frac{\exp(a)}{2}(\eta_t-a)^2-\frac{1}{2\sigma^2}\left(\eta_t-\alpha-\beta x_t\right)^2\right)\\ =\,&\exp\left(-\frac{1}{2}\eta_t^2(\sigma^{-2}+\exp(a))+\eta_t\left(y_t-\exp(a)+a\exp(a)+\sigma^{-2}(\alpha+\beta x_t)\right)\right) \end{align*} \]

Das ist eine Normalverteilung für \(\eta_t\):

• \(\eta_t\sim N\left(\frac{m}{v},\frac{1}{v}\right)\)

• \(m=y_t-\exp(a)+a\exp(a)+\sigma^{-2}(\alpha+\beta x_t)\)

• \(v=\sigma^{-2}+\exp(a)\)

• Diese Normalverteilung approximiert uns die full conditional von \(\eta_t|\cdot\)

• Diese Normalverteilung benutzen wir als Vorschlagsverteilung im Metropolis-Hastings-Schritt

• Die Verteilung gilt für alle \(t=1,\ldots,T\) und die \(\eta_t\) sind voneinander unabhängig.

\(\eta_t|\cdot=\eta_t|\alpha,\beta,\eta_{-t};y\) meint hier die Verteilung von \(\eta_t\) gegeben aller anderen Parameter und den Daten.

Wir setzen zuerst wieder die Hyperprameter fest:

ma = mb = 0
va = vb = 10^3
sigma2 = 10^{-3}

und Startwerte:

alpha = 10
beta = 1
epsilon = rep(0,T)
eta = alpha + beta*x +epsilon

Vektoren, in denen wir die Ziehungen abspeichern:

nriter=50000
alpha.sample = beta.sample = vector(length=nriter)
eta.sample = matrix(ncol = T, nrow=nriter)
epsilon.sample = matrix(ncol = T, nrow=nriter)

und noch einige Vorberechungen

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

Algorithmus

for (i in 1:nriter)
{
  # Ziehe beta
  v.beta <- 1/vb+x2sum/sigma2
  m.beta <- mb/vb+sum((eta-alpha)*x)/sigma2
  beta <- rnorm(1, m.beta/v.beta, sqrt(1/v.beta))

  # Ziehe alpha
  v.alpha <- 1/va+T/sigma2
  m.alpha <- ma/va+sum(eta-beta*x)/sigma2
  alpha <- rnorm(1, m.alpha/v.alpha, sqrt(1/v.alpha))
  
  # Vorschlag eta, als Vektor
  v <- 1/sigma2+exp(eta)
  m <- y-exp(eta)+eta*exp(eta)+(alpha+beta*x)/sigma2 
  eta.stern <- rnorm(T, m/v, sqrt(1/v)) 
  
  # Berechnung Log-Akzeptanzwahrscheinlichkeit
  v2 <- 1/sigma2+exp(eta.stern)
  m2 <- y-exp(eta.stern)+eta*exp(eta.stern)+(alpha+beta*x)/sigma2 
  log.akzeptanz <- eta.stern*y-exp(eta.stern)-0.5*(eta.stern-alpha-beta*x)^2/sigma2
  log.akzeptanz <- log.akzeptanz - (eta*y - exp(eta)-0.5*(eta-alpha-beta*x)^2/sigma2)
  log.akzeptanz <- log.akzeptanz - dnorm(eta.stern, m/v, sqrt(1/v), log=TRUE)
  log.akzeptanz <- log.akzeptanz + dnorm(eta, m2/v2, sqrt(1/v2), log=TRUE)

  # Akzeptanz
  ja <- exp(log.akzeptanz)>runif(T)   
  eta[ja] <- eta.stern[ja]
  
  # Abspeichern
  alpha.sample[i] <- alpha
  beta.sample[i] <- beta
  eta.sample[i,] <- eta
  epsilon.sample[i,] <- eta-alpha-beta
}

Wir verwerfen wieder 5000 Iterationen als burn-in.

alpha.sample<-alpha.sample[-(1:5000)]
beta.sample<-beta.sample[-(1:5000)]
epsilon.sample<-epsilon.sample[-(1:5000),]