Zur Erinnerung:

\[ p(\mu|x)\propto\frac{f(x|\mu,\tau)p(\mu)}{p(\tau|\mu,x)} \]

• \(f(x|\mu,\tau)\), die Datendichte einer Normalverteilung N(\(\mu,1/\tau)\)
• \(p(\mu)\) die Prioridichte (N\((\mu_0,\sigma^2_0)\))
• \(p(\tau|\mu;x)\) die Dichte einer Gammaverteilung \(Ga\left(a+\frac{n}{2},b+\frac{1}{2}\sum_{i=1}^n(x_i-\mu)^2\right)\).

Unsere Schokohasen hatten folgende Gewichte:

## 100.37 103.04 99.34 101.5 97.44 100.26 99.19 100.12 99.76 98.51

Setzen wir erstmal die Parameter der Prioris, die sogenannten Hyperparameter:

a<-b<-1
sigma20<-10^9
mu0<-0

Dann brauchen wir eine Funktion, die uns die marginale Dichte von \(\mu\) bis auf die Konstante berechnet. Wir benutzen dazu die bekannten Dichten. Weil Dichten oft sehr klein sind, ist es numerisch besser, mit den Log-Dichten zu arbeiten.

p.mu.x<-function(mu, tau, x, a, b, sigma20)
{
  n <- length(x)
  logf <- 0
  for (i in 1:n)
  logf <- logf + dnorm(x[i], mu, sqrt(1/tau), log=TRUE)
  logf <- logf + dnorm(mu, mu0, sqrt(sigma20), log=TRUE)
  aa <- a+n/2
  for (i in 1:length(mu))
    {
    bb <- b+0.5*sum((x-mu[i])^2)
    logf[i] <- logf[i] - dgamma(tau, aa, bb, log=TRUE)
  }
  logf <- logf - mean(logf)
  return(logf)
}

• Nun können wir auf einem Gitter von \(\mu=95.00, 95.01, 95.02, \ldots, 104.99, 105.00\) die marginale Dichte berechnen.
• Unsere Formel und unsere Funktion hängt noch von \(\tau\) ab.
• Da \(p(\mu|x)\) aber nicht von \(\tau\) abhängt, können wir einen beliebigen Wert für \(\tau\) wählen! (Extreme Werte könnten aber zu numerischen Probleme führen).

mu.x<-seq(95,105,length=1001)
tau <- 42 # Hier kann beliebiger Wert > 0 stehen
log.p<-p.mu.x(mu.x, tau, gewicht, a, b, sigma20)

Nun kommt numerische Integration. Durch unsere Diskretisierung wird das Integral zur Summe. Die Summe über die Wahrscheinlichkeiten muss 1 sein:

expint<-sum(exp(log.p))
p<- exp(log.p-log(expint))
plot(mu.x,p,type="l")

Wir haben die marginale Posteriori von \(\mu\) damit schon gut approximiert. Daraus lässt sich zum Beispiel der Erwartungswert berechnen:

\[ E(\mu|x)=\sum_\mu \mu\cdot p(\mu|x) \]

sum(mu.x*p)

## [1] 99.953

und die marginale Posteriori von \(\mu\) zeichnen:

Vorsicht: p ist hier die Wahrscheinlichkeit, das \(\mu\) in einem Intervall der Länge 0.01 liegt. Für die Dichte müssen wir daher p durch 0.01 teilen, damit insgesamt das Integral über die Dichte gleich 1 ist.

Nun ziehen wir \(m=10000\) mal aus der marginalen Posteriori von \(\mu\):

m<-10000
mu<-sample(mu.x,m,replace=TRUE,prob=p)

Der Mittelwert der Ziehungen sollte dem obig numerisch berechneten Erwartungswert entsprechen, zumindest in den ersten paar Stellen.

mean(mu)

## [1] 99.94348

sum(mu.x*p)

## [1] 99.953

Die Posteriori-Varianz von \(\mu\) ist hier höher als im eindimensionalen Fall, als wir \(\sigma^2\) als bekannt angenommen hatten.

var(mu)

## [1] 0.2716631

Wichtiger ist nun das Ziehen von \(\tau\), wobei wir die Ziehungen von \(\mu\) einsetzen:

bb<-vector(length=m)
for (i in 1:m)bb[i] <- b+0.5*sum((gewicht-mu[i])^2)
tau <- rgamma(m,a+n/2,bb)

Mittels Kerndichteschätzer können wir die marginalen Posterioris der beiden Parameter zeichnen:

plot(density(mu), main="", xlab=expression(mu), xlim=c(95,105))

plot(density(tau), main="", xlab=expression(tau))

und aus den Posteriori-Ziehungen die Punkt- und Intervallschätzer gewinnen:

## [1] 0.4636742

## [1] 0.4363333

##      2.5%     97.5% 
## 0.1598051 0.9329517

## [1] 0.01716948

Eigentlich interessiert uns aber nicht \(\tau\), sondern die Varianz \(\sigma^2\). Die Ziehungen für \(\sigma^2\) erhalten wir einfach durch \(\sigma^2=1/\tau\).

sigma2<-1/tau

Damit haben wir also die Posteriori von \(\sigma^2\) ohne weitere Rechnungen:

plot(density(sigma2), main="", xlab=expression(sigma^2), xlim=c(0,500))