Zuerst setzen wir die Priori-Parameter (Hyper-Parameter) fest

## Priori-Parameter
ma = mb = 0
va = vb = 10^3
a = b = 10^{-3}

Einige Werte, die wir brauchen, können wir im Vorraus berechnen

a2 = a + n/2
x2 <- sum(x^2)

Wir bereiten Vektoren vor, in denen wir die Ziehungen abspeichern:

AnzahlIterationen = 1000
alpha.sample = vector(length=AnzahlIterationen)
beta.sample = vector(length=AnzahlIterationen)
sigma2.sample = vector(length=AnzahlIterationen)

Für den Beginn des Algorithmus brauchen wir Startwerte:

alpha = 0
beta = -10
sigma2 = 1

• Wir ziehen ja einen Parameter gegeben den anderen, bei der ersten Ziehung müssen die anderen also gegeben sein.
• Das heißt, dass wir uns einen Startwert eigentlich sparen können: den Parameter, denn wir zuerst ziehen.

Jetzt der Algorithmus. Wir ziehen abwechselnd aus den full conditionals, also der

• Normalverteilung für \(\alpha|\beta,\sigma^2\)
• Normalverteilung für \(\beta|\alpha,\sigma^2\)
• Invers-Gamma-Verteilung für \(\sigma^2|\alpha,\beta\)

Nachdem alle Parameter einmal gezogen wurden, speichern wir alle im vorbereiteten Vektor.

for (i in 1:AnzahlIterationen)
{
  v = 1/va + n/sigma2
  m = ma/va + sum(y-beta*x)/sigma2 
  alpha <- rnorm(1, m/v, sqrt(1/v))
  
  v = 1/vb + x2/sigma2
  m = mb/vb + sum(x*(y-alpha))/sigma2
  beta <- rnorm(1, m/v, sqrt(1/v))
  
  b2 <- b + 0.5*sum((y-alpha-beta*x)^2)
  sigma2 <- LearnBayes::rigamma(1, a2, b2)
  
  alpha.sample[i] <- alpha
  beta.sample[i] <- beta
  sigma2.sample[i] <- sigma2
}