Zur Erinnerung, konjugierte Prioris im eindimensionalen waren:

• \(\mu\sim \text{N}(\mu_0,\sigma^2_0)\)
• \(\tau\sim \text{Ga}(a,b)\)

Gemeinsame Priori

Als ersten Ansatz folgende Idee: Wir wollen die konjugierten Prioris aus dem eindimensionalen benutzen und gehen außerdem *a priori** von Unabhängigkeit der Parameter aus. Die gemeinsam Priori-Dichte ist dann das Produkt der einzelnen

\[ p(\mu,\tau)=p(\mu)\cdot p(\tau) \]

Posteriori

Die Posteriori lautet dann bis auf Konstanten:

\[\begin{eqnarray*} p(\mu,\tau|x) &\propto& \exp\left(-\frac{1}{2\sigma_0^2} (\mu-\mu_0)^2\right)\cdot\\ &\cdot&\tau^{n/2}\exp\left(-\tau/2\sum_{i=1}^n(x-\mu)^2\right)\cdot\tau^{a-1}\exp(-b\tau) = \\ &=&\tau^{n/2+a-1} \exp\left(-\frac{1}{2\sigma_0^2} (\mu-\mu_0)^2)-\tau/2\sum_{i=1}^n(x-\mu)^2-b\tau\right) \end{eqnarray*}\]

Diese Posteriori hat nicht die Form der Priori. Die beiden Parameter \(\mu\) und \(\tau\) sind a posteriori voneinander abhängig, a priori aber nicht.

Das Produkt der beiden konjugierten Priori-Dichten ergibt also keine konjugierte Priori! Vielmehr ist die entstandene Posteriori keine bekannte zweiparametrische Verteilung.

Wir werden später zu diesem Ansatz zurückkommen, versuchen aber erstmal, eine konjugierte zweidimensionale Priori zu finden.

Für \(\tau\) geben wir eine marginale Priori-Verteilung vor:

\[ \tau \sim \text{Ga}(a,b) \]

Wir übernehmen also die konjugierten Verteilungen (Normalverteilung bzw. Gammaverteilung), aber nehmen keine Unabhängigkeit der Parameter an.

Die Begründung für diesen Ansatz ist wie folgt: Bei größerer Varianz (kleinerer Präzision) in den Daten brauchen wir auch größere Varianz (kleinere Präzision) in der Priori von \(\mu\), da sonst die Priori zu informativ wird.

Gemeinsame Priori

Die gemeinsame Priori \(p(\mu, \tau)\) ergibt sich nach der Definition der bedingten Dichte ebenfalls als Produkt der bedingten Dichte von \(\mu\) gegeben \(\tau\) und der marginalen Dichte \(\tau\):

\[ p(\mu, \tau) = p(\mu|\tau)p(\tau) \]

Die gemeinsame Priori-Dichte hat die Form

\[ p(\mu,\tau) = \frac{b^a \sqrt{\lambda}}{\Gamma(a)\sqrt{2\pi}} \tau^{a-\frac{1}{2}}\exp\left(-b\tau\right)\exp\left( -\frac{ \lambda \tau (\mu- \mu_0)^2}{2}\right) \]

Konjugierte Posteriori

Diese sogenannte Normal-Gamma-Verteilung mit Parametern \((\mu_0,\lambda,a,b)\) ist in der Tat die konjugierte Priori zur zweiparametrischen Normalverteilung. Als Posteriori ergibt sich:

\[ p(\mu,\tau|x) = \text{NormalGamma}\left(\frac{\lambda \mu_0 + n \bar{x}}{\lambda + n}, \lambda + n, a+\frac{n}{2}, b+ \frac{1}{2}\left(n s^2 + \frac{\lambda n (\bar{x} - \mu_0 )^2}{\lambda +n} \right) \right) \]

mit \(n\) der Anzahl an Beobachtungen, \(\bar{x}= \frac{1}{n}\sum_{i=1}^n x_i\) dem Mittelwert der Beobachtungen und \(s^2 = \frac{1}{n} \sum_{i=1}^n(x_i-\bar{x})^2\) der empirischen Varianz.

Interpretation der Priori-Parameter

Die Priori-Parameter können wie folgt interpretiert werden:

• \(\mu_0\) ist der Priori-Erwartungswert von \(\mu\)
• \(a/b\) ist der Priori-Erwartungswert von \(\tau\)
• \(a/b^2\) ist die Priori-Varianz von \(\tau\)
• \(\lambda\) ist die “Anzahl von Beobachtungen” in der Priori, vergleiche dazu das Binomialmodell

Für \(\lambda \to 0\) geht die bedingte Priori von \(\mu\) gegen die Gleichverteilung auf ganz \(\mathbb{R}\).

Posteriori-Schätzer

Posteriori-Erwartungswert und Posteriori-Modus lassen sich direkt aus den Parametern ablesen.

Der Posteriori-Erwartungswert ist

\[ E(\mu,\tau|x) = \left(\frac{\lambda \mu_0 + n \bar{x}}{\lambda + n}, \frac{a+\frac{n}{2}}{b+ \frac{1}{2}\left(n s^2 + \frac{\lambda n (\bar{x} - \mu_0 )^2}{\lambda +n} \right)}\right) \]

Der Posteriori-Modus (Maximum-A-Posteriori-Schätzer) ist

\[ (\hat{\mu},\hat{\tau})_{MAP}= \left(\frac{\lambda \mu_0 + n \bar{x}}{\lambda + n}, \frac{a+\frac{n-1}{2}}{b+ \frac{1}{2}\left(n s^2 + \frac{\lambda n (\bar{x} - \mu_0 )^2}{\lambda +n} \right)}\right) \]

Im vorliegenden Fall haben wir drei unterschiedliche Zufallsgrößen: \(\mu\), \(\tau\) und die Daten \(x\).

Betrachten wir die bedingte Verteilung von \(\mu\) gegeben \(\tau\) und \(x\), dann nennen wir dies bedingte Posteriori des Parameters \(\mu\). Posteriori bezeichnet immer die Bedingung auf die Beobachtungen \(x\), der Zusatz bedingt bezieht sich hier also darauf, dass wir bezüglich eines anderen Parameters, hier also \(\tau\) bedingen. Alle Verteilungen, die wir im Folgenden betrachten, sind immer gegeben der Daten \(x\).

Bedingte Posteriori

Nach der Definition der bedingten Dichte gilt

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

Nach dem Satz von Bayes gilt weiterhin

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

und da \(\mu\) und \(\tau\) a priori unabhängig:

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

Gegeben \(\tau\) ist also die Normalverteilungspriori quasi die konjugierte Priori zur bedingten Priori von \(\mu\). Die Herleitung kennen wir aus dem eindimensionalen Fall, bei dem ja \(\tau\) gegeben war.

Semikonjugierte Priori

Wie bereits beschrieben ist die Unabhängigkeitspriori keine konjugierte Priori. Allerdings ist die Priori von \(\mu\) eine semikonjugierte Priori zur bedingten Posteriori.

Definition:

Eine Familie \(\mathcal{F}\) von Verteilungen auf \(\Theta\) heißt semikonjugiert wenn für jede Priori \(p(\theta)\) auf \(\mathcal{F}\) die vollständig bedingte Posteriori \(p(\theta_i|\theta_1,\ldots\theta_{i-1},\theta_{i+1},\ldots,\theta_p;x)\) ebenfalls zu \(\mathcal{F}\) gehört.

Hier ist die bedingte Posteriori-Dichte von \(\mu\) gegeben \(\tau\) also die Normalverteilung. Die bedingte Posteriori hilft uns aber erstmal nicht weiter, weil wir den Parameter \(\tau\) ja nicht kennen.

Einfacher ist es, die Definition der bedingten Dichte zu benutzen und umzustellen:

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

Durch einsetzen erhalten wir:

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

Die einzelnen Terme sind bekannt

• \(f(x|\mu,\tau)\), die Datendichte einer Normalverteilung N(\(\mu,1/\tau)\)
• \(p(\tau)\) die Prioridichte einer Gammaverteilung Ga(\(a,b\))
• \(p(\mu|\tau;x)\) die Dichte einer Normalverteilung N(\(\tilde{\mu},1/\tilde{\tau}\)) mit \(\tilde\tau=\frac{n}{\sigma^2}+\frac{1}{\sigma_0^2}\) und \(\tilde{\mu}=\tilde{\tau}\cdot\left(\frac{\sum_{i=1}^nx_i}{\sigma^2}+\frac{\mu_0}{\sigma^2_0}\right)\) (siehe einparametrische Normalverteilung mit bekannter Varianz).

CDF-Sampler

Hier können wir den CDF-Sampler benutzen. CDF steht für cumulative distribution function, also (kumulative) Verteilungsfunktion.

Idee dabei ist, den Träger der Verteilung zu diskretisieren.

• Diskretisiere den Träger der zu simulierenden Verteilung in eine Menge von \(N\) Punkten \(x_1 \leq \ldots \leq x_N\).
• Evaluiere die bis auf Proportionalität bekannte Dichte an \(x_1 \leq \ldots \leq x_N\), um Werte \(f_1,\ldots,f_N\) zu erhalten.
• Schätze die Proportionalitätskonstante \(c\) über \(c = f_1 + \ldots + f_N\).
• Ziehe Zufallszahlen aus \(x_1 \leq \ldots \leq x_N\) gemäß den Wahrscheinlichkeiten \(f_1/c,\ldots,f_N/c\).

Die Implementierung des CDF-Samplers in R finden Sie hier

Mittels des CDF-Sampler, wie mit anderen Monte Carlo-Methoden, erhalten wir Ziehungen aus der Posteriori. Aus diesen können wir die Posteriori-Dichte der Parameter \(\mu\) und \(\tau\) schätzen:

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

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

Aus den Posteriori-Ziehungen lassen sich leicht Punkt- und Intervallschätzer gewinnen:

mean(mu) # Posteriori-Erwartungswert

## [1] 99.92407

median(mu) # Posteriori-Median

## [1] 99.92

quantile(mu, c(0.025, 0.975)) # symmetrisches 95%-KI

##   2.5%  97.5% 
##  98.26 101.55

mean(tau) # Posteriori-Erwartungswert

## [1] 0.184269

median(tau) # Posteriori-Median

## [1] 0.1724794

quantile(tau, c(0.025, 0.975)) # symmetrisches 95%-KI

##      2.5%     97.5% 
## 0.0645881 0.3681661

cor(mu, tau) # Posteriori-Korrelation

## [1] 0.01882033

Beachten Sie dabei: Wir beschränken uns auf Posteriori-Erwartungswert und Posteriori-Median, das sich diese leicht berechnen lassen. Ebenso ein symmatrisches Kredibilitätsintervall.

Eigentlich interessiert uns aber nicht \(\tau\), sondern die Varianz \(\sigma^2\). Die Ziehungen aus der Posterioriverteilung von \(\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))

Punkt- und Intervallschätzer von einzelnen Parametern erhalten wir also immer aus den marginalen Posteriori-Verteilungen.

Ist nun der Mittelwert des Geschichts der Schokohasen kleiner als 101 Gramm. Und ist die Varianz größer als 0.5? Die Posteriori-Wahrscheinlichkeiten schätzen wir über die relativen Häufigkeiten dieser Fälle unter unseren Ziehungen. Dabei ist \(m\) die Anzahl der Ziehungen.

sum(mu<101)/m

## [1] 0.9137

sum(sigma2>0.5)/m

## [1] 1

Die Posteriori-Wahrscheinlichkeit, dass \(\mu\) 101 Gramm oder mehr ist, liegt also gerade mal bei 8.6 %: eine geringe Wahrscheinlichkeit, aber eventuell ist das Ergebnis doch einfach auf die zufällige Auswahl der Schokohasen zurückzuführen.

Die Posteriori-Wahrscheinlichkeit, dass \(\sigma^2\) größer als 0.5 ist, ist dagegen 100 % (in allen Ziehungen ist \(\sigma^2\) größer als 0.5). Wir sind uns also sehr sicher, dass die Varianz in den Gewichten zu hoch ist.

(Zur Transparenz: es handelt sich nicht um reale Messungen, die Gewichte wurden aus einer N(101, 2)-Verteilung simuliert.)

Alternativ lassen sich die marginalen Posterioris beider Parameter über Diskretisierung approximieren und man kann auf das Ziehen komplett verzichten.

Dabei kann man folgenden Weg gehen: Wir hatten oben \(p(\mu|x)\) hergeleitet. Nun gilt

\[ p(\tau|x)=\int p(\mu,\tau|x)d\mu = \int p(\tau|\mu;x)p(\mu|x) d\mu \]

Für jedes feste \(\tau\) ist

\[ E_{\mu|x}(p(\tau|\mu;x)) = \int p(\tau|\mu;x)p(\mu|x) d\mu = p(\tau|x) \]

Die marginale Posteriori-Dichte an einem Punkt \(\tau\) erhalten wir also als Erwartungswert der bedingten Dichte am Punkt \(\tau\), wobei der Erwartungswert bezüglich der marginalen Posteriori von \(\mu|x\) gebildet wird.

Da wir \(\mu\) diskretisiert haben, können wir den Erwartungswert leicht berechnen (siehe oben). Dies machen wir auf einem Gitter für \(\tau\):

Priori

Für die Konstruktion von Prioris haben wir zwei Wege benutzt:

• Unabhängigkeit der Parameter, gemeinsame Priori ist dann Produkt der einzelnen marginalen Prioris.

• Bedingte Priori eines Parameters gegeben einem anderen. Dann ergibt sich die gemeinsame Prior als Produkt von bedingter und marginaler Priori.

• Beide Ansätze lassen sich auf mehr als zwei Dimensionen verallgemeinern.

• Später werden wir noch Modelle sehen, bei denen wir Abhängigkeiten in die Priori-Verteilungen zulassen.

Posteriori

Die Berechnung der Posteriori erfolgte im letzten Ansatz nicht mehr analytisch, sondern numerisch. Zwei Ansätze haben wir dabei erfolgt:

• Monte Carlo-Verfahren: Wir ziehen aus der gemeinsamen Posteriori-Verteilung, über das Gesetz der großen Zahlen erhalten wir dann Punktschätzer, Intervallschätzer und Posteriori-Wahrscheinlichkeiten.
• Numerische Approximation der Posteriori. Hier vor allem über Diskretisierung der marginale Posterioris.

Beide Ansätze werden in höherdimensionalen Modellen auch angewandt.