Kapitel 4

Bayes-Inferenz

4.1 Uberblick

”Definition” bayesianischer Inferenz: Anpassen eines Wahrscheinlichkeitsmodells an eine Men-

ge von Daten.

Ergebnis: Wahrscheinlichkeitsverteilung fur die Parameter des Modells (und andere unbeob-achtete Großen, zum Beispiel Vorhersagen fur neue Beobachtungen).

Idealisierter Prozess bayesianischer Datenanalyse:

1. Stelle ein volles Wahrscheinlichkeitsmodell oder eine gemeinsame Wahrscheinlichkeits-verteilung fur alle beobachtbaren und unbeobachtbaren Großen auf. Dabei ist Wissenuber das zugrundeliegende wissenschaftliche Problem und den datengenerierenden Pro-zess hilfreich.

2. Berechnung der Posterioriverteilung der unbeobachtbaren Großen (Parameter, missingdata, . . .): bedingte Wahrscheinlichkeitsverteilung der unbeobachtbaren Großen gegebendie beobachteten Daten.

3. Modelldiagnose: Fit, Sensitivitat (bezuglich der Annahmen in 1.).

Ergebnis:”koharentes System”

4.2 Exchangeability

Exchangeability (”Austauschbarkeit”) ist ein wichtiges Konzept fur die statistische Modell-

bildung. Es geht auf de Finetti zuruck.

Definition 4.1 (Finite Exchangeability). Die Zufallsgroßen X1, . . . , Xn sind exchangeablebezuglich des Wahrscheinlichkeitsmaßes P, wenn

P(x1, . . . , xn) = P(xπ(1), . . . , xπ(n))

84

fur alle Permutationenπ : {1, . . . , n} → {1, . . . , n}

gilt. Existiert eine Dichte f zu P, so gilt entsprechend:

f(x1, . . . , xn) = f(xπ(1), . . . , xπ(n)).

Definition 4.2 (Infinite Exchangeability). Die unendliche Folge X1, X2, . . . ist exchangeable,wenn jede endliche Teilfolge exchangeable ist.

Bemerkung. Analog zu obigen Definitionen kann auch bedingte Exchangeability definiertwerden, etwa im Regressionsfall fur Y1|x1, . . . , Yn|xn.

Satz 4.3 (Darstellungssatz fur 0-1 Zufallsvariablen). Sei X1, X2, . . . eine unendliche Fol-ge binarer Zufallsvariablen, die exchangeable sind, mit zugrundeliegendem Wahrscheinlich-keitsmaß P. Dann existiert eine Verteilungsfunktion Q, so dass die gemeinsame Dichtef(x1, . . . , xn) folgende Gestalt hat:

f(x1, . . . , xn) =

∫ 1

0

n∏i=1

θxi(1− θ)1−xi dQ(θ)

mitQ(θ) = lim

n→∞P(yn/n ≤ θ)

und

yn =n∑i=1

xi, θ = limn→∞

yn/n.

Interpretation:

1. Bedingt auf θ sind X1, X2, . . . unabhangige, identisch verteilte Bernoulli-Zufallsgroßen.

2. θ wird eine Verteilung zugeordnet.

3. Q ist der”Glaube”(

”Belief”) uber den Grenzwert der relativen Haufigkeit der Einsen.

Konventionelle Schreibweise:

f(x1, . . . , xn|θ) =n∏i=1

θxi(1− θ)1−xi .

Satz 4.4. Wenn die benotigten Dichten existieren und X1, X2, . . . eine (unendliche) Folgereellwertiger Zufallsgroßen ist, dann gilt

f(x1, . . . , xn) =

∫Θ

n∏i=1

f(xi|θ) dQ(θ).

85

Wir betrachten nun die a posteriori pradiktive Verteilung oder bedingte pradiktive Vertei-lung von zukunftigen (unbeobachteten) Daten xm+1, . . . , xn gegeben die beobachteten Da-ten x1, . . . , xm:

f(xm+1, . . . , xn|x1, . . . , xm) =f(x1, . . . , xn)

f(x1, . . . , xm)(Satz von Bayes)

=

∫Θ

n∏i=1

f(xi|θ) dQ(θ)

∫Θ

m∏i=1

f(xi|θ) dQ(θ)

(Darstellungssatz)

=

∫Θ

m∏i=1

f(xi|θ)n∏

i=m+1f(xi|θ) dQ(θ)

∫Θ

m∏i=1

f(xi|θ) dQ(θ)

=

∫Θ

n∏i=m+1

f(xi|θ) ·

m∏i=1

f(xi|θ) dQ(θ)∫Θ

m∏i=1

f(xi|θ) dQ(θ)

.

Dabei istm∏i=1

f(xi|θ) dQ(θ)∫Θ

m∏i=1

f(xi|θ) dQ(θ)

= dQ(θ|x1, . . . , xm)

die Posterioriverteilung fur θ gegeben Daten x1, . . . , xm. Hier haben wir aus”vergangenen”,

beobachteten Daten fur zukunftige Beobachtungen gelernt. Eine Erweiterung auf andere Zu-fallsgroßen ist moglich:

Satz 4.5 (Allgemeiner Darstellungssatz). Sei X1, X2, . . . eine unendliche Folge reellwertigerZufallsvariablen, die exchangeable sind, mit zugrundeliegendem Wahrscheinlichkeitsmaß P.Dann existiert ein Wahrscheinlichkeitsmaß Q uber F , dem Raum aller VerteilungsfunktionenF auf R, so dass

P(x1, . . . , xn) =

∫F

n∏i=1

F (xi) dQ(F ),

wobeiQ(F ) = lim

n→∞P(Fn),

wobei Fn die zu x1, . . . , xn gehorende empirische Verteilungsfunktion bezeichnet.

Man beachte, dass obige Aussage sich (auch) auf nichtparametrische Inferenz bezieht. Sosteht Q(F ) fur eine Prioriverteilung auf dem Raum aller Verteilungsfunktionen.

86

4.3 Bayes-Inferenz im Schnelldurchlauf

Notation:

• X: beobachtete Daten

• X: unbeobachtete Daten

• θ: Parameter

Ziel:

• Wahrscheinlichkeitsaussagen bedingt auf beobachtete Daten

• Vorhersage / pradiktive Inferenz

Basiskomponenten in der Bayes-Inferenz:

• p(θ) Prioriverteilung

• f(x|θ) Datenverteilung

• f(θ|x) Posterioriverteilung

• f(x|x) pradiktive Verteilung

Nach dem Satz von Bayes ist die gemeinsame Verteilung von (θ, x) gleich

f(θ, x) = f(x|θ) · p(θ),

deshalb

f(θ|x) =f(θ, x)

f(x)=f(x|θ)p(θ)f(x)

,

wobei

f(x) =∑θ∈Θ

f(x|θ)p(θ), falls θ diskret,

f(x) =

∫Θf(x|θ)p(θ) dθ, falls θ stetig.

Unnormalisierte Posteriori:f(θ|x) ∝ f(x|θ)p(θ)

A priori pradiktive Verteilung (vor Beobachtung der Daten):

f(x) =

∫Θf(θ, x) dθ =

∫Θf(x|θ) p(θ) dθ

A posteriori pradiktive Verteilung (nach Beobachtung der Daten x):

f(x|x) =

∫Θf(x, θ|x) dθ =

∫Θf(x|θ, x) f(θ|x) dθ =

∫Θf(x|θ) f(θ|x) dθ,

da x bedingt unabhangig von x gegeben θ ist.

87

Likelihood und Odds Ratios

Die Likelihoodfunktion ist f(x|θ) als Funktion von θ nach Beobachtung von x. Die Datenbeeinflussen die Posteriori-Inferenz also nur uber die Likelihood. Die Posteriori-Odds von θ1

verglichen mit θ2 sind

fθ|x(θ1|x)

fθ|x(θ2|x)=

f(x|θ1)f(θ1)f(x)

f(x|θ2)f(θ2)f(x)

=f(x|θ1)

f(x|θ2)· f(θ1)

f(θ2),

es gilt alsoPosteriori-Odds = Priori-Odds× Likelihoodquotient.

4.4 Wiederholung: Modelle mit einem Parameter

• Gemeint ist: θ ist skalar.

• Prioriverteilung kann mehr als einen Parameter haben.

• Hier funktionieren folgende Konzepte gut:

– Konjugierte Prioriverteilungen, zum Beispiel bei (einparametrischen) Expo-nentialfamilien.

Vorteil: Analytische Berechenbarkeit, keine Simulation notig.

Nachteil: Fur komplexe Modelle meist nicht verfugbar, deshalb eher als Bausteinin komplizierteren Modellen verwendet.

– Referenzprioris/Referenzanalyse:

∗ Idee: Priori so wahlen, dass die Daten auch im Fall geringen Stichproben-umfangs die Posterioriverteilung dominieren (

”let the data speak for themsel-

ves”). Dies benotigt entscheidungs- und informationstheoretische Grundlagen.Suche nach nicht-informativen Prioriverteilungen: Im skalaren Fall zum Bei-spiel

0 < θ < 1 → ψ = log

(θ

1− θ

).

∗ Jeffreys’ Priorip(θ) ∝

√I(θ)

ist invariant gegenuber bijektiven Transformationen von θ.

Beispiel 4.1 (Binomial- und Negative Binomialverteilung).

1. Binomialverteilung: Die Likelihood lautet

f(x|θ) =n∏i=1

θxi(1− θ)1−xi .

Als Referenzpriori kann Jeffreys’ Priori, Beta(

12 ,

12

), verwendet werden:

p(θ) ∝ θ−1/2(1− θ)−1/2.

88

Sei y =∑n

i=1 xi. Dann ist die Referenzposteriori:

f(θ|x) ∝ f(x|θ)p(θ)∝ θy(1− θ)n−yθ−1/2(1− θ)−1/2

= θy−1/2(1− θ)n−y−1/2.

Dies entspricht dem Kern der Dichte einer Beta(

12 + y, 1

2 + n− y)–Verteilung. f(θ|x)

ist auch fur die Extremfalle y = 0 oder y = n noch”

proper”. Verwendet man dagegenHaldane’s Priori

p(θ) ∝ θ−1(1− θ)−1,

die eine uneigentliche Priori ’Beta(0, 0)’ darstellt, ist die Posteriori Beta(y, n− y) furdie Extreme y = 0 oder y = n nicht proper.

2. Negative Binomialverteilung: Sei X die Anzahl an Versuchen bis y ≥ 1 Erfolge eintre-ten. Dann lautet die Likelihood

f(x|θ) ∝(x− 1

y − 1

)θy(1− θ)x−y fur x ≥ y.

Die Referenzpriori ist durch Jeffreys’ Priori fur die geometrische Verteilung gegeben(das entspricht y = 1):

p(θ) ∝ θ−1(1− θ)−1/2,

woraus die Referenzposteriori

f(θ|x) ∝ θy−1(1− θ)x−y−1/2,

also eine Beta(y, x−y+1/2)–Verteilung, resultiert. Da y ≥ 1 und x ≥ y, ist auch diesea posteriori stets proper.

Bemerkung. Konzepte fur eindimensionale Modelle sind im mehrdimensionalen Fall imAllgemeinen schwierig umzusetzen bzw. umstritten (zum Beispiel Verwendung von Referenz-prioris). Man geht daher oft zu sogenannten hierarchischen Modellen uber: Fuge zusatzlicheStufen in das Modell ein mit dem Ziel, die Posteriori-Analyse starker von Priori-Annahmenzu entkoppeln.

4.5 Mehr-Parameter-Modelle

4.5.1 Normalverteilung

Wir betrachten in diesem Abschnitt Daten x1, . . . , xn|µ, σ2 i.i.d.∼ N(µ, σ2) mit µ, σ2 unbe-kannt.

(i) Gemeinsame Posterioriverteilung von µ, σ2|x:

f(µ, σ2|x) ∝ f(x|µ, σ2) · p(µ, σ2)

89

(ii) Bedingte Posterioriverteilungen von µ|σ2, x bzw. σ2|µ, x:

f(µ|σ2, x) bzw. f(σ2|µ, x)

(iii) Marginale Posterioriverteilung von µ|x:

f(µ|x) =

∫f(µ, σ2|x) dσ2 =

∫f(µ|σ2, x)f(σ2|x) dσ2

I. Nichtinformative Prioriverteilung

Ist nur einer der beiden Parameter unbekannt, so wahlt man oft folgende Prioriverteilungen(Jeffreys’ Prioris):

σ2 bekannt: p(µ) ∝ const,

µ bekannt: p(σ2) ∝ (σ2)−1.

Eine Moglichkeit, daraus eine mehrdimensionale Priori zu konstruieren, ist:

p(µ, σ2) = p(µ) · p(σ2) ∝ (σ2)−1,

d.h. wir nehmen unabhangige Prioris fur µ und σ2 an. Die gemeinsame Posteriorivertei-lung f(µ, σ2|x) lautet dann:

f(µ, σ2|x) ∝ Likelihood × Priori

=

[n∏i=1

1√2πσ−1 exp

(− 1

2σ2(xi − µ)2

)]· (σ2)−1

∝ σ−n−2 exp

(− 1

2σ2

n∑i=1

(xi − µ)2

)

= σ−n−2 exp

(− 1

2σ2

((n− 1)s2 + n(x− µ)2

))mit s2 =

∑ni=1(xi − x)2/(n−1). Die bedingte Posteriori von µ, f(µ|σ2, x), kann auf den

Fall mit bekannter Varianz σ2 zuruckgefuhrt werden. Aus Statistik IV ist bekannt, dassf(µ|σ2, x) ∼ N(x, σ2/n). Fur die marginale Posteriori f(σ2|x) hat man

f(σ2|x) =

∫f(µ, σ2|x) dµ

∝∫σ−n−2 exp

(− 1

2σ2

((n− 1)s2 + n(x− µ)2

))dµ

∝ σ−n−2 exp

(− 1

2σ2(n− 1)s2

)∫exp

(− 1

2σ2n(x− µ)2

)dµ.

Es gilt ∫exp

(− 1

2σ2n(x− µ)2

)dµ =

√2π

σ2

n

90

und damit

f(σ2|x) ∝ σ−n−2 exp

(− 1

2σ2(n− 1)s2

)√2π

σ2

n

∝ (σ2)−n−2

2 (σ2)12 exp

(− 1

2σ2(n− 1)s2

)= (σ2)−

n+12 exp

(− 1

2σ2(n− 1)s2

).

Der Kern dieser Dichte gehort zur inversen Gamma-Verteilung mit den Parametern (n−1)/2und (n− 1)s2/2.

Wegenf(µ, σ2|x1, . . . , xn) = f(µ|σ2, x1, . . . , xn) · f(σ2|x1, . . . , xn)

kann die gemeinsame Posterioriverteilung von µ, σ2|x1, . . . , xn nun wie folgt simuliert werden:

Algorithmus 2 : Direkte Simulation der gemeinsamen Posterioriverteilungbei nichtinformativer Priori

Wiederhole fur s = 1, . . . , S:

Schritt 1: Ziehe (σ2)(s) aus IG(n−1

2 , n−12 s2

).

Schritt 2: Ziehe (µ)(s) aus N(x, 1n(σ2)(s)).

Man erhalt Paare[(µ(1), (σ2)(1)), . . . , (µ(S), (σ2)(S))

].

σ2 als Nuisance-Parameter

Interessiert nur µ, so gibt es (mindestens) zwei Moglichkeiten zur Simulation:

1. Simuliere die gemeinsame Posteriori f(µ, σ2|x) gemaß obigem Algorithmus und betrach-te nur die Ziehungen µ(1), . . . , µ(S).

2. Berechne direkt die marginale Posteriori f(µ|x):

f(µ|x) =

∫ ∞0

f(µ, σ2|x) dσ2.

Fuhrt man die Substitution z = A/(2σ2) mit A = (n−1)s2 +n(µ−x)2 durch, so erhaltman wegen

σ2 =1

2Az−1 und dσ2 = −2A−1σ4dz = −1

2Az−2dz

fur f(µ|x) ∫ ∞0

f(µ, σ2|x) dσ2 ∝∫ ∞

0A−

n+22 z

n+22 exp(−z) A z−2dz

=

∫ ∞0

A−n2 z

n−22 exp(−z) dz

= A−n2

∫ ∞0

zn−22 exp(−z) dz.

91

Allgemein gilt fur a > 0 und m > −1:∫ ∞0

xm exp(−ax) dx =Γ(m+ 1)

am+1.

Daraus folgt, dass das Integral konstant bezuglich µ ist und somit

f(µ|x) ∝ A−n2

=[(n− 1)s2 + n(µ− x)2

]−n2

=

[1 +

(µ− x)2

(n− 1)s2/n

]−n2

=

1 +1

n− 1

(µ− x

s√n

)2−n2

was der Kern einer skalierten nichtzentralen t-Verteilung mit Skalenparameterm = s/

√n, Lokationsparameter l = x und ν = n− 1 Freiheitsgraden ist. Allgemein hat

der Kern der Dichte einer solchen allgemeinen t-Verteilung die Gestalt

kern(f(θ)) =

[1 +

1

ν

(θ − lm

)2]− ν+1

2

.

Bemerkung.

1. Statt σ2 lasst sich auch die sogenannte Prazision κ = (σ2)−1 verwenden. Bei Verwen-dung von p(µ, κ) ∝ (κ)−1 folgt, dass κ|x ∼ Gamma

(n−1

2 , n−12 s2

).

2. Statt inverser Gammaverteilung wird haufig der Spezialfall einer sogenannten skalierteninversen χ2-Verteilung inv-χ2 verwendet (siehe unten).

II. Konjugierte Prioriverteilung

Verwende gemaß Bemerkung 2 die skalierte inverse χ2(ν0, σ20)-Verteilung als Priori.

Vorteil: Bessere Interpretation (das werden wir allerdings erst dann verstehen, wenn wirinformative Prioriverteilungen in Form von (inversen) Gammaverteilungen betrachten).

Nachteil: Diese Vorgehensweise ist”non-standard”.

Zufallszahlen aus einer skalierten inversen χ2 -Verteilung kann man wie folgt simulieren:

Algorithmus 3 : Simulation von θ ∼ inv-χ2(ν0, σ20)

1. Ziehe X∗ ∼ χ2(ν0).

2. Setze θ =ν0σ2

0X∗ .

92

Es gilt:

inv-χ2(ν0, σ20) = IG

(ν0

2,ν0

2σ2

0

).

Dies lasst sich mit dem Transformationssatz fur Dichten verifizieren: Definiere α = ν0/2 undβ = 1/2, so dass X∗ ∼ Gamma(α, β). Die Umkehrfunktion der Transformation in Schritt 2lautet

X∗ = g−1(θ) =ν0σ

20

θund die zugehorige Ableitung nach θ

(g−1)′(θ) = −ν0σ20

θ2.

Man erhalt somit:

f(θ) = fX∗(g−1(θ)) · |(g−1)′(θ)| = (βν0σ

20)α

Γ(α)θ−α−1 exp

(−βν0σ

20/θ).

Dies ist die Dichte einer inversen Gammaverteilung mit Parametern (α, βν0σ20), welche gerade

der gewunschten inversen χ2-Verteilung entspricht. Eine mogliche Parametrisierung ist nun

p(µ, σ2) = p(µ|σ2) · p(σ2)

mit

µ|σ2 ∼ N

(µ0,

σ2

κ0

)und σ2 ∼ inv-χ2(ν0, σ

20)

Man schreibt hierfur kurz: N-inv-χ2(µ0,

σ20κ0

; ν0, σ20

). Die Prioris sind nunmehr voneinander

abhangig.

Sei nun also a priori (µ, σ2) ∼ N-inv-χ2(µ0,

σ20κ0

; ν0, σ20

). Die Prioridichte lautet dann

p(µ, σ2) =1√

2πσ2κ−10

exp

(− 1

2σ2

κ0

(µ− µ0)2

)×

(12ν0σ

20)ν0/2

Γ(ν0/2)(σ2)−( ν02 +1) exp

(−1

2

ν0σ20

σ2

)

∝ (σ2)−12 (σ2)−( ν02 +1) exp

(− 1

2σ2

[ν0σ

20 + κ0(µ− µ0)2

]).

Die gemeinsame Posteriori bei gegebenen Daten x = (x1, . . . , xn) aus N(µ, σ2) ergibt sichzu:

f(µ, σ2|x) ∝ (σ2)−12 (σ2)−( ν02 +1) exp

(− 1

2σ2

[ν0σ

20 + κ0(µ− µ0)2

])×(σ2)−

n2 exp

(− 1

2σ2

((n− 1)s2 + n(x− µ)2

)).

Man kann zeigen (vgl. Ubung), dass die Posteriori wieder N-inv-χ2-verteilt ist mit Parametern

µn =

(κ0

κ0 + n

)µ0 +

(n

κ0 + n

)x,

κn = κ0 + n,

νn = ν0 + n,

νnσ2n = ν0σ

20 + (n− 1)s2 +

κ0n

κ0 + n(x− µ0)2.

Die Interpretation der Parameter ist wie folgt:

93

• µn ist gewichteter Mittelwert aus Stichprobenmittel und Priori-Erwartungswert. In denGrenzfallen κ0 →∞ ist µn = µ0 bzw. fur n→∞ ist µn = x.

• νn sind die Posteriori-Freiheitsgrade als Summe von Priori-Freiheitsgraden und Stich-probenumfang.

• Die Posteriori-Quadratsumme νnσ2n lasst sich partitionieren in die Priori-Quadratsum-

me ν0σ20, die Quadratsumme (n − 1)s2 der Stichprobe und einen Term, der die Unsi-

cherheit, die durch die Differenz zwischen Stichprobenmittel und Priori-Erwartungswertensteht, quantifiziert.

Die bedingte Posteriori von µ|σ2, x ist

µ|σ2, x ∼ N

(µn,

σ2

κn

)ˆ= N

(κ0

κ0 + nµ0 +

n

κ0 + nx,

σ2

κ0 + n

)ˆ= N

( κ0σ2µ0 + n

σ2xκ0σ2 + n

σ2

,1

κ0σ2 + n

σ2

).

Die Gewichte κ0/σ2 und n/σ2 entsprechen der Priori- bzw. Datenprazision. Die marginale

Posteriori von σ2|x istσ2|x ∼ inv-χ2(νn, σ

2n).

Dies ermoglicht die Simulation der gemeinsamen Posteriori Verteilung:

Algorithmus 4 : Direkte Simulation der gemeinsamen Posterioriverteilungbei konjugierter Priori

Schritt 1: Ziehe (σ2)∗ aus inv-χ2(νn, σ2n).

Schritt 2: Ziehe µ∗ aus N(µn,

(σ2)∗

κn

).

Die marginale Posteriori von µ|x lautet

f(µ|x) ∝[1 +

κn(µ− µn)2

νnσ2n

]−(νn+1)/2

.

Dies entspricht einer t-Verteilung mit νn Freiheitsgraden, Lokationsparameter µn und Ska-lenparameter σ2

n/κn.

III. Semi-konjugierte Prioriverteilung

Die Parameter µ und σ2 sollen nun a priori unabhangig sein. Wir wahlen deshalb a priori

µ|σ2 ∼ N(µ0, τ20 ) und σ2 ∼ inv-χ2(ν0, σ

20).

94

Der einzige Unterschied zum Fall der konjugierten Priori ist also, dass wir τ20 statt σ2

0/κ0

verwenden und so die Prioris entkoppeln. Es folgt:

p(µ, σ2) = p(µ) · p(σ2).

Die resultierende gemeinsame Posteriori hat allerdings keine Form, die einer bekannten Ver-teilung zugeordnet werden kann. Allerdings ist f(µ|σ2, x) explizit berechenbar und f(σ2|x)einfach zu simulieren. Die bedingte Posteriori ist µ|σ2, x ∼ N(µn, τ

2n) mit

µn =

1τ20µ0 + n

σ2x

1τ20

+ nσ2

und τ2n =

11τ20

+ nσ2

.

Zur Herleitung der Posteriori f(σ2|x) benutzt man, dass

f(σ2|x) =f(µ, σ2|x)

f(µ|σ2, x)

und (salopp), dass

f(µ, σ2|x) ∝ N(µ|µ0, τ20 )× inv-χ2(σ2|ν0, σ

20)×

n∏i=1

N(xi|µ, σ2).

Die marginale Posteriori hat dann die Struktur

f(σ2|x) ∝N(µ|µ0, τ

20 ) · inv-χ2(σ2|ν0, σ

20) ·∏ni=1N(xi|µ, σ2)

N(µ|µn, τ2n)

.

Da f(σ2|x) nicht von µ abhangen kann, konnen die entsprechenden Terme ignoriert werden,die nur von µ abhangen. Man beachte jedoch, dass der Nenner uber die Parameter µn, τ

2n

noch von σ2 abhangt. Setzen wir µ = µn, so erhalten wir

f(σ2|x) ∝ τn exp

(− 1

2τ20

(µn − µ0)2

)· inv-χ2(σ2|ν0, σ

20) ·

n∏i=1

1√2πσ2

exp

(− 1

2σ2(xi − µn)2

).

Diese Verteilung lasst sich beispielsweise uber einen empirischen CDF-Sampler simulieren.Voraussetzung hierfur ist, dass wir die Dichte einer (univariaten) Zufallsvariable bis auf eineProportionalitatskonstante c kennen.

Algorithmus 5 : Empirischer CDF-Sampler

• Diskretisiere den Trager der zu simulierenden Verteilung in eine Menge von NPunkten x1 ≤ . . . ≤ xN .

• Evaluiere die bis auf Proportionalitat bekannte Dichte an x1 ≤ . . . ≤ xN , um Wertef1, . . . , fN zu erhalten.

• Schatze die Proportionalitatskonstante c uber c = f1 + . . .+ fN .

• Ziehe Zufallszahlen aus x1 ≤ . . . ≤ xN gemaß den Wahrscheinlichkeitenf1/c, . . . , fN/c.

95

Dies fuhrt zu folgendem Algorithmus zur Simulation aus f(µ, σ2|x):

1. Ziehe (σ2)∗ aus der marginalen (approximativen) Posteriori gemaß CDF-Sampler.

2. Ziehe µ aus f(µ|(σ2)∗, x).

Abschließend betrachten wir noch die pradiktive Posterioriverteilung fur zukunftige Be-obachtungen x gegeben Daten x. Diese lautet

f(x|x) =

∫ ∫f(x|µ, σ2, x)f(µ, σ2|x) dµ dσ2 =

∫ ∫f(x|µ, σ2)f(µ, σ2|x) dµ dσ2.

Simulation:

1. Ziehe (µ∗, (σ2)∗) aus der Posteriori wie oben beschrieben.

2. Ziehe x ∼ N(µ∗, (σ2)∗).

4.5.2 Dirichlet-Multinomial Modell

Im Dirichlet-Multinomial-Modell wird fur die Daten y1, . . . , yn eine Multinomialverteilungangenommen, also die Verallgemeinerung der Binomialverteilung auf mehr als zwei mogli-che Ereignisse bei festem Stichprobenumfang n. Beispielsweise konnte eine fest vorgegebeneAnzahl an Personen nach ihrer Parteipraferenz befragt werden.

Eine multinomialverteilte Zufallsvariable Y kann k mogliche Auspragungen annehmen (zumBeispiel CDU/CSU, SPD, FDP, Grune, Linke, andere). Die Zufallsvariable Xj , 1 ≤ j ≤ k,

bezeichnet die Anzahl der j-ten Auspragung in der Stichprobe; es gilt∑k

j=1Xj = n. Der

Parameter θj = P(Y = j) ∈ [0, 1] fur Y ∈ {1, . . . , k} mit∑k

j=1 θj = 1 bezeichnet dieWahrscheinlichkeit fur die Auspragung j.

Die Likelihood von θ = (θ1, . . . , θk)> bei Beobachtungen x = (x1, . . . , xk)

> lautet

L(θ) = f(x|θ) ∝k∏j=1

θxjj .

Wegen der Restriktion∑k

j=1 θj = 1 liegen faktisch nur k − 1 Parameter vor, denn der k-teParameter lasst sich deterministisch durch θk = 1 − θ1 − . . . θk−1 berechnen. Die Likelihoodlasst sich daher auch in der Form

L(θ) ∝

k−1∏j=1

θxjj

(1− θ1 − . . .− θk−1)xk

schreiben.

Die zur Multinomialverteilung konjugierte Verteilung ist die sogenannte Dirichletverteilung,eine Verallgemeinerung der Beta-Verteilung, geschrieben

θ = (θ1, . . . , θk)> ∼ Dirichlet(α1, . . . , αk) = D(α),

96

mit Dichtefunktion

p(θ) =Γ(α1 + . . .+ αk)

Γ(α1) · . . . · Γ(αk)θα1−1

1 · . . . · θαk−1k ,

wobei αj > 0 fur alle j = 1, . . . , k und wieder θj ∈ [0, 1] mit∑k

j=1 θj = 1. Die Dirichlet-Verteilung spezifiziert also eine Verteilung auf einem (k − 1)-dimensionalen offenen Simplex.

Eigenschaften der Dirichletverteilung:

Definiere α0 =∑k

j=1 αj .

• Momente:

E(θj) =αjα0

, Var(θj) =αj(α0 − αj)α2

0(α0 + 1), Cov(θi, θj) = − αiαj

α20(α0 + 1)

,

wobei die Restriktion∑k

j=1 θj = 1 die negative Korrelation impliziert.

• Modus:

Modus(θ)j =αj − 1

α0 − kist die j-te Komponente des k-dimensionalen Modus.

• Jede Randverteilung ist wieder eine Dirichletverteilung, zum Beispiel

(θi, θj , 1− θi − θj) ∼ Dirichlet(αi, αj , α0 − αi − αj).

Insbesondere istθj ∼ Beta(αj , α0 − αj).

• Die bedingten Verteilungen sind ebenfalls Dirichletverteilt. Setzt man

θ′i =θm

1−∑m−1

r=1 θr, m ≤ i ≤ k,

gegeben die Realisationen θ1, . . . , θm−1, so ist

(θ′m, . . . , θ′k)> ∼ Dirichlet(αm, . . . , αk).

97

Algorithmus 6 : Simulation aus der Dirichletverteilung

• Simulation 1:

1. Ziehe Z1, Z2, . . . , Zk aus (unabhangigen) Gammaverteilungen mitParametern (α1, 1), . . . , (αk, 1).

2. Setze

θj =Zj∑ki=1 Zi

.

• Simulation 2 (”Stick Breaking Prior”):

1. Ziehe θ1 ∼ Beta(α1, α0 − α1).

2. Fur j = 2, . . . , k − 1:

(i) Ziehe Zj ∼ Beta(αj ,∑k

i=j+1 αi).

(ii) Setze θj =(

1−∑j−1

i=1 θi

)Zj .

3. Setze θk = 1−∑k−1

i=1 θi.

Fur x|θ ∼ Multinomial(n; θ1, . . . , θk) und θ|α ∼ Dirichlet(α1, . . . , αk) lautet die Posteriori-verteilung von θ:

f(θ|x) ∝ L(θ) · p(θ|α)

∝k∏j=1

θxjj ·

k∏j=1

θαj−1j

=k∏j=1

θxj+αj−1j ,

d.h. die Posteriori ist Dirichlet(x1 + α1, . . . , xk + αk)-verteilt.

Interpretation der Posteriori

Der Posteriori-Erwartungswert

E[θj |x] =xj + αj∑k

i=1 xi +∑k

i=1 αi=xj + αjn+ α0

lasst sich umschreiben zu

E[θj |x] =α0

α0 + n· αj

α0︸︷︷︸Priori-

Erwartungswert

+n

α0 + n· xj

n︸︷︷︸MLE

.

Der Parameter α0 lasst sich als”a priori Anzahl Beobachtungen” und αj als

”a priori Erfolge”

fur Kategorie j interpretieren.

98

Bemerkung.

1. Die Wahl α1 = . . . = αk = 0 entspricht einer Gleichverteilung fur {log θj}kj=1. Indiesem Fall ist die Posteriori nur dann proper, wenn xj ≥ 1, j = 1, . . . , p.

2. Die Wahl α1 = . . . = αk = 1/2 entspricht Jeffreys’ Priori.

3. Die Wahl α1 = . . . = αk = 1 entspricht einer Priori-Gleichverteilung auf dem Simplex.

Bemerkung. Die Dirichlet-Verteilung eignet sich auch als Priori bei der Analyse von Kon-tigenztafeln mit multinomialem Erhebungsschema:

x1 x2

x3 x4

(n = x1+x2+x3+x4). Erweiterungen auf restringierte Multinomialverteilungen sind moglich(loglineare Modelle).

Ad-hoc Prozedur: Addiere 1/2 zu jedem Eintrag der Kontingenztabelle und berechne dann denMaximum-Likelihood-Schatzer; das entspricht α1 = . . . = αk = 3/2 und der Posteriori-ModusSchatzung

Modus(θ|x)j =xj + αj − 1∑k

i=1 xi +∑k

i=1 αi − k

=xj + 1

2∑ki=1 xi + 1

2k.

4.5.3 Multivariate Normalverteilung

Notation:

• X = (X1, . . . , Xp)> ist p-dimensionaler Zufallsvektor.

• µ = (µ1, . . . , µp) ist p-dimensionaler Erwartungswertvektor.

• Die symmetrische und positiv definite (Notation: Σ > 0) Matrix

Σ =

σ11 σ12 . . . σ1p

σ21 σ22 . . . σ2p...

.... . .

...σp1 σp2 . . . σpp

ist p× p-dimensionale Kovarianzmatrix.

• Eine Beobachtung x = (x1, . . . , xp)> ist MVN(µ,Σ) (multivariat normalverteilt), wenn

f(x|µ,Σ) ∝ |Σ|−1/2 exp

(−1

2(x− µ)>Σ−1(x− µ)

).

99

• Die Likelihood fur n i.i.d. Realisationen x1, . . . ,xn lautet

L(µ,Σ) = f(x1, . . . ,xn|µ,Σ) ∝ |Σ|−n/2 exp

(−1

2

n∑i=1

(xi − µ)>Σ−1(xi − µ)

)

= |Σ|−n/2 exp

(−1

2tr(Σ−1S0)

)mit S0 =

∑ni=1(xi−µ)(xi−µ)>. Die zweite Identitat ergibt sich uber die Umformungen

n∑i=1

(xi − µ)>Σ−1(xi − µ) = tr

(n∑i=1

(xi − µ)>Σ−1(xi − µ)

)

= tr

(n∑i=1

Σ−1(xi − µ)(xi − µ)>

)= tr(Σ−1S0).

I. Konjugierte Prioriverteilung bei unbekanntem µ und bekanntem Σ

Konjugierte Prioriverteilung fur µ bei bekanntem Σ ist

µ ∼ MVN(µ0,Λ0)

mit Λ0 > 0. Die Posteriori fur µ ist

f(µ|x,Σ) ∝ exp

(−1

2(µ− µ0)>Λ−1

0 (µ− µ0)− 1

2

n∑i=1

(xi − µ)>Σ−1(xi − µ)

).

Der Term im Exponenten ist eine quadratische Form in µ, die sich uber eine quadratischeErganzung und Vernachlassigung von Konstanten ergibt. Man erhalt

f(µ|x,Σ) ∼ MVN(µn,Λn)

mit

µn =(Λ−1

0 + nΣ−1)−1 (

Λ−10 µ0 + nΣ−1x

),

Λ−1n = Λ−1

0 + nΣ−1

und x = (∑n

i=1 xi)/n in Analogie zum univariaten Fall.

Die bedingte und marginale Posteriori fur Subvektoren von µ folgen aus den Eigenschaf-ten der multivariaten Normalverteilung: Betrachte die Partitionierungen

µ =

(µ(1)

µ(2)

), µn =

(µ

(1)n

µ(2)n

), Λn =

(Λ

(11)n Λ

(12)n

Λ(21)n Λ

(22)n

).

Dann gilt fur die bedingten Verteilungen

µ(1)|µ(2),x ∼ MVN(µ(1)n + β1|2

(µ(2) − µ(2)

n

),Λ1|2

n

)100

mit

β1|2 = Λ(12)n

(Λ(22)n

)−1,

Λ1|2 = Λ(11)n −Λ(12)

n

(Λ(22)n

)−1Λ(21)n

und fur die marginalen Verteilungen

µ(1) ∼ MVN(µ(1)n ,Λ(11)

n

).

Die pradiktive Posterioriverteilung ist (informell)

f(x,µ|x) = MVN(x|µ,Σ) · MVN(µ|µn,Λn).

Im Exponenten erhalt man eine quadratische Form in (x,µ). (x,µ) sind gemeinsam multiva-riat normalverteilt, und daher folgt x als Randverteilung der beiden Komponenten ebenfallseiner multivariaten Normalverteilung. Die erforderlichen Kenngroßen lassen sich uber dieiterierte Erwartung und Varianz berechnen:

E[x|x] = E[E[x|µ,x]|x] = E[µ|x] = µn

und

Var(x|x) = E[Var(x|µ,x)|x] + Var[E(x|µ,x)|x]

= E[Σ|x] + Var[µ|x]

= Σ + Λn.

II. Konjugierte Prioriverteilung bei unbekanntem µ und unbekanntem Σ

In Abschnitt 4.5.1-II (Seite 92) hatten wir als konjugierte Prioriverteilungen fur die Parame-ter µ und σ2 der univariaten Normalverteilung

µ|σ2 ∼ N

(µ0,

σ2

κ0

)und σ2 ∼ inv-χ2(ν0, σ

20),

kurz

µ, σ2 ∼ N-inv-χ2

(µ0,

σ20

κ0; ν0, σ

20

),

verwendet. Hier nun verwenden wir die multivariaten Analoga

µ|Σ ∼ MVN

(µ0,

1

κ0Σ

)und Σ ∼ inv-Wishartν0(Λ−1

0 ),

kurz

µ,Σ ∼ MVN-inv-Wishart

(µ0,

1

κ0Λ0; ν0,Λ0

).

101

Die gemeinsame Prioridichte ist dann

p(µ,Σ) ∝ |Σ|−(ν0+p

2+1)· exp

(−1

2

(tr(Λ0Σ

−1)− κ0(µ− µ0)>Σ−1(µ− µ0)))

.

Dabei bezeichnet ν0 die a priori Anzahl der Freiheitsgrade, κ0 die a priori Anzahl an Mes-sungen auf der Σ-Skala. Die gemeinsame Posterioriverteilung von µ und Σ lautet

µ,Σ|x ∼ MVN-inv-Wishart

(µn,

1

κnΛn; νn,Λn

).

mit

µn =κ0

κ0 + nµ0 +

n

κ0 + nx,

κn = κ0 + n,

νn = ν0 + n,

Λn = Λ0 + S +κ0n

κ0 + n(x− µ0)(x− µ0)>,

wobei S =∑n

i=1(xi − x)(xi − x)> in Analogie zum univariaten Fall. Die Interpretationder Parameter von Seite 93 lasst sich direkt ubertragen: Der Posteriori-Erwartungswert istein gewichtetes Mittel aus Stichprobenmittelwertvektor und Priori-Erwartungswert. Die Ge-samtstreuungsmatrix Λn lasst sich in Priori-Streuungsmatrix, empirische Streuungsmatrixund Streuung zwischen Priori-Erwartungswert und Stichprobenmittel partitionieren.

Die marginale Posteriori fur µ folgt einer multivariaten t-Verteilung mit Parametern µnund Λn/(κn · (νn − p+ 1)), die marginale Posteriori fur Σ einer inversen Wishart-Verteilungmit Parametern νn und Λ−1

n . Zur Simulation aus der gemeinsamen Posteriori oder aus derpradiktiven Verteilung ist folgender Algorithmus anwendbar:

Algorithmus 7 : Simulation aus der gemeinsamen Posteriori und der pradik-tiven Verteilung bei konjugierter Priori

Fur s = 1, . . . , S:

1. Ziehe Σ(s)|x ∼ inv-Wishartνn(Λ−1n

).

2. Ziehe µ(s)|Σ(s),x ∼ MVN(µn,

1κn

Σ(s))

.

3. Ziehe x(s)|µ(s),Σ(s),x ∼ MVN(µ(s),Σ(s)

).

Dann ist (µ(s),Σ(s)) eine Ziehung aus der gemeinsamen Posterioridichte, x eineZiehung aus der pradiktiven Verteilung.

III. Nichtinformative Prioriverteilung bei unbekanntem µ und unbekanntem Σ

Als nichtinformative Prioriverteilung bei unbekanntem µ und Σ eignet sich die multivariateJeffreys’ Priori

p(µ,Σ) ∝ |Σ|−(p+1)/2.

102

Diese entspricht dem Grenzfall κ0 → 0, ν0 → 1, |Λ0| → 0 bei der konjugierten Priori. Fur diePosteriori-Kenngroßen ergibt sich in diesem Fall fur die bedingte Verteilung von µ

µ|Σ,x ∼ MVN

(x,

1

nΣ

),

und fur die marginalen Verteilungen

Σ|x ∼ inv-Wishartn−1(S),

µ|x ∼ mv-tn−p

(x,

1

n(n− p)S

).

Als Beispiel wird im Folgenden die bivariate Normalverteilung betrachtet. Es dient auchdazu, die folgende fur die Bayes-Inferenz wichtige Simulationsstrategie, das sogenannte Gibbs-Sampling, zu illustrieren.

Algorithmus 8 : Gibbs-Sampler

Gegeben: Ein mehrdimensionaler stetiger Zufallsvektor X mit Verteilung π. DerEinfachheit halber seien alle Komponenten von X stetig.

Wir erzeugen im Folgenden eine Markovkette X(0),X(1), . . . mit Startwert X(0) undstationarer Verteilung π. Sei X(t) der aktuelle Zustand der Markovkette. X(t) lasse

sich in k Subvektoren X(t) = (X(t)•1 ,X

(t)•2 , . . . ,X

(t)•k ) partitionieren. Definiere

X(t)−s =

(X

(t)•1 ,X

(t)•2 , . . . ,X

(t)•(s−1),X

(t−1)•(s+1), . . . ,X

(t−1)•k

)fur s = 1, . . . , k. Ferner seien die vollstandig bedingten Verteilungen (

”full

conditionals”)

π(t)•s|−s = π

(X

(t)•s∣∣X(t)−s

)=

π(X

(t)•s ,X

(t)−s

)∫π(X

(t)•s ,X

(t)−s

)dX

(t)•s

gegeben und simulierbar.Dann wird der nachste Zustand X(t+1) komponentenweise wie folgt erzeugt:

Schritt 1: Ziehe X(t+1)•1 ∼ π

(t+1)•1|−1.

Schritt 2: Ziehe X(t+1)•2 ∼ π

(t+1)•2|−2.

...

Schritt k: Ziehe X(t+1)•k ∼ π

(t+1)•k|−k.

Wiederhole diese Schritte ausreichend oft.

Nach einer gewissen Zahl von Wiederholungen kann X(t) als Ziehung aus π angesehen werden.Im Gegensatz zu obigen

”direkten” Simulationsalgorithmen liegen nun allerdings abhangige

Realisationen vor.

103

Beispiel 4.2. (Bivariate Normalverteilung) Sei x bivariat normalverteilt mit Erwartungs-wertvektor (µ1, µ2)> und Kovarianzmatrix(

σ21 σ12

σ12 σ22

)=

(σ2

1 ρρ σ2

2

)=

(1 ρρ 1

)und ρ bekannt.

Bei einer nichtinformativen Priori p(µ1, µ2) ∝ const fur µ1, µ2 reduziert sich die Posterioriauf die Likelihood bei gegebenen Daten x = ((x11, x12)>, . . . , (xn1, xn2)>):

L(µ1, µ2) =

(1

2π

)n (1− ρ2

)−n2 exp

(− 1

2(1− ρ2)A

),

wobei

A =n∑i=1

[(xi1 − µ1)2 − 2ρ(xi1 − µ1)(xi2 − µ2) + (xi2 − µ2)2

]Wir mochten nun die vollstandig bedingten Verteilungen µ1|µ2,x und µ2|µ1,x berechnen.Naturlich ist es aus Symmetriegrunden ausreichend, nur µ1|µ2,x zu ermitteln. Wegen

f(µ1|µ2, x) =f(µ1, µ2|x)

f(µ2|x)=f(µ1, µ2,x)

f(µ2,x)∝ f(µ1, µ2|x)

genugt es, aus der gemeinsamen Posteriori lediglich die Terme zu betrachten, die von derjeweiligen Variablen in der bedingten Verteilung abhangen. Man erhalt dann

f(µ1|µ2,x) ∝ exp

(− 1

2(1− ρ2)n[µ2

1 − 2µ1(x1 + ρ(µ2 − x2))])

mit xj = (∑n

i=1 xij)/n fur j = 1, 2. Eine quadratische Erganzung des Terms in eckigenKlammern um x1 + ρ(µ2 − x2) liefert schließlich das Endresultat

p(µ1|µ2,x) ∝ exp

(− 1

21−ρ2n

(µ1 − [x1 + ρ(µ2 − x2)]

)2).

Dies entspricht dem Kern einer N(x1 + ρ(µ2 − x2), (1 − ρ2)/n)-Verteilung. Der zugehorigeGibbs-Sampler hat die Gestalt:

1. Wahle einen Startwert µ(0)2 .

2. Fur s = 1, . . . , S:

(a) Ziehe µ(s)1 |µ

(s−1)2 ∼ N

(x1 + ρ

(µ

(s−1)2 − x2

), 1−ρ2

n

).

(b) Ziehe µ(s)2 |µ

(s)1 ∼ N

(x2 + ρ

(µ

(s)1 − x1

), 1−ρ2

n

).

104

4.6 Bayesianisches lineares Modell

Modell:y = Xβ + ε ,

wobei y ∈ Rn,X ∈ Rn×p,β ∈ Rp, ε ∈ Rn

Annahmen und Notation:

p = rang(X)

ε = (ε1, . . . , εn)>, εii.i.d∼ N(0, σ2)

Bayesianisch:y|β, σ2,X ∼ MVN(Xβ, σ2I)

Likelihood:

f(y|X,β, σ2) ∝ |σ2I|−1/2 exp

(−1

2(y −Xβ)>(σ2I)−1(y −Xβ)

)= (σ2)−n/2 exp

(− 1

2σ2(y −Xβ)>(y −Xβ)

)

4.6.1 Nichtinformative Prioriverteilung

Die nichtinformative Priorip(β, σ2) ∝ (σ2)−1

ist insbesondere im Fall p� n nutzlich. Fur die gemeinsame Posteriori folgt:

p(β, σ2|y,X) ∝ (σ2)−(n2

+1) exp

(− 1

2σ2(y −Xβ)>(y −Xβ)

).

Sei

β = (X>X)−1X>y,

y = Xβ = X(X>X)−1X>︸ ︷︷ ︸H

y = Hy,

ε = (I −H)y = y − y.

Aus der Theorie linearer Modelle ist bekannt, dass

X>ε = 0, y>ε = 0.

Daraus ergeben sich folgende Umformungen:

(y −Xβ)>(y −Xβ) = [(y − y) + (y −Xβ)]>[(y − y) + (y −Xβ)]

= ε>ε+ (y −Xβ)>(y −Xβ) + 2(y −Xβ)>ε

= ε>ε+ (y −Xβ)>(y −Xβ)

y=Xβ= ε>ε+ (β − β)>X>X(β − β),

105

so dass sich die Posteriori schreiben lasst als

p(β, σ2|y,X) ∝ (σ2)−(n2

+1) exp

(− 1

2σ2

(ε>ε+ (β − β)>X>X(β − β)

)).

Die bedingte Posteriori von β|σ2,y,X ist

p(β|σ2,y,X) ∝ exp

(− 1

2σ2(β − β)>X>X(β − β)

),

da ε = y − Xβ nicht von β abhangt. Man identifiziert obigen Ausdruck als Kern einermultivariaten Normalverteilung, genauer ist

p(β|σ2,y,X) ∼ MVN(β, σ2(X>X)−1),

ein bekanntes Resultat aus der Theorie linearer Modelle.

Die marginale Posteriori von σ2 erhalt man uber Herausintegrieren von β bzw. einfacheruber den Satz von Bayes

f(σ2|y,X) =f(β, σ2|y,X)

f(β|σ2,y,X).

Die Normalisierungskonstante fur die bedingte Posteriori von β ist σ−p/2, also

f(σ2|y,X) ∝ (σ2)−(n2

+1)σp2 exp

(− 1

2σ2ε>ε

)= (σ2)−(n−p

2+1) exp

(− 1

2σ2ε>ε

).

Dies ist der Kern einer inv-χ2(n− p, ε>εn−p

)bzw. IG

(n−p

2 , ε>ε2

)-Verteilung. Es gilt:

E[σ2|y,X] =n− p

n− p− 2· ε>ε

n− p=

ε>ε

n− p− 2.

Algorithmus 9 : Direkte Simulation von β und σ2 im bayesianischen linearenModell

Fur t = 1, . . . , T :

1. Ziehe(σ2)(t)

aus f(σ2|y,X

), d.h. aus inv-χ2

(n− p, ε>εn−p

).

2. Ziehe β(t) aus f(β|(σ2)(t)

,y,X)

, d.h. aus MVN(β,(σ2)(t)

(X>X)−1)

, wobei

β = (X>X)−1X>y.

Eine Alternative zur direkten Simulation besteht in der Verwendung von Gibbs-Sampling,indem zusatzlich zur vollstandig bedingten Dichte von β die vollstandig bedingte Dichtevon σ2,

f(σ2|β,y,X) ∝ p(β, σ2|y,X) ∝ (σ2)−(n2

+1) exp

(− 1

2σ2

(ε>ε+ (β − β)>X>X(β − β)

)),

106

zur Simulation verwendet wird. Dies ist fur festes β der Kern einer skalierten

inv-χ2(n, ε

>ε+(β−β)>X>X(β−β)n

)-Verteilung. Damit lasst sich auch die marginale Poste-

riori von β herleiten:

f(β|y,X) =f(β, σ2|y,X)

f(σ2|y,X,β)=f(β|σ2,y,X) · f(σ2|y,X)

f(σ2|β,y,X)

∝exp

(− 1

2σ2 (β − β)>X>X(β − β))

[ε>ε+(β−β)>X>X(β−β)

n

]n/2exp

(− 1

2σ2 [ε>ε+ (β − β)>X>X(β − β)]) .

Damit:

f(β|y,X) ∝[ε>ε+ (β − β)>X>X(β − β)

]−n/2.

Setzt man

σ2ε =

ε>ε

n− p⇔ ε>ε = (n− p)σ2

ε ,

so ist

f(β|y,X) ∝[(n− p)σ2

ε + (β − β)>X>X(β − β)]−n/2

=

((n− p)σ2

ε ·

[1 +

(β − β)>X>X(β − β)

(n− p)σ2ε

])−n2

∝

[1 +

(β − β)>X>X(β − β)

(n− p)σ2ε

]− (n−p)+p2

.

Dies entspricht dem Kern einer multivariaten t-Verteilung mit n− p Freiheitsgraden, Lokati-onsparameter β und Skalenparameter σ2

ε(X>X)−1, also

β|y,X ∼ mv-tn−p

(β, σ2

ε(X>X)−1

).

Abschließend betrachten wir noch die pradiktive Verteilung fur y|y,X, X. Seien

• m die Anzahl neuer Beobachtungen,

• X neue Beobachtungen von Regressoren der Dimension m× p,

• y der Vektor der Prognosen der Dimension m× 1.

Zur Simulation konnen wir Algorithmus 9 wie folgt erweitern:

107

Algorithmus 10 : Direkte Simulation der pradiktiven Verteilung im bayesia-nischen linearen Modell

Fur t = 1, . . . , T :

1. Ziehe(σ2)(t)

aus f(σ2|y,X

), d.h. aus inv-χ2

(n− p, ε>εn−p

).

2. Ziehe β(t) aus f(β|(σ2)(t)

,y,X)

, d.h. aus MVN(β,(σ2)(t)

(X>X)−1)

, wobei

β = (X>X)−1X>y.

3. Fur i = 1, . . . ,m:

Ziehe yi ∼ MVN(xi>β(t), (σ2)(t)

), wobei xi

> die i-te Zeile von X bezeichnet.

Es ist sogar eine analytische Berechnung moglich:

f(y|y,X, X) ∼ mv-tn−p

[Xβ, σ2

ε

(X(X>X)−1X> + I

)]in Analogie zur Berechnung von Prognose und Prognoseintervallen fur lineare Modelle ausfrequentistischer Sicht.

4.6.2 Konjugierte Prioriverteilung

Im Falle der konjugierten Priori

σ2 ∼ inv-χ2(κ0, σ20) , β|σ2 ∼ MVN(β0, σ

2Σ0)

bzw.β, σ2 ∼ MVN-inv-χ2(β0, σ

20Σ0;κ0, σ

20)

ergibt sich die gemeinsame Posteriori

σ2|y,X ∼ inv-χ2(κn, σ2n) , β|σ2,y,X ∼ MVN(βn, σ

2Σn)

bzw.β, σ2|y,X ∼ MVN-inv-χ2(βn, σ

2nΣn;κn, σ

2n) ,

wobei

βn = (Σ−10 +X>X)−1(Σ−1

0 β0 +X>y) ,

Σn = (Σ−10 +X>X)−1 ,

κn = κ0 + n ,

σ2n = (β>0 Σ−1

0 β0 − β>n Σ−1n βn + y>y + κ0σ

20)/(κ0 + n) .

Als bedingte Posteriori von β ergibt sich

β |σ2,y,X ∼ MVN(βn, σ2Σn),

als marginale Posteriori von σ2

σ2|y,X ∼ inv-χ2(κn, σ2n).

108

4.6.3 Spezialfalle und Erweiterungen

1. Ridge-Regression

Hinweis: Es ist im Allgemeinen sinnvoll, die”echten” Kovariablen (ohne Intercept)

zu standardisieren, um Unterschiede in der Skala zu beseitigen. Ferner geht man zumzentrierten Response uber, so dass der Intercept entfallt. Man erhalt X∗,y∗. Betrachtenun

y∗ = X∗β∗ + ε, ε ∼ N(0, σ2I).

Der Ridge-Schatzer ist durch

[(X∗)>X∗ + λI]−1(X∗)>y∗

mit λ > 0 gegeben. Dieser lasst sich wie folgt bayesianisch interpretieren: Sei

p(β∗) ∼ N(0, τ2I),

d.h. die Komponenten von β∗ sind a priori unkorreliert (also wegen der Normalvertei-lung auch unabhangig). Dann ist die bedingte Posteriori f

(β∗|y∗,X∗, σ2, τ2

)gleich

MVN

([(X∗)>X∗ +

σ2

τ2I

]−1

(X∗)>y∗, σ2

((X∗)>X∗ +

σ2

τ2I

)−1).

Der Parameter λ enstpricht dabei dem Quotienten σ2/τ2.

2. Ungleiche Varianzen der Storvariablen / abhangige Storvariablen

Allgemein:y = Xβ + ε, ε ∼ N(0,Σε)

y ∼ MVN(Xβ,Σy), Σy = Σε

Problem: Spezifikation der Prioriverteilung fur Σε.

Mogliche Auswege sind:

(a) ParametrisiereΣy = σ2Qy

mit Qy bekannt und p(β, σ2) ∝ (σ2)−1. Dieser Fall ist auf das Modell aus Ab-schnitt 4.6.1 reduzierbar, indem man das Modell

Q−1/2y︸ ︷︷ ︸y∗

= Q−1/2X︸ ︷︷ ︸X∗

β +Q−1/2ε︸ ︷︷ ︸ε∗

betrachtet. Man erhalt dann wieder ein homoskedastisches Regressionsmodell inden Großen y∗,X∗, ε∗.

(b) Gewichtete Regression:Σy = diag(σ2w−1

i )1≤i≤n

lasst sich als Spezialfall von (a) auffassen.

109

(c) Korrelationen: Schreibe

Σy = SRS mit S = diag(σ1, . . . , σp)

mit beispielweise

p(σ21, σ

22, . . . , σ

2p) =

p∏j=1

p(σ2j ) und p(σ2

j ) ∼ inv-χ2(νj , σ20j).

Es bleibt die Spezifikation der Korrelationsmatrix. Priori-Spezifikationen musseninsbesondere positive Definitheit gewahrleisten. Eine einfache Variante besteht inder Verwendung von (positiver)

”Aqui-Korrelation”, was zum Beispiel bei Cluster-

daten eine vernunftige Annahme darstellt:

R =

1 ρ · · · ρρ 1 · · · ρ...

.... . .

...ρ ρ · · · 1

,wobei ρ ∼ U(0, 1) eine positive Korrelation erzwingt. Bei Messwiederholungengreift man oft auf eine autoregressive Kovarianzstruktur zuruck. Fur in 1. Ordnungautokorrelierte Residuen

εt = ρεt−1 + Zt, Zt ∼ N(0, σ2),

erhalt man

R =

1 ρ ρ2 · · · ρp−1

ρ 1 ρ · · · ρp−2

ρ2 ρ 1 · · · ρp−3

......

.... . .

...ρp−1 ρp−2 ρp−3 · · · 1

.Andere Zerlegungen basieren auf der Cholesky- oder Spektralzerlegung und sindrelativ komplex.

(d) Ubergang zu Modellen mit Zufallseffekten: Modelle mit Zufallseffekten (linear mi-xed models, generalized linear mixed models) fuhren zu strukturierten, meist pa-rametersparsamen Kovarianzmatrizen.

Aber : Die Modellgleichung andert sich und Σy 6= Σε, d.h. man kommt in eineandere Modellklasse.

4.7 Bayesianisches generalisiertes lineares Modell

Struktur von GLMs: Der Response folgt einer Verteilung aus einer einfachen Exponentialfa-milie (i = 1, . . . , n)

f(yi|θi) = exp

(yiθi − b(θi)

φi

)· c(yi, φi) (4.1)

110

oder

f(yi|θi) = exp

(yiθi − b(θi)a(φi)

)· c(yi, φi),

wobei in vielen Fallen φi ≡ φ (Bernoulli-, Poissonverteilung). Es ist

µi = E[yi|θi] = b′(θi) , Var(yi|θi) = b′′(θi)φi

und θi der kanonische Parameter. Mit einer Linkfunktion g bzw. Responsefunktion h = g−1

gelteg(µi) = ηi = x>i β. (4.2)

Beispiel 4.3 (Logit-Modell). Mit µi = P(yi = 1) ist

f(yi|µi) = µyii (1− µi)1−yi

= exp(yi log(µi) + (1− yi) log(1− µi)

)= exp

(yi log

(µi

1− µi

)︸ ︷︷ ︸

θi

+ log(1− µi)

)

mit

θi = log

(µi

1− µi

)⇔ µi =

exp(θi)

1 + exp(θi).

Dies entspricht (4.1) mit φi = 1, c(yi, φi) = 1,

b(θi) = − log

(1− exp(θi)

1 + exp(θi)

)= − log

(1

1 + exp(θi)

)= log (1 + exp(θi))

und

b′(θi) =1

1 + exp(θi)· exp(θi) = µi.

Als Prioriverteilung fur β in (4.2) eignet sich

β ∼ MVN(β0,B0)

mit B0 > 0 (vgl. Abschnitt 4.5.3 zur multivariaten Normalverteilung bei bekannter Kovari-anzmatrix). β beeinflusst µ = (µ1, . . . , µn)> uber den Pradiktor µ(β) = h(Xβ), wobei dieAuswertung komponentenweise zu verstehen ist. Fur B−1

0 → 0 erhalt man eine nichtinfor-mative Priori.

Uber die Darstellung (4.1) als Exponentialfamile mit kanonischen Parametern erhalt man alsPosterioriverteilung

f(β|y,X) ∝ exp

(−1

2(β − β0)>B−1

0 (β − β0)

)·n∏i=1

exp

(yiθi − b(θi)

φi

)

= exp

(−1

2(β − β0)>B−1

0 (β − β0)

)exp

(n∑i=1

yiθi − b(θi)φi

)

= exp

(−1

2(β − β0)>B−1

0 (β − β0) +

n∑i=1

yiθi − b(θi)φi

).

111

Beispiel 4.4 (Logit-Modell). Im Falle des Logit-Modells erhalt man als Posteriori

f(β|y,X) ∝ exp

(−1

2(β − β0)>B−1

0 (β − β0)

)×

n∏i=1

µi(β)yi(1− µi(β))1−yi

= exp

(−1

2(β − β0)>B−1

0 (β − β0)

)×

n∏i=1

h(x>i β)yi(1− h(x>i β))1−yi

= exp

(−1

2(β − β0)>B−1

0 (β − β0)

)×

n∏i=1

(exp(x>i β)

1 + exp(x>i β)

)yi( 1

1 + exp(x>i β)

)1−yi.

Problem: Der Posteriori-Kern entspricht keinem Kern einer bekannten Verteilung. Die Posteriori-Verteilung ist demnach nicht analytisch zuganglich. Mogliche Auswege sind

1. Approximation oder

2. Exploration der Posteriori durch Generierung von Samples aus der Posteriori.

Wir betrachten im Folgenden Losung 2. Hier gibt es mehrere Moglichkeiten; sehr etabliert istein Vorschlag von Gamerman (1997)1, eine Variante des Metropolis-Hastings-Algorithmus.

4.7.1 Ein MCMC-Algorithmus: Metropolis-Hastings

Zunachst folgt eine Darstellung des Grundproblems, ohne naher auf die zugrundeliegende

mathematische Theorie einzugehen. Bekannt ist, dass fur Xii.i.d∼ π, i = 1, . . . , n, wobei π

eine Verteilung bezeichnet, interessierende Kennzeichen dieser Verteilung (Momente, Dichteetc.) — Existenz vorausgesetzt — durch Simulation von Zufallszahlen gemaß π als Monte-Carlo-Schatzung gewonnen werden konnen, z.B.

E[X] =1

n

n∑i=1

xi.

Dies ist wenig”spannend”, da, wenn π bekannt ist, in der Regel auch der Erwartungswert

zuganglich ist. Angenommen jedoch, man betrachtet eine (nichtlineare) Funktion von X,zum Beispiel g(X) = X2. Dann ist moglicherweise die Dichte der Transformation g(X) nochanalytisch bestimmbar, aber der Erwartungswert komplex zu berechnen. Im Fall, dass Xmehrdimensional ist, kann die analytische Bestimmung derartiger Kenngroßen analytischunmoglich und bei hoherer Dimension mittels numerischer Integration zu instabil sein.

Unter geeigneten Voraussetzungen lasst sich obige Monte-Carlo Schatzung erweitern zu

E[g(X)] =1

n

n∑i=1

g(xi).

(Dies ist ein allgemeines Prinzip, also nicht notwendigerweise bayesianisch, solange es sichbei π nicht zum Beispiel um eine Posterioriverteilung handelt.) Es sei allerdings bemerkt,

1Gamerman (1997): Sampling from the posteriori distribution in generalized linear models. Statistics andComputing 7, pp. 57-68.

112

dass dieses Vorgehen im Vergleich zur exakten Losung mit einem Monte-Carlo Fehler be-haftet ist. Wesentliche Voraussetzung ist zudem, dass Zufallszahlen aus π gezogen werdenkonnen. Verfahren zur Generierung von i.i.d. Zufallszahlen sind beispielsweise das Inversions-verfahren, Rejection Sampling oder Importance Sampling (vgl. Vorlesung ComputerintensiveMethoden). Gerade bei hoherer Dimension sind diese jedoch zum Teil nicht oder nur sehrkompliziert anwendbar.

Eine Alternative stellen Markov Chain Monte Carlo (MCMC)-Verfahren dar. Ziel ist dieGenerierung einer Markov-Kette (X0, . . . , Xn) von (abhangigen!) Zufallszahlen, deren Vertei-lung gegen die interessierende Verteilung konvergiert, d.h. π ist die stationare oder invarianteVerteilung der Markov-Kette. Der Ergodensatz erlaubt dann Schatzungen der Form

E[X] =1

n− burnin

n∑i=burnin+1

xi bzw. E[g(X)] =1

n− burnin

n∑i=burnin+1

g(xi),

wobei x0, . . . , xburnin Werte am Anfang der Sequenz bezeichnen, bevor sich die Kette in derstationaren Verteilung befindet, und die deshalb

”weggeworfen” werden.

Praktische Umsetzung: Starte mit einem Startwert x0 und ziehe dann fur i = 1, . . . , n WerteXi ∼ P (·|Xi−1), wobei P den Markov-Ubergangskern bezeichnet, der nur vom aktuellen Zu-stand der Kette abhangt. An ihn bzw. die Markov-Kette werden die folgenden Anforderungengestellt:

1. Die Markov-Kette ist homogen.

2. Die Markov-Kette ist irreduzibel.

3. Die Markov-Kette ist aperiodisch.

4. Die Markov-Kette ist positiv rekurrent.

Wir betrachten hier Markov-Ketten in diskreter Zeit bei diskretem oder stetigem Zustands-raum, gewohnlich eine Teilmenge des Rp. Fur allgemeine Zustandsraume ist

”mehr Technik”

erforderlich, aber keine”neuen Ideen”. Fur den hier betrachteten Fall ist die Zielverteilung π

immer gegeben.

Univariater Metropolis-Hastings

Wir beschreiben nun den Metropolis-Hastings-Algorithmus (kurz: MH ) zur Generierung einerwie oben beschriebenen Markov-Kette fur den univariaten Fall; dieser Algorithmus enthaltden Gibbs-Sampler als Spezialfall.

Sei π die Dichte der Zielverteilung, aus der wir simulieren mochten, und q eine geeigneteVorschlagsdichte, aus der neue Zustande der Kette generiert werden, d.h.

Xi ∼ q(·|xi−1),

zum Beispiel qxi|xi−1= N(xi−1, 1) oder qxi|xi−1

= U(xi−1 − c, xi−1 + c). Die Vorschlagewerden nicht immer, sondern nur mit einer gewissen Akzeptanzwahrscheinlichkeit α(xi−1, xi)

113

akzeptiert. Fur den MH-Algorithmus hat diese die Gestalt

α(xi−1, xi) = min

(1,

π(xi) · q(xi−1|xi)π(xi−1) · q(xi|xi−1)

).

Wird xi nicht akzeptiert, so setzt man xi ← xi−1, d.h. der alte Zustand wird beibehalten.

Ein wesentlicher Vorteil dieses Verfahrens besteht darin, dass sich die (meist unbekannte)Normalisierungskonstante von π herauskurzt, d.h. der MH-Algorithmus kann auch (bzw.gerade) fur diese Falle angewendet werden. Die Konstruktion von α gewahrleistet, dass dieBedingungen 1. bis 4. eingehalten werden.

Fur q(xi−1|xi) = q(xi|xi−1) reduziert sich der MH-Algorithmus auf den Metropolis-Algorithmusmit

α(xi−1, xi) = min

(1,

π(xi)

π(xi−1)

),

d.h. wenn die Zieldichte an der Stelle xi großer als an xi−1 ist, wird der neue Vorschlag stetsakzeptiert, andernfalls nur im Verhaltnis π(xi)/π(xi−1). Setzt man die Akzeptanzwahrschein-lichkeit konstant gleich eins, erhalt man den Gibbs-Sampler.

Der MH-Algorithmus akzeptiert tendenziell neue Werte in Bereichen mit hoher Dichte (rele-vante Bereiche). Die Akzeptanzwahrscheinlichkeit sollte nicht zu gering sein, um regelmaßigeZustandsanderungen in der Kette zu erhalten. Sie sollte allerdings auch nicht zu hoch sein,d.h. die Varianz der Vorschlagsverteilung sollte nicht zu niedrig sein, damit der Trager von πausreichend gut exploriert wird.

Algorithmus 11 : Univariater Metropolis-Hastings-Algorithmus

Setze Startwert X0.

Fur i = 1, . . . , n:

1. Ziehe Xi aus q(·|xi−1).

2. Ziehe U ∼ U(0, 1); akzeptiere, wenn

U ≤ α(xi−1, xi),

ansonsten setze xi ← xi−1.

Multivariater Metropolis-Hastings

Die Verallgemeinerung auf den multivariaten Fall ist im Prinzip einfach, zum Beispiel mit

q(xi|xi−1) = MVN(xi−1,Σ).

Ein Problem stellt hier die Wahl der”Tuning-Matrix” Σ dar, die die Akzeptanzwahrschein-

lichkeit steuert. Meist ist Σ = diag(σ21, . . . , σ

2p); man startet mehrere Laufe und berechnet die

Akzeptanzrate. Die Varianzen der Vorschlagsdichte werden dann solange variiert, bis”ange-

messene” Akzeptanzraten erreicht werden.

Im Fall bayesianischer GLMs existiert eine Variante, die automatisch brauchbare Vorschlags-dichten berechnet, wie im folgenden Abschnitt beschrieben wird.

114

4.7.2 Metropolis-Hastings mit IWLS-Vorschlagsdichte

Aus der Vorlesung Generalisierte Regression ist das Fisher-Scoring bekannt:

Beispiel 4.5 (Fisher-Scoring beim Logit-Modell). Die Scorefunktion im Logit-Modell (beikanonischer Linkfunktion) hat die Form

s(β) = X>(y − µ(β)).

Die Fisher-Information ist

F (β) = X>diag(µi(β)(1− µi(β))

)X.

Bezeichnet β den ML-Schatzer, so ist

Cov(β) =[X>diag

(µi(β)(1− µi(β))

)X]−1

.

Der Fisher-Scoring Algorithmus hat dann die Form

β(k+1) = β(k) + F−1(β(k))s(β(k))

= β(k) +[X>diag

(µi(β

(k))(1− µi(β(k))))X]−1

X>(y − µ(β(k))).

Allgemein lasst sich das Fisher-Scoring wie folgt umschreiben: Definiere Pseudo-Beobachtungeny = (y1(β), . . . , yn(β))>, wobei

yi(β) = x>i β +D−1i (yi − µi)

mit

Di(β) =∂h(ηi)

∂ηi=∂h(x>i β)

∂x>i βund ηi = xTi β.

Im Spezialfall des Logit-Modells ist

Di(β) = µi(1− µi) = µi(β)(1− µi(β)).

Fasse diese Eintrage zu D = diag(D1, . . . , Dn) zusammen. Definiere weiter

wi(β) = D2i (β)[Var(yi)]

−1 und W = diag(w1(β), . . . , wn(β)).

Im Logit-Modell:

wi(β) =[µi(1− µi)]2

µi(1− µi)= µi(1− µi).

Dann lasst sich das Fisher-Scoring als iterierte kleinste Quadrate-Schatzung (IWLS, iteratively(re)-weighted least squares) schreiben:

β(k+1) = (X>W (k)X)−1X>W (k)y(k).

Aus der Analogie von kleinster Quadrate- und Maximum-Likelihood-Schatzung im Normal-verteilungsfall lasst sich dies interpretieren als

y(k) ∼ MVN(Xβ, (W−1)(k)

).

Bayesianische Version: Kombiniere das Ganze mit der Prioriverteilung β ∼ MVN(β0,B0).Iteriere dazu:

115

1. Aktueller Zustand sei β(t); berechne

y(t) = X>β(t) +D−1(β(t))(y − µ(β)(t)).

2. Ziehe β∗ ∼ MVN(β(t+1),C(t+1)) mit

β(t+1) = (B−10 +X>W (β(t))X)−1 · [B−1

0 β0 +X>W (β(t))y(β(t))],

C(t+1) = (B−10 +X>W (β(t))X)−1.

3. Akzeptiere β∗ mit Wahrscheinlichkeit

α(β(t),β∗) = min

(1,f(β∗|X)

f(β(t)|X)× q(β(t)|β∗)q(β∗|β(t))

),

wobei q(β(t)|β∗) dem Wert der Dichte von

MVN

((B−1

0 +X>W (β∗)X)−1(

B−10 β0 +X>W (β∗)y(β∗)

),(B−1

0 +X>W (β∗)X)−1)

an der Stelle β(t) entspricht und analog q(β∗|β(t)) dem Wert der Dichte von

MVN

((B−1

0 +X>W (β(t))X)−1(

B−10 β0+X>W (β(t))y(β(t))

),(B−1

0 +X>W (β(t))X)−1)

an der Stelle β∗.

4.8 Bayesianische generalisierte lineare gemischte Modelle

Der Pradiktor des GLM aus dem vorherigen Abschnitt wird hier erweitert. Im Folgendenkonzentrieren wir uns auf Cluster- und Longitudinaldaten. Bei letzteren lassen sich die Datenfur ein Individuum i wie folgt strukturieren:

Response Kovariablen Zeitpunkt der Beobachtung

yi1 xi11 . . . xip1 ti1...

...yiTi xi1Ti . . . xipTi tiTi

Dabei sind yi1, . . . , yiTi korreliert, Ti kann variieren und die Beobachtungszeitpunkte konnenvon Individuum zu Individuum variieren. Die Beobachtungszeitpunkte sollten jedoch nichtinformativ fur den Response sein. Die folgende Abbildung zeigt schematisch eine solche Da-tensituation, wie sie zum Beispiel bei einer kontrollierten Studie auftauchen konnte.

116

Diese Art der Daten stellt nicht nur wegen der Abhangigkeit der Beobachtungen eine Her-ausforderung dar. Haufig sind hier zum Beispiel Drop-Outs und damit fehlende Werte. Ofttreten Longitudinaldaten zudem in Kombination mit Survival-Daten auf, zum Beispiel kannyiti , . . . , yiTi der (mit Messfehlern behaftete) Verlauf eines Biomarkers sein. Die Frage ist dann,ob der Biomarkerverlauf prognostisch fur die Uberlebenszeit ist. Dies fuhrt zu sogenanntenjoint models (siehe Modelle fur Longitudinal- und Uberlebenszeitdaten). GLMMs konnen mitdiesem Datentyp gut umgehen, wenn man die sogenannte

”bedingte Unabhangigkeitsannah-

me” akzeptiert:

1. Erweitere den Pradiktor zu

ηit = x>it β︸︷︷︸feste Effekte

+ z>it αi︸︷︷︸zufallige Effekte

in der Annahme, dassαi ∼ MVN(0,Σ).

Dabei ist x>it = (xi1, . . . , xipt) der Kovariablenvektor, und z>it ∈ R1×q kann Kovariablenaus xit enthalten und zum Beispiel die Zeit t selbst.

Beispiel 4.6 (Random Intercept Modell). Sei

ηiti = β0 + β1ti + αi, αi ∼ N(0, σ2α),

dann haben wir fur ein Individuum i: ηiti1...

ηitTi

=

1 ti1...

...1 tiTi

( β0

β1

)+

1...1

αi.

117

2. Wir treffen die bedingte Unabhangigkeitsannahme

yit |= yit′∣∣ αi,β

fur alle t 6= t′. Diese erlaubt die Darstellung der gemeinsamen Verteilung von (yi1, . . . , yiTi)als Produkt der bedingten Verteilungen

f(yi1, . . . , yiTi) =

Ti∏t=1

f(yit|αi).

Ohne diese bedingte Unabhangigkeitsannahme verlieren GLMMs deutlich an Attraktivitat.

Das volle Setup bei n Individuen sieht wie folgt aus:

yi =

yi1...yiTi

, Xi =

xi11 . . . xip1...

xi1Ti . . . xipTi

, Zi =

zi11 . . . ziq1...

zi1Ti . . . ziqTi

, α =

α1...αn

und

ηi = Xiβ +Ziαi.

Zusammenfassend in Matrixnotation ergibt sich:

y1...yn

=

X1...Xn

β+

Z1 0 · · · 00 Z2 · · · 0...

.... . .

...0 · · · 0 Zn

α1

...αn

=

X1 Z1 · · · 0

X2 0 · · ·...

......

. . . 0Xn 0 · · · Zn

βα1...αn

.

Die bayesianische Herangehensweise ist hier im Prinzip wie im GLM mit

β ∼ MVN(β0,B0) und α ∼ MVN(0,diag(Σ, . . . ,Σ)︸ ︷︷ ︸(n·q)×(n·q)

)

118

bzw. (βα

)∼ MVN

β0

0...0

,diag(B0,Σ, . . . ,Σ)︸ ︷︷ ︸(p+nq)×(p+nq)

.

Dabei bezeichnen diag(Σ, . . . ,Σ) bzw. diag(B0,Σ, . . . ,Σ) Blockdiagonalmatrizen.

Bemerkung.

(i) Bei komplexeren Situationen, zum Beispiel bei nicht unabhangigen Individuen, kanndiag(B0,Σ, . . . ,Σ) durch eine nicht-blockweise Matrix ersetzt werden.

(ii) GLMMs sind hochdimensional, wenn n groß ist. Spezielle Algorithmen zur Optimierungsind notwendig.

Zusatzlich wird eine (Hyper-) Prioriverteilung fur Σ konstruiert, weil die αi unbeobachtete,latente Variablen sind, d.h. die αi sind zu behandeln wie die εi im linearen Modell; auch dorthaben wir fur die Varianz eine Priori angenommen.

Die Priori konnte zum Beispiel Σ ∼ inv-Wishartν0(Λ−10 ), also

p(Σ) ∝ |Σ|−(ν0+q+1)/2 exp

(−1

2tr(Σ−1Λ0)

),

mit resultierender Posteriori

f(β,α,Σ|y,X) ∝

n∏i=1

Ti∏j=1

f(yitj |·)

×exp

(−1

2(β − β0)>B−1

0 (β − β0)

)×

|Σ|−n/2 exp

(−1

2

n∑i=1

α>i Σ−1αi

)×

|Σ|−(ν0+q+1)/2 exp

(−1

2tr(Σ−1Λ0)

).

Ein moglicher Algorithmus zur Simulation der Posteriori ist dann ein blockweiser Gibbs-Sampler.

(i) Full-Conditional fur den β-Block:

f(β|α,Σ,y,X) ∝

n∏i=1

Ti∏j=1

f(yitj |·)

· exp

(−1

2(β − β0)>B−1

0 (β − β0)

)lasst sich wie im bayesianischen GLM behandeln, wenn man zusatzlich einen OffsetzTi αi verwendet:

yi(β(t−1)|αi) = x>i β

(t−1) + z>i αi +D−1i [yi − µi(β(t−1),αi)]

119

bzw.yi(β

(t−1)|αi)− z>i αi = x>i β(t−1) +D−1

i [yi − µi(β(t−1),αi)]

Definiert man≈yi(β

(t−1)|αi) = yi(β(t−1)|αi)− z>i αi, so lasst sich der IWLS-Metropolis-

Hastings-Algorithmus aus 4.7.2 zum Ziehen aus dieser Full-Conditional anwenden.

(ii) Full-Conditional fur den αi-Block: Fur i = 1, . . . , n erhalt man

f(αi|β,Σ,y,X) ∝

n∏i=1

Ti∏j=1

f(yitj |·)

exp(−α>i Σ−1αi

).

Dies lasst sich wieder wie ein GLM mit Offset x>i β interpretieren und mit dem IWLS-MH-Algorithmus mit Proposal-Verteilung

MVN(

[Σ−1 + ziWi(α(t−1)i )z>i ]−1ziWi(α

(t−1)i )[yi(α

(t−1)i )− x>i β], [Σ−1 + ziWi(α

(t−1)i )z>i ]−1

)behandeln.

(iii) Full-Conditional fur Σ:

f(Σ|β,α,y,X) ∝ |Σ|−(n+ν0+q+1)/2 exp

(−1

2tr

(Σ−1Λ0 +

n∑i=1

αiα>i

))

entspricht dem Kern einer inv-Wishartν0+n(Λ0 +∑n

i=1αiα>i )-Verteilung (implizite Di-

mension q × q); diese lasst sich direkt mit einem geeigneten Zufallszahlengenerator si-mulieren.

Insgesamt hat man also einen blockweisen Gibbs-Sampler mit

1 + n + 1 = n+ 2β {αi}ni=1 Σ

Blocken. Verwendet man beim Update von β und αi einen Akzeptanzmechanismus, dannhandelt es sich bei dem Gibbs-Sampler genauer gesagt um einen Metropolis-Hastings-within-Gibbs-Algorithmus, da die einzelnen Blocke jeweils mit Metropolis-Hastings erzeugt, am Endedes Durchgangs aber die n+2 Blocke noch einmal (mit Wahrscheinlichkeit 1) formal akzeptiertwerden.

4.9 Hierarchische Modelle

Dieser Abschnitt behandelt ein Beispiel aus Gelman, Carlin, Stern und Rubin (2003) zumTumorrisiko bei Ratten. Durchgefuhrt wurden 71 Experimente i mit jeweils ni Ratten j:

yij =

{1, Ratte entwickelt Tumor,

0, Ratte entwickelt keinen Tumor.

Also ist yi =∑ni

j=1 yij die Anzahl an Ratten in Experiment i, die einen Tumor entwickeln.

Ideen:

120

1. Experiment-spezifische Wahrscheinlichkeiten θi, i = 1, . . . , n, fur Tumorentwicklung be-trachten, potentiell zuruckzufuhren auf Heterogenitat der Ratten, unterschiedliche expe-rimentelle Bedingungen usw.

2. Alle θi stammen aus einer Population, zum Beispiel einer Beta(α, β)-Verteilung.

3. Anstatt 71 Parameter θ1, . . . , θ71 direkt aus den Daten zu schatzen, nehmen wir eineVerteilung fur die θi an.

4. Ohne weitere Information sind θ1, . . . , θ71 als exchangeable zu betrachten, d.h. fur diegemeinsame Priori von θ = (θ1, . . . , θ71) gilt

p(θ|φ) =71∏i=1

p(θi|φ),

wobei φ die Hyperparameter bezeichnet.”Averaging” uber φ liefert die marginale Priori

p(θ) =

∫ ( n∏i=1

p(θi|φ)

)p(φ) dφ.

Da φ nicht bekannt ist, erhalt es eine eigene (Hyper-) Prioriverteilung p(φ), im Beispieleine Priori fur α und β.

Struktur des hierarchischen Modells:

φ = (α, β) ∼ p(α, β)↙ ↙ ↘ ↘θ1 θ2 . . . θ70 θ71

↓ ↓ ↓ ↓y1|θ1, n1 y2|θ2, n2 . . . y70|θ70, n70 y71|θ71, n71

0/20 0/20 . . . 9/24 4/14

Hierarchisches Modell des Beispiels in der top-down Schreibweise:

tauchen in derLikelihood auf

{yij |ni, θi ∼ Binomial(ni, θi)θi|α, β ∼ Beta(α, β)

taucht nicht in derLikelihood auf

{(α, β) ∼ p(α, β)

Die Posteriori fur alle Parameter lautet

f(θ, α, β|y1, . . . , y71)

∝

(71∏i=1

θyii (1− θi)ni−yi)(

71∏i=1

Γ(α+ β)

Γ(α)Γ(β)θα−1i (1− θi)β−1

)· p(α, β).

Fur feste α, β ist die Posteriori von (θ1, . . . , θ71) das Produkt unabhangiger Posterioris, diejeweils einer Beta(αi, βi)-Verteilung folgen mit αi = α+ yi, βi = β + (ni − yi). Genauer:

f(θ|α, β, {yij}) =

71∏i=1

Γ(α+ β + ni)

Γ(α+ yi)Γ(β + ni − yi)· θα+yi−1i (1− θi)β+ni−yi−1.

121

Die marginale Posteriori f(θ|α, β, {yij}) von (α, β) ergibt sich uber die bereits vielfach ver-wandte Formel

f(α, β|{yij}) =f(θ, α, β|{yij})f(θ|α, β, {yij})

.

Im Beispiel erhalt man

f(α, β|{yij}) ∝ p(α, β) ·71∏i=1

Γ(α+ β)

Γ(α)Γ(β)

Γ(α+ yi)Γ(β + ni − yi)Γ(α+ β + ni)

.

”Knackpunkt” ist die nicht-triviale Wahl der Hyperpriori p(α, β): Oft sind keine Anhalts-

punkte fur die Verteilung der Hyperparameter vorhanden. Die Hyperpriori sollte in diesemFall moglichst wenig informativ sein, allerdings gleichzeitig so, dass die Posteriori proper ist.Fur das Beispiel schlagen Gelman, Carlin, Stern und Rubin (2003) vor, die Hyperparameterin den (kompletten) R2 zu transformieren, zum Beispiel durch

(α, β) 7→[logit

(α

α+ β

), log(α+ β)

],

wobei hier logit(α/(α + β)) = log(α/β). Zur Interpretation sei daran erinnert, dassα/(α + β) dem (Priori-) Erwartungswert einer Betaverteilung entspricht und α + β sich alsPriori-Stichprobengroße auffassen lasst. Eine Gleichverteilung auf der transformierten Skalafuhrt jedoch zu einer uneigentlichen Posteriori. Alternativ schlagen obige Autoren folgendePriori vor:

p

(log

(α

β

), log(α+ β)

)∝ αβ(α+ β)−5/2. (4.3)

Diese ergibt sich aus einer Gleichverteilung fur eine Approximation der Priori-Standardab-weichung (α+β)−1/2, die unabhangig mit einer Gleichverteilung auf dem Priori-Erwartungs-wert kombiniert wird, d.h.

p

(α

α+ β, (α+ β)−

12

)∝ 1.

Die Simulation der Posterioriverteilung f(θ|α, β, {yij}) ist dann wie folgt:

1. Ziehe Zufallszahlen aus f(log(α/β), log(α+β)|{yij}): Dazu wird f(α, β|{yij}) gemaß Dich-tetransformationssatz transformiert (unter anderem wird also p(α, β) durch die Hyperprio-ri (4.3) ersetzt) und dann auf einem feinen Gitter berechnet. (log(α/β), log(α+ β))|{yij}kann nun unter Verwendung des CDF-Samplers (Algorithmus 5 auf Seite 95) simuliertwerden.

2. Transformiere die in Schritt 1 gezogenen Zufallszahlen auf die ursprungliche (α, β)-Skalazuruck.

3. Ziehe dann fur i = 1, . . . , 71 θi|α, β, {yij} gemaß einer Beta(α+yi, β+n−yi)-Verteilung.

4.10 Konvergenzdiagnostik

Schwierigkeiten bei statistischer Inferenz entstehen durch iterative Simulation:

122

1. Simulierte Zufallszahlen aus der Posteriori reprasentieren die Zielverteilung eventuell un-zureichend (Problem der Startwertewahl, Einfluss der Startwerte auf spatere Ziehungen).

2. MCMC: Zufallszahlen sind korreliert — die Inferenz ist ungenauer, als wenn die gleicheAnzahl unabhangiger Zufallszahlen verwandt wurde. Stichwort:

”effektive Sitchproben-

große”.

Strategien:

1. Multiple Sequenzen, die stark uber den Parameterraum streuen.

2. Konvergenz-Monitoring.

3. Falls das”Mixing” schlecht ist (der Parameterraum wird unzureichend exploriert, vor-

wiegend Bewegung entlang weniger lokaler Maxima der Zielverteilung usw.), sollte derAlgorithmus geandert werden.

Bezuglich Punkt 2 schlagen Gelman, Carlin, Stern und Rubin (2003) fur skalaren Parametervor:

• Generiere m parallele Sequenzen der Lange n ≥ 2 nach Entfernen der Burnin-Werte,also

{ψij} fur i = 1, . . . , n, j = 1, . . . ,m.

• Berechne die Varianz zwischen den Sequenzen und innerhalb jeder Sequenz: Setze

ψ·j =1

n

n∑i=1

ψij , ψ·· =1

m

m∑j=1

ψ·j

und definiere

B = n · 1

m− 1

m∑j=1

(ψ·j − ψ··)2

sowie

W =1

m

m∑j=1

s2j mit s2

j =1

n− 1

n∑i=1

(ψij − ψ·j)2.

Betrachte dann folgende Schatzung fur die marginale Posteriori-Varianz:

Var+

(ψ|y) =n− 1

nW +

1

nB.

Die Varianz fur ψ wird in der Regel uberschatzt. Die Schatzung ist jedoch furn→∞ unverzerrt. W allein unterschatzt in der Regel die Varianz, aber

limn→∞

E[W ] = Var(ψ|y).

123

• Fuhre nun ein Monitoring fur die Große

R =

√Var

+(ψ|y)

W

mit R→ 1 fur n→∞ durch, d.h. wenn R groß ist, erhalt man potentiell eine Verbes-serung der Inferenz, wenn weitere Simulationen durchgefuhrt werden.

• Also:

(a) Wenn R groß: Lasse Simulation weiter laufen.

(b) Wenn R”nahe” 1: Verwende alle m ·n Werte fur Posteriori-Inferenz. Die optimale

Große von R ist wiederum ein Tuning-Parameter, zum Beispiel R < 1.1.

Erfolg ist allerdings nicht garantiert. Die Methode funktioniert gut bei approximativnormaler marginaler Posteriori, ist jedoch weniger geeignet, wenn Interesse an den ex-tremen Quantilen besteht.

• Effektive Anzahl unabhangiger Ziehungen:

neff = m · n · Var+

(ψ)

W.

Dabei sollte m nicht zu klein sein, da sonst B wiederum schlecht geschatzt wird.

Obige Methode ist im R-Paket coda (convergence diagnostics and output analysis) implemen-tiert.

4.11 Modellwahl und Modellkritik

Modellwahl ist ein weites Feld, auch in Likelihood-basierter Inferenz. Generelle Strategiensind schwierig zu finden, da in der Regel eine Abhangigkeit vom Kontext (Substanzwissen-schaft) und Randbedingungen (zum Beispiel randomisierte Studie oder Beobachtungsstudie)oder auch dem signal-to-noise-ratio (kontrolliertes experimentelles Umfeld oder starke Hete-rogenitat) besteht.

Sensitivitatsanalyse:

• Einfluss der Priori auf Ergebnisse.

• Einfluss des Modells (Likelihood) auf Gute der Vorhersagen. Welche Falle werdenschlecht durch das Modell beschrieben?

Modellwahl:

Typischerweise, zum Beispiel bei einem Regressionsmodell, erfolgt eine Auswahl der Kova-riablen;

”nested versus non-nested” Situation.

124

Ein popularer Vorschlag zur Modellwahl ist das DIC (Deviance Information Criterion). Wiedas AIC und BIC ist das DIC eine asymptotische Approximation und nur anwendbar, fallsdie Posteriori approximativ multivariat normal ist.

Allgemein ist die Devianz fur Daten y und Parameter(-vektor) θ definiert als

D(y, θ) = −2 log f(y|θ) + C(y).

Beurteilt man Modelle nach der Devianz, so ist ein Modell umso besser, je kleiner diese ist.Das DIC kann aus den generierten Samples der MCMC-Simulation berechnet werden. Seiθ1, . . . , θL eine generierte Sequenz. Die erwartete Devianz bezuglich der Posterioriverteilungvon θ ist

E[D(y, θ)|y]

und wird durch

D =1

L

L∑l=1

D(y, θl)

geschatzt. Dies ist ein Maß dafur, wie gut das Modell an die Daten angepasst ist (je kleiner,umso besser).

Die effektive Anzahl der Parameter in einem bayesianischen Modell ist

pD = D −D(θ),

wobei zum Beispiel

θ =1

L

L∑l=1

θl

die Schatzung des Posteriori-Mittelwerts von θ und D(θ) die Devianz, ausgewertet an θ, ist.Fur ein lineares Modell mit Normalverteilungsannahme entspricht pD der Anzahl unrestrin-gierter Parameter im Modell. Das DIC ist dann definiert als

DIC = 2D −D(θ) = D + pD.

Hinweis:

Dprediction

= E[D(yrep, θ)]

= E

[1

n

n∑i=1

(yrepi − E[yrep

i |y])2

],

wobei der Erwartungswert bezuglich der a posteriori pradiktiven Verteilung zu verstehen ist.

Es ist Dprediction ≈ DIC.

125