Kapitel 5 Bayes’sche Ansätze

5.1 Einführung

In der Statistik gibt es zwei verschiedene Lehrmeinungen. Es sind dies

1. die Frequentisten und
2. die Bayesianer.

Alle bisher vorgestellten statistischen Konzepte, so zum Beispiel Least Squares, Maximum Likelihood, REML und BLUP stammen aus dem Lager der Frequentisten.

Die Unterschiede zwischen Frequentisten und Bayesianern bestehen hauptsächlich in

• deren Verständnis von Wahrscheinlichkeiten
• deren Unterteilung von Modell- und Datenkomponenten
• deren Techniken zur Schätzung von Parametern

Die folgende Tabelle gibt eine Übersicht über die Unterschiede.

5.2 Das Lineare Modell

Die Bayes’sche Art der Parameterschätzung soll an einem einfachen linearen Modell gezeigt werden. Angenommen, wir betrachten das Modell

\[\begin{equation} y_i = \beta_0 + \beta_1 x_{i1} + \epsilon_i \label{eq:BayLinMod} \end{equation}\]

wobei \(y_i\) die \(i\)-te Beobachtung einer Zielgrösse ist, \(\beta_0\) für den Achsenabschnitt steht, \(x_1\) eine erklärende Variable ist und \(\epsilon_i\) für den Restterm steht. Für den Restterm nehmen wir an, dass deren Varianz konstant gleich \(\sigma^2\) ist.

5.2.1 Bekannte und Unbekannte

Unter der Annahme, dass wir für die Zielgrösse \(y_i\) und die erklärende Variable \(x_1\) keine fehlenden Daten haben, dann machen wir als Bayesianer folgende Einteilung in bekannte und unbekannte Grössen.

und als bekannte Grössen

5.2.2 Vorgehen bei Parameterschätzung

Bayesianer basieren Schätzungen von unbekannten Grössen auf der sogenannten a posteriori Verteiung der unbekannten Grössen gegeben die bekannten Grössen. Die a posteriori Verteilung wird mithilfe des Satzes von Bayes aufgrund der a priori Verteilung der unbekannten und aufgrund der Likelihood berechnet.

Die Bezeichnungen “a priori” und “a posteriori” beziehen sich immer auf den Zeitpunkt der Beobachtung der analysierten Daten. Die jeweiligen Verteilungen quantifizieren den Informationsstand zu den Unbekannten um jeweiligen Zeitpunkt. Dieses Konzept soll anhand der folgenden Grafik verdeutlicht werden.

Für unser Beispiel des einfachen linearen Modells, definieren wir zuerst den Vektor \(\mathbf{\beta}\) als

\[\beta = \left[\begin{array}{c} \beta_0 \\ \beta_1 \end{array} \right].\]

Die Beobachtungen \(y_i\) fassen wir ebenfalls zu einem Vektor \(y\) zusammen. Für den Moment nehmen wir an, dass \(\sigma^2\) bekannt sei. Eine Bayes’sche Parameterschätzung für \(\beta\) basiert dann auf der a posteriori Verteilung \(f(\beta | y, \sigma^2)\) der Unbekannten \(\beta\) gegeben die Bekannten \(y\) und \(\sigma^2\). Diese a posteriori Verteilung lässt sich mit dem Satz von Bayes, wie folgt berechnen

\[\begin{eqnarray} f(\beta | y, \sigma^2) & = & \frac{f(\beta, \sigma^2, y)}{f(y, \sigma^2)} \nonumber \\ & = & \frac{f(y | \beta, \sigma^2)f(\beta)f(\sigma^2)}{f(y, \sigma^2)} \nonumber \\ & \propto & f(y | \beta, \sigma^2)f(\beta)f(\sigma^2) \label{LinModAPostProb} \end{eqnarray}\]

In Gleichung () konnten wir die a posteriori Verteilung \(f(\beta | y, \sigma^2)\) als Produkt der a priori Verteilungen (\(f(\beta)\) und \(f(\sigma^2)\)) der unbekannten Grössen \(\beta\) und \(\sigma^2\) und der Likelihood \(f(y | \beta, \sigma^2)\) ausdrücken. Der Faktor \(f(y, \sigma^2)^{-1}\) (Term im Nenner) entspricht der sogenannten Normalisierungskonstanten und ist nicht weiter von Interesse. Somit wird die a posteriori Verteilung oft als Proporzionalitätsbeziehung angegeben.

Die a posteriori Verteilung \(f(\beta | y, \sigma^2)\) ist in vielen Fällen nicht explizit darstellbar. Das war lange ein Problem, welches die Anwendung von Bayes’schen Analysen sehr einschränkte. Zwei Entwicklungen haben dieses Problem beseitigt.

1. In seinem Paper (Besag (1974)) zeigte Julian Besag, dass jede posteriori Verteilung durch eine Serie von Zufallszahlen aus den voll-bedingten Verteilungen bestimmt ist. Für unser Beispiel lauten die voll-bedingten Verteilungen: Bedingte Verteilung von \(\beta_0\) gegeben alle anderen Grössen: \(f(\beta_0 | \beta_1, \sigma^2, y)\), bedingte Verteilung von \(\beta_1\) gegeben alle anderen Grössen: \(f(\beta_1 | \beta_0, \sigma^2, y)\).
2. Die Entwicklung von effizienten Pseudo-Zufallszahlen-Generatoren auf dem Computer.

5.3 Gibbs Sampler

Die Umsetzung der beiden oben aufgelisteten Punkte führt zu einer Prozedur, welche als Gibbs Sampler bezeichnet wird. Wenden wir den Gibbs Sampler auf einfaches lineares Regressionsmodell an, dann resultiert das folgende Vorgehen bei der Analyse. Unabhängig vom verwendeten Modell läuft die Konstruktion einer Gibbs Sampling Prozedur immer in den folgenden Schritten ab. Diese Schritte können für die meisten Analysen wie ein Kochbuchrezept verwendet werden.

1. Bestimmung der a priori Verteilungen für die unbekannten Grössen.
2. Bestimmung der Likelihood
3. Bestimmung der voll-bedingten Verteilungen

5.3.1 A priori Verteilungen

In unserem Bespiel handelt es sich dabei um \(f(\beta)\) und \(f(\sigma^2)\). In den meisten Fällen, wenn man das erste Mal eine bestimmte Art von Daten analysisern soll, empfielt es sich eine sogenannte uniformative a priori Verteilung zu wählen. Eine uninformative a priori Verteilung bedeutet einfach, dass deren Dichtewert überall gleich, also eine Konstante ist. Wenden wir zum Beispiel für die Unbekannte \(\beta\) eine uninformative a priori Verteilung an, dann bedeutet das, dass wir \(f(\beta) = c\).

Alternativ zu der uniformativen a priori Verteilung gibt es auch a priori Verteiungen für bestimmte unbekannte Grössen, welche als de-facto Standard aktzeptiert sind. Ein Bespiel dafür ist die a priori Verteilung der unbekannten Restvarianz, welche üblicherweise als Inverse-Chi-Quadrat Verteilung angenommen wird.

5.3.2 Likelihood

Die Likelihood ist wie bei den Frequentisten als begingte Verteilung (\(f(y | \beta, \sigma^2)\)) der Daten \(y\) gegeben die Parameter (\(\beta\) und \(\sigma^2\)). Falls keine Daten fehlen, dann ist die Bayes’sche Likelihood und die frequentistische Likelihood gleich.

5.3.3 Vollbedingte Verteilungen

Mit vollbedingten Verteilungen ist gemeint, dass für jede unbekannte Grösse die bedingte Verteilung gegeben alle anderen Grössen bestimmt wird. In unserem Bespiel des linearen Regressionsmodells haben wir zwei unbekannte Grössen \(\beta_0\) und \(\beta_1\). Somit haben wir auch zwei vollbedingte Verteilungen. Unter der Annahme, dass unsere Daten (\(y\)) einer Normalverteilung folgen, resultieren die folgenden vollbedingten Verteilungen.

Aufgrund von Berechnungen, welche hier nicht gezeigt sind, können wir die oben aufgelisteten vollbedingten Verteilungen bestimmen. Die entsprechenden Verteilungen sind in der Kolonnen ganz rechts, welche mit “resultierende Verteilung” überschrieben ist, aufgelistet. Dabei steht \(\mathcal{N}()\) für die Normalverteilung. Für die Erwartungswerte und Varianzen wird das Modell in Gleichung () leicht umformuliert.

\[\begin{equation} \mathbf{y} = \mathbf{1}\beta_0 + \mathbf{x}\beta_1 + \mathbf{\epsilon} \label{eq:BayLinModReform} \end{equation}\]

Aus dem obigen Modell bilden wir ein neues Modell, welches auf der rechten Seite der Gleichung nur von \(\beta_0\) und \(\mathbf{\epsilon}\) abhängt. Da wir wissen, dass die Verteilung der Least Squares Schätzer eine Normalverteilung ist, werden wir diese für die Bestimmung der vollbedingten Verteilungen verwenden.

\[\begin{equation} \mathbf{w}_0 = \mathbf{1}\beta_0 + \mathbf{\epsilon} \label{eq:BayLinModW0} \end{equation}\]

wobei \(\mathbf{w}_0 = \mathbf{y} - \mathbf{x}\beta_1\). Aufgrund des Modells in Gleichung () können wir den Least Squares Schätzer für \(\beta_0\) aufstellen. Dieser lautet:

\[\begin{equation} \hat{\beta}_0 = (\mathbf{1}^T\mathbf{1})^{-1}\mathbf{1}^T\mathbf{w}_0 \label{eq:Beta0LsEst} \end{equation}\]

Die Varianz des Least Squares Schätzers für \(\beta_0\) lautet:

\[\begin{equation} var(\hat{\beta}_0) = (\mathbf{1}^T\mathbf{1})^{-1}\sigma^2 \label{eq:VarBeta0LsEst} \end{equation}\]

Analog zu \(\beta_0\) berechnen wir den Least Squares Schätzer für \(\beta_1\) und dessen Varianz.

\[\begin{equation} \hat{\beta}_1 = (\mathbf{x}^T\mathbf{x})^{-1}\mathbf{x}^T\mathbf{w}_1 \label{eq:Beta1LsEst} \end{equation}\]

wobei \(\mathbf{w}_1 = \mathbf{y} - \mathbf{1}\beta_0\)

\[\begin{equation} var(\hat{\beta}_1) = (\mathbf{x}^T\mathbf{x})^{-1}\sigma^2 \label{eq:VarBeta1LsEst} \end{equation}\]

5.3.4 Umsetzung des Gibbs Samplers

Der Gibbs Sampler wird durch wiederholtes ziehen von Zufallszahlen aus den oben angegebenen vollbedingten Verteilungen umgesetzt. Das heisst, wir setzen für alle unbekannten Grössen sinnvolle Startwerte ein. Für \(\beta_0\) und \(\beta_1\) wählen wir \(0\) als Startwert. Dann berechnen wir den Erwartungswert und die Varianz für die vollbedingte Verteilung von \(\beta_0\). Aus dieser Verteilung ziehen wir einen neuen Wert für \(\beta_0\). In einem zweiten Schritt berechnen wir den Erwartungswert und die Varianz für die vollbedingte Verteilung von \(\beta_1\), wobei wir für \(\beta_0\) schon den neuen Wert einsetzen. Aus der Verteilung für \(\beta_1\) ziehen wir einen neuen Wert für \(\beta_1\). Danach beginnen wir die Schritte wieder bei \(\beta_0\). Diese Schrittabfolge wiederholen wir \(10000\) mal und speichern alle gezogenen Werte für \(\beta_0\) und \(\beta_1\). Die Bayes’schen Parameterschätzungen entsprechen dann den Mittelwerten der gespeicherten Werte.

Der folgende R-Codeblock soll die Umsetzung des Gibbs Samplers für \(\beta_0\) und \(\beta_1\) als Programm zeigen. Der Einfachheit halber wurde \(\sigma^2\) konstant \(\sigma^2=1\) angenommen.

# ### Startwerte für beta0 und beta1
beta <– c(0, 0)
# ### Bestimmung der Anzahl Iterationen
niter <– 10000
# ### Initialisierung des Vektors mit Resultaten
meanBeta <– c(0, 0)
for (iter in 1:niter) {
  # Ziehung des Wertes des Achsenabschnitts beta0
  w <– y - X[, 2] * beta[2]
  x <– X[, 1]
  xpxi <– 1/(t(x) %*% x)
  betaHat <– t(x) %*% w * xpxi
  # ### neue Zufallszahl fuer beta0
  beta[1] <– rnorm(1, betaHat, sqrt(xpxi))
  # Ziehung der Steigung beta1
  w <– y - X[, 1] * beta[1]
  x <– X[, 2]
  xpxi <– 1/(t(x) %*% x)
  betaHat <– t(x) %*% w * xpxi
  # ### neue Zufallszahl fuer beta1
  beta[2] <– rnorm(1, betaHat, sqrt(xpxi))
  meanBeta <– meanBeta + beta
}
# ### Ausgabe der Ergebnisse
cat(sprintf("Achsenabschnitt = %6.3f \n", meanBeta[1]/iter))
cat(sprintf("Steigung = %6.3f \n", meanBeta[2]/iter))

References

Besag, J. 1974. “Spatial Interaction and the Statistical Analysis of Lattice Systems.” Journal of the Royal Statistical Society. Series B (Methodological), no. 36: 192–236. http://www.jstor.org/stable/2984812.