Die Schweiz trat Ende des 19. Jahrhunderts in eine Periode ein, die als demographischer Übergang bezeichnet wird, d.h. ihre Fruchtbarkeit begann von dem für unterentwickelte Länder typischen hohen Niveau abzufallen. Francine Vanderwalle sammelte Daten über die Fruchtbarkeit und verschiedene Kovariablen in der Schweiz in dieser Zeit. Wir benutzten die Daten aus dem Jahr 1888 für die 47 französischsprachigen Bezirke.

Wir wollen die Abhängigkeit des Fertilitätsindex vom Anteil des Agrarsektors (genauer: Anteil der Arbeitsplätze (von Männern) in der Landwirtschaft untersuchen. Wir gehen von einem linearen Zusammenhang aus und benutzten lineare Regression.

Das übliche lineare Regressionsmodell für \(y\) mit einer Kovariablen \(x\) ist

\[ y_i=\alpha+\beta x_i+\epsilon_i, \quad i=1,\ldots,n \]

Dabei werden in der Regel folgende Annahmen an den Fehler \(\epsilon_i\) getroffen:

• E\((\epsilon_i)=0\)
• Var\((\epsilon_i)=\sigma^2\)
• Cov\((\epsilon_i,\epsilon_j)=0\)

Zusätzlich nimmt man oft an

• \(\epsilon_i\sim\)N\((0,\sigma^2)\)

Bayesianisches lineares Regressionsmodell

Für einen Bayesianischen Ansatz brauchen wir eine Datendichte, also eine Verteilungsannahme.

Aus dem linearen Modell \(y_i=\alpha+\beta x_i+\epsilon_i\) und der Annahme \(\epsilon_i\sim\)N\((0,\sigma^2)\) ergibt sich eine Verteilung für \(y\)

\[ y_i|\alpha,\beta, \sigma^2 \sim \text{N}(\alpha+\beta x_i,\sigma^2) \]

• Genauer gesagt ist dies die bedingte Verteilung von \(y_i\) gegeben den Parametern \(\alpha,\beta\) und \(\sigma^2\).
• In der Bayesianischen Statistik sind die Parameter ja Zufallsvariablen und haben eine Verteilung: Die Priori-Verteilung, die wir noch festlegen müssen, bzw. die Posteriori-Verteilung gegeben der Daten.

Die Datendichte aller \(y=(y_1,\ldots,y_n)\) gegeben den Parametern \(\theta=(\alpha,\beta,\sigma^2)\) ist damit:

\[ \begin{aligned} f(y|\theta) = & \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{1}{2\sigma^2}(y_i-\alpha-\beta x_i)^2\right)\\ = & \frac{1}{(2\pi\sigma^2)^{n/2}}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\alpha-\beta x_i)^2\right) \end{aligned} \]

Priori der Regressionsparameter

Bei der Wahl der Priori lehnen wir uns an der Wahl der Priori im Normalverteilungsmodell an. Dort hatten wir eine Normalverteilungspriori für den Erwartungswert benutzt. Nehmen wir an, dass

• \(\alpha\sim\text{N}(m_\alpha,v_\alpha^2)\)
• \(\beta\sim\text{N}(m_\beta,v_\beta^2)\)

und \(\alpha\) und \(\beta\) a priori unabhängig. Dann gilt

\[ (\alpha+\beta x_i)\sim\text{N}(m_\alpha+x_in_\beta, v_\alpha^2+x_i^2v_\beta^2) \]

• Wie wir später zeigen werden, handelt es sich hierbei um die semi-konjugierten Prioris für \(\alpha\) und \(\beta\).

• Für eine möglicht wenig informative Priori setzten wir die Priori-Varianzen \(v_\alpha^2\) und \(v_\beta^2\) sehr hoch (flache Priori) oder lassen sie gar gegen \(+\infty\) gehen. Dann haben wir uneigentliche Verteilungen:

• \(p(\alpha)\propto\)const.

• \(p(\beta)\propto\)const.

Priori der Varianz

Analog benutzen wir für die Varianz \(\sigma^2\) die (semi-)konjugierte Priori:

\[ \sigma^2 \sim IG(a,b), \]

oder, basierend auf der eindimensionalen Jeffreys’ Priori

\[ p(\sigma^2)\propto \sigma^{-2}, \]

welche wieder eine uneigentliche Priori ist.

Gemeinsame Priori

Zusätzlich gehen wir a priori davon aus, das \(\sigma^2\) von \(\alpha\) und \(\beta\) stochastisch unabhängig ist. Damit erhalten wir die gemeinsame Priori-Dichte aller Parameter:

\[ p(\alpha, \beta, \sigma^2)=p(\alpha)p(\beta)p(\sigma^2) \]

Die Annahme der Priori-Unabhängigkeit der Parameter ist Standard und wird implizit immer unterstellt, solange keine konkrete Abhängigkeit in die Parameter modelliert wird.

Posteriori

Leiten wir nun die Posteriori der Parameter gegeben der Daten her. Es gilt mit Daten \(y=(y_1,\ldots,y_n)\) und Parametern \(\theta=(\alpha,\beta,\sigma^2)\):

\[\begin{eqnarray*} p(\theta|y)&\propto&f(y|\theta)p(\theta)\\ &=& \prod_{i=1}^n f(y_i|\theta)p(\alpha)p(\beta)p(\sigma^2)\\ &\propto & \sigma^{-n}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\alpha-\beta x_i)^2\right)\cdot\\ &&\cdot v_\alpha^{-1} \exp \left(-\frac{1}{2v_\alpha^2}(\alpha-m_\alpha)^2\right)\\ &&\cdot v_\beta^{-1} \exp \left(-\frac{1}{2v_\beta^2} (\beta-m_\beta)^2\right)\\ &&\cdot \sigma^{-2(a+1)}\exp\left(b\sigma^{-2}\right) \end{eqnarray*}\]

Diese Posterioridichte ist wohl keine Standard-Verteilung mehr. Wie können wir die gemeinsame Posteriori erhalten, wie die marginalen Posterioris, aus denen wir Schätzer ableiten?

Im Folgenden benutzen wir den sogenannten Gibbs-Sampler. Hierbei handelt es sich wieder um einen Monte-Carlo-Algorithmus. Ziel ist es, Zufallszahlen aus der Posterioriverteilung zu ziehen und:

• den Posteriori-Erwartungswert über den Mittelwert der gezogenen Zufallszahlen schätzen;
• die Posteriori-Wahrscheinlichkeit eines Intervalls über die relative Häufigkeit zu schätzen, mit der die gezogenen Zufallszahlen in das Intervall fallen;
• die Posteriori-Dichte über ein Histogramm oder mittels des Kern-Dichte-Schätzers zu schätzen.

Idee

Die Grundidee des Gibbs-Samplers ist es, abwechselnd aus der bedingten (Posteriori-)Verteilung eines Parameters gegeben der anderen Parameter zu ziehen.

Algorithmus

Gegeben sei die Dichte \(f(\theta)\) mit \(\theta=(\theta_1,\ldots,\theta_p)\). Der Algorithmus des Gibbs-Samplers ist wie folgt:

1. Lege Startwerte \(\theta_i^{(0)}\) für alle \(i=1,\ldots,p\) fest.
2. Setze Iteration \(k=1\).
3. Ziehe einen Wert \(\theta_1^{(k)}\) gegeben \(\theta_2^{(k-1)},\ldots,\theta_p^{(k-1)}\) aus der Dichte \(f(\theta_1|\theta_2,\ldots,\theta_p)\).
4. Ziehe analog für \(j=2,\ldots,p\) einen Wert \(\theta_j^{(k)}\) gegeben \(\theta_1^{(k)},\ldots,\theta_{j-1}^{(k)},\theta_{j+1}^{(k-1)},\ldots\theta_p^{(k-1)}\) aus der Dichte \(f(\theta_j|\theta_{-j})\).
5. Erhöhe \(k\) um 1.
6. Wiederhole Schritte 3 bis 5 so oft wie nötig.

Dabei ist \(\theta_{-j}=\theta_1,\ldots,\theta_{j-1},\theta_{j+1},\ldots\theta_p)\), also alle Parameter außer \(\theta_j\). \(\theta_j^{(k)}\) steht für den Wert, den \(\theta_j\) in der Iteration \(k\) hat.

Einfaches Beispiel

Wir wollen aus einer zweidimensionalen Normalverteilung ziehen:

\[ x\sim N_2\left({0\choose 0},\left(\begin{array}{cc}1&0.5\\0.5&1\end{array}\right)\right) \]

Die bedingten Verteilungen sind wieder Normalverteilungen:

\[ x_1|x_2 \sim N(0.5\cdot x_2,1-0.5^2)\\ x_2|x_1 \sim N(0.5\cdot x_1,1-0.5^2) \]

x1<-x2<-vector(length=100)
x2[1] <- 0 # Startwert von x2
# Den Startwert für x1 können wir und sparen!

for (i in 2:1000)
{
  x1[i]<-rnorm(1,0.5*x2[i-1],sqrt(.75))
  x2[i]<-rnorm(1,0.5*x1[i],sqrt(.75))
}

Damit erhalten wir folgende 1000 Ziehungen aus der zweidimensionalen Normalverteilung:

Die eindimensionalen marginalen Dichten können wir über Kerndichteschätzung erhalten:

d1<-density(x1)
d2<-density(x2)

Schliessliche können wir die Erwartungswert über die Mittelwerte schätzen, die marginalen Varianzen über die empirischen Varianzen und die Korrelation über die empirische Korrelation:

mean(x1)

## [1] -0.001457514

mean(x2)

## [1] -0.02283133

var(x1)

## [1] 0.9530554

var(x2)

## [1] 0.9788648

cor(x1,x2)

## [1] 0.5059267

• Die Ziehungen sind beim Gibbs-Sampler voneinander abhängig: \(\theta_1^{(k+1)}\) hängt von \(\theta_{-1}{(k)}\) ab, die wiederum von \(\theta_1^{(k)}\) abhängen. Damit hängt \(\theta_1^{(k+1)}\) also vom vorherigen Wert \(\theta_1^{(k)}\) ab.

• Trotz der Abhängigkeit der Ziehungen kann man zeigen, dass der Algorithmus Zufallszahlen aus der gemeinsamen Dichte \(f(\theta)\) zieht.

• Die Schätzung ist allerdings weniger effizient, als wenn wir unabhängige Ziehungen hätten. Dass heißt, wir müssen insgesamt öfter ziehen.

• Zu Beginn des Algorithmus hängen die Ziehungen zudem von den im Schritt 1 gewählten Startwerten ab. Wir müssen die Ziehungen am Anfang also ignorieren – die sogenannte burn-in-Phase.

• Nach dem burn-in können wir Erwartungwert, Varianz, Wahrscheinlichkeiten, Dichte etc. schätzen (trotz der Abhängigkeit der Ziehungen).

• Dafür braucht man eine allgemeinere Version des Gesetzes der großen Zahlen, den Ergodensatz.

Vollständig bedingte Posterioris

• Kommen wir zur Anwendung des Gibbs-Samplers in unserem Regressionsmodell.
• Wir brauchen für den Algorithmus die Verteilung - in unserem Fall die Posteriori-Verteilung - eines Parameter (\(\theta_i\)) gegeben den restlichen.
• Diese Verteilung nennen wir vollständig bedingte Posteriori oder englisch full conditional.

Wie können wir die full conditionals (man benutzt den Begriff sowohl für die Verteilung als auch für die Dichte) herleiten? Aus der Definition der bedingten Dichte ergibt sich:

\[ f(\theta_i|\theta_{-i})=\frac{f(\theta)}{f(\theta_{-i})}\propto f(\theta) \]

Die Proportionalität ergibt sich dabei, da \(f(\theta_{-i})\) nicht von \(\theta_i\) abhängt, also bezüglich \(\theta_i\) eine Konstante ist.

Die vollständige bedingte Dichte ist also einfach proportional zur gemeinsamen Dichte.

• Praktisch nehmen wir uns aus der Posterioridichte einfach die Terme, die von \(\theta_i\) abhängen (und hoffentlich kommt dann eine Verteilung heraus, aus der wir ziehen können).
• Um aus \(f(\theta_i|\theta_{-i})\) zu ziehen, können wir entweder:
  • aus \(f(\theta)\) auf den Kern einer (Standard-)Verteilung schließen
  • oder ein approximatives Verfahren wie den CDF-Sampler benutzen (ist in der Regel aber langsamer)

Full conditional von \(\beta\)

Die Posterioridichte von \(\theta=(\alpha,\beta,\sigma^2)\) ist

\[\begin{eqnarray*} p(\theta|y)&\propto & \sigma^{-n}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\alpha-\beta x_i)^2\right)\cdot\\ &\cdot& v_\alpha^{-1} \exp \left(-\frac{1}{2v_\alpha^2}(\alpha-m_\alpha)^2\right)\\ &\cdot& v_\beta^{-1} \exp \left(-\frac{1}{2v_\beta^2} (\beta-m_\beta)^2\right)\\ &\cdot& \sigma^{-2(a+1)}\exp\left(b\sigma^{-2}\right) \end{eqnarray*}\]

Nehmen wir nun die Terme, die von \(\beta\) abhängen:

\[\begin{eqnarray*} p(\beta|\alpha,\sigma^2;y) &\propto& \exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\alpha-\beta x_i)^2\right) \exp \left(-\frac{1}{2v_\beta^2} (\beta-m_\beta)^2\right)\\ &\propto& \exp\left(-\frac{1}{2\sigma^2}\left(\beta^2\sum_{i=1}^nx_i^2-2\beta\sum_{i=1}^n(x_i(y_i-\alpha))\right)-\frac{1}{2v_\beta^2} (\beta^2-2\beta m_\beta)\right)\\ &\propto& \exp\left(-\frac{1}{2}\beta^2\left(v_\beta^{-2}+\sigma^{-2}\sum_{i=1}^nx_i^2\right)+\beta\left(\frac{m_\beta}{v_\beta^2}+\sigma^{-2}\sum_{i=1}^n(x_i(y_i-\alpha))\right)\right) \end{eqnarray*}\]

Durch quadratische Ergänzung erhält man:

• \(\beta|\alpha,\sigma^2;y\) ist normalverteilt mit

• Inverser Varianz \(\tilde{v}_{\beta}^{-2}=v_\beta^{-2}+\sigma^{-2}\sum_{i=1}^nx_i^2\) und

• Erwartungwert \(\tilde{m}_\beta=\tilde{v}_\beta^{-2}\left(v_{\beta}^{-2}m_\beta +\sigma^{-2}\sum_{i=1}^n(x_i(y_i-\alpha))\right)\).

Full conditional von \(\alpha\)

Analog erhält man

\[ p(\alpha|\beta,\sigma^2;y) \propto \exp\left(-\frac{1}{2}\alpha^2(v_\alpha^{-2}+n\sigma^{-2})+\alpha\left(m_\alpha+\sum_{i=1}^n(y_i-\beta x_i)\right)\right) \]

Die Full Conditional ist also eine Normalverteilung mit

• Inverser Varianz \(\tilde{v}_{\alpha}^{-2}=v_\alpha^{-2}+n\sigma^{-2}\) und
• Erwartungwert \(\tilde{m}_\alpha=\tilde{v}_\alpha^{-1}\left(v_\alpha^{-2}m_\alpha +\sigma^{-2}\sum_{i=1}^n(y_i-\beta x_i)\right)\).

Full conditional von \(\sigma^2\)

Zuletzt leiten wir noch die Full Conditional von \(\sigma^2\) her:

Die Posterioridichte ist

\[\begin{eqnarray*} p(\theta|y)&\propto & \sigma^{-n}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\alpha-\beta x_i)^2\right)\cdot\\ &&\cdot v_\alpha^{-1} \exp \left(-\frac{1}{2v_\alpha^2}(\alpha-m_\alpha)^2\right)\\ &&\cdot v_\beta^{-1} \exp \left(-\frac{1}{2v_\beta^2} (\beta-m_\beta)^2\right)\\ &&\cdot \sigma^{-2(a+1)}\exp\left(b\sigma^{-2}\right) \end{eqnarray*}\]

Folgende Terme hängen von \(\sigma^2\) ab:

\[ p(\sigma^2|\alpha,\beta;y)\propto \sigma^{-n}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\alpha-\beta x_i)^2\right)\cdot \sigma^{-2(a+1)}\exp\left(b\sigma^{-2}\right)\\ \propto \sigma^{2(a+n/2-1)}\exp\left(\sigma^{-2}\left(b+\frac{1}{2}\sum_{i=1}^n(y_i-\alpha-\beta x_i)^2\right)\right) \]

Das ist der Kern einer Invers-Gammaverteilung mit Parametern

• \(\tilde{a}=a+n/2\)
• \(\tilde{b}=b+0.5\sum_{i=1}^n(y_i-\alpha-\beta x_i)^2\)

Semi-Konjugiertheit

• Wir erinnern uns, dass wir (angeblich) semi-Konjugierte Prioris benutzt haben.
• Semi-Konjugiertheit hatten wir wie folgt definiert:

Die vollständig bedingte Posteriori eines Parameters hat die selbe Verteilung wie die Priori des Parameters, bis auf andere Parameter.

Die Prioris sind also in der Tat normalverteilt, die full conditionals entsprechen den Priori_Verteilungen ist es nicht verwunderlich, dass wir auf diese full conditionals gekommen sind:

• \(\alpha\to\) Normalverteilung
• \(\beta\to\) Normalverteilung
• \(\sigma^2\to\) Invers-Gamma-Verteilung

Implementierung

Die Implementierung des Gibbs-Samplers in R finden Sie hier

Ergebnisse

Sehen wir uns die Ergebnisse der Ziehungen an. Die folgenden Graphiken zeigen die Ziehungen der drei Parameter über 20 Iteration, die sogenannten Trace Plots:

Burn-In

• Wir sehen, dass die Zufallszahlen in den ersten paar Iterationen einen sehr starken Trend aufweisen.
• Danach sind die Ziehungen annähernd stabil aus einem bestimmten Bereich.
• Diesen Bereich am Anfang nennen wir den burn-in.
• Der burn-in hängt von verschiedenen Eigenschaften des Algorithmus ab, insbesondere aber von den Startwerten.
• Theoretisch können die Startwerte beliebig gewählt werden, für ein große Anzahl von Iterationen erhalten wir Ziehungen aus der richtigen Verteilung. Ein “falscher” Startwert kann aber dazu führen, dass der burn-in länger wird.
• Im Beispiel wurde \(\beta_0=-10\) gesetzt, so dass der Algorithmus einige Iterationen braucht, um in den “richtigen Bereich” zu kommen.
• Die korrekte Bestimmung des burn-in, also ab wann der Algorithmus wirklich Ziehungen aus der gemeinsamen Verteilung liefert, ist nicht trivial.
• Wir werden später sogenannten Kovergenzdiagnostiken besprechen.
• Vorerst bestimmen wir denn burn-in visuell und großzügig.
• Hier legen wir ihn auf 100 Iterationen fest, sprich wir löschen die ersten 100 Iterationen:

alpha.sample<-alpha.sample[-(1:100)]
beta.sample<-beta.sample[-(1:100)]
sigma2.sample<-sigma2.sample[-(1:100)]
AnzahlIterationen<-AnzahlIterationen-100

Trace plots ohne burn-in

Mittelwert und Anzahl der Iterationen

Eine andere Frage ist, wieviele Iterationen wir brauchen. Die Antwort auf diese Frage hängt stark davon ab, welche Schlüsse man ziehen will.

Den Posteriori-Erwartungswert kann man über den Mittelwert der Ziehungen schätzen. Eine Möglichkeit ist, sich den Mittelwert über die bisherige Iterationen im Verlauf der Iterationen zu visualisieren (ohne den burn-in):

plot(cumsum(beta.sample)/(1:AnzahlIterationen), type="l",
     xlab="Iteration", ylab=expression(beta))

Wir sehen, dass der Mittelwert “konvergiert”, also mit steigender Anzahl von Iterationen stabil auf einer Wert zustrebt. Nach 900 Iterationen (ohne den burn-in) ändern sich (bei \(\beta\)) die ersten drei Nachkommastellen des Mittelwerts nicht mehr. Die Anzahl der Iterationen ist hier (für den Posteriori-Erwartungswert) groß genug.

Schätzer

Aus den Ziehungen können wir nun Punkt- und Intervallschätzer berechnen. Konkret schätzen wir z.B. den Posteriori-Erwartungswert über den Mittelwert der Ziehungen und das symmetrische 95%-Kredibilitätsintervall über die 2.5%- und 97.5%-Quantile der Ziehungen:

	Posteriori-Erwartungswert	Kredibilitätsintervall
\(\alpha\)	59.437	(51.129, 68.654)
\(\beta\)	0.208	(0.048, 0.359)
\(\sigma^2\)	146.915	(96.358, 221.647)

Interpretation

Die Parameter lassen sich wie im linearen Regressionsmodell üblich interpretieren:

• \(\beta\) ist der lineare Einfluß der Kovariable “Anteil des Agrarsektors”. Pro Punkt, den der Anteil der Agrarsektors höher ist, ist der Fertilitätsindex im Mittel um 0.208 Punkte höher
• Der Intercept liegt bei 59.437, also bei einem Anteil des Agrarsektors von 0 liegt der Fertilitätsindex im Mittel bei 59.437.
• Das 95%-Kredibilitätsintervall von \(\beta\) umfasst nicht die 0, dass heißt der Anteil des Agrarsektors hat “signifikanten” Einfluß auf den Fertilitätsindex.
• Für diese Frage kann man auch \(P(\beta\leq 0|y)\) durch die relative Häufigkeit von \(\beta\leq 0\) schätzen; diese liegt bei 0.0066667. Dies kann man als Bayesianische Entsprechung des p-Wertes interpretieren.
• Der Punktschätzer der Fehlervarianz liegt bei 146.915, bzw. des Standardfehlers bei 0.4560702. Diese sagt aus, wie stark die Beobachtungen vom Modell abweichen.

Dichteschätzer

Schliesslich können wir uns auch die marginalen Posterioris der drei Parameter ansehen. Dazu berechnen wir einen Kerndichteschätzer aus den Ziehungen. Für die Schätzung der Dichte braucht man sehr viel mehr Ziehungen, daher basieren folgende Graphiken auf 20000 Ziehungen man

plot(density(alpha.sample[-(1:1000)]), xlab=expression(alpha), ylab="", main="Geschätzte Posteriori-Dichte")

plot(density(beta.sample[-(1:1000)]), 
     xlab=expression(beta), ylab="", 
     main="Geschätzte Posteriori-Dichte")

Bei der Posteriori von \(\beta\) sehen wir nochmal, dass die 0 (sprich: kein Einfluß der Kovariable) “am Rand” der Posteriori-Dichte liegt.