Numerische Methoden der Bayes-Inferenz...

Numerische Methoden der Bayes-Inferenz:Simulationsbasierte Methoden

Helga Wagner Bayes Statistik WS 2010/11 369

Monte-Carlo-Integration

Ist p(x) eine Dichtefunktion, so konnen Integrale der Form

E(g(x)) =

∫g(x)p(x)dx,

mit einer Stichprobe x(1), . . . , x(M) aus p(x) durch den Stichprobenmittelwert

gM =1

M

M∑

i=1

g(x(i))

approximiert werden.

Aus dem starken Gesetz der grossen Zahlen folgt

lim1

M

M∑

i=1

g(x(i)) →a.s.

∫g(x)p(x)dx



Die Varianz des Monte-Carlo-Schatzers gM ist gegeben als

Var(gM) =1

M

∫(g(X)− E(g(x)))2p(x)dx =

1

MVar(g). (39)

soferne Var(g) existiert.

Durch Erhohung von M kann der Approximationsfehler also beliebig verkleinertwerden.

Aus dem zentralen Grenzwertsatz folgt

√M (gM − E(g(y))) ∼ N (0,Var(g)) . (40)



Einen Schatzer fur Var(gM) erhalt man mit

Var(gM) =1

M − 1

M∑

m=1

(g(x(m))− gM)2 ≈ 1

M

M∑

m=1

(g(x(m))− gM)2


Methoden zum Ziehen aus einer Verteilung

Ziehungen x(m), . . . , x(M) aus einer Verteilung P mit Verteilungsfunktion F (x)und Dichtefunktion f(x) konnen mit verschiedenen Verfahren, z.B.

• Inversionsmethode

• Verwerfungsmethode (Accept-Reject-Sampling)

erzeugt werden.


Inversionsmethode

Sei u ∼ U [0, 1], so ist

x = F−1(u) = inf{x : F (x) ≥ u} ∼ P

Beispiel: Zufallszahlen aus der auf das Intervall [a, b] gestutzen Normalverteilung

F (x) =

0 x < aΦ(x)−Φ(a)Φ(b)−Φ(a) a ≤ x ≤ b

1 x > b

Damit istx = Φ−1

(u(Φ(b)− Φ(a)

)+Φ(a)

)∼ F (x)


Verwerfungsmethode

Ziehung aus einer Verteilung mit Dichtefunktion fX(x), indem

• Z aus einer anderen Verteilung mit Dichtefunktion gZ(z) gezogen wird

• Z mit Wahrscheinlichkeit

p =fX(z)

gZ(z)c

akzeptiert wird

Voraussetzung: Es existiert c ≥ 1, sodass

gZ(z)c ≥ fX(z)

fur alle z mit fX(z) > 0.


Verwerfungsmethode

−3 −2 −1 0 1 2 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

fX

c.gZ

akzeptiere

verwirf


Verwerfungsmethode

Fur m = 1, . . . ,M :

• ziehe z ∼ gZ

• ziehe u ∼ U [0, 1].

• Wenn

u ≤ fX(z)

gZ(z)c,

akzeptiere z, d.h. setze x(m) = z und erhohe m um 1.

Sonst verwirf z (d.h. z wird nicht als Ziehung aus der Zielverteilung akzeptiert).


Verwerfungsmethode

Es gilt

P (X ≤ x∗) = P (Z ≤ x∗|U ≤ f(z)

gZ(z)c) =

P ({Z ≤ x∗} ∩ {U ≤ fX(z)cgZ(z)})

P (U ≤ fX(z)cgZ(z))

Wegen

P (U ≤ fX(z)

cgZ(z)) =

∫ ∞

−∞

∫ fX(z)

cgZ(z)

0

dugZ(z)dz =

∫ ∞

−∞

fX(z)

cgZ(z)gZ(z)dz =

1

c

P ({X ≤ x∗} ∩ {U f(z)

gZ(z)c}) =

∫ x∗

−∞

∫ fX(z)

cgZ(z)

0

dugz(z)dz =1

c

∫ x∗

−∞

fX(z)dz

istP (X ≤ x∗) =

∫ x∗

−∞

fX(z)dz


Verwerfungsmethode

• Wahl von gZ, sodass einfaches Ziehen moglich ist und die Hulle cgZ eng ander Zieldichte fX liegt.

• Die Verwerfungsmethode kann auch angewendet werden, wenn die Proportio-nalitatskonstante einer Verteilung nicht bekannt ist:

Sei f∗(x) = kfX(x) und cgZ(x) eine Hulle fur f∗(x), d.h.

cgZ(x) > f∗(x)

fur f∗(x) > 0.

z wird angenommen, wenn

u ≤ f∗(x)(z)

cgZ(z).


Verwerfungsmethode: Beispiel

Ziehen aus der tν (0, 1)-Verteilung:

• Die t1-Verteilung ist die Cauchy (0, 1)-Verteilung, d.h. f(x) = 1π(1+x2)

.

Zufallszahlen konnen mit der Inversionsmethode einfach erzeugt werden:

F (x) =1

2πarctan(y) = u

=⇒x = tan(π(u− 0.5)) ∼ Cauchy (0, 1)

• Fur ν > 1: Verwerfungsmethode mit Z ∼ Cauchy (0, 1)


Verwerfungsmethode: BeispielBestimmen von c, sodass

fX(z)

gZ(z)=

Γ(ν+12 )π(1 + z2)

Γ(ν2)√νπ(1 + z2

ν

)(ν+1)/2≤ c

=⇒ ln(1 + z2)− ν + 1

2ln(1 +

z2

ν) ≤ c∗

Ableiten ergibt2z

1 + z2− (ν + 1)/(2ν)

2z

1 + z2/ν= 0

und damit

(1 + z2/ν)2ν = (ν + 1)(1 + z2) ⇒ ν − 1 = z2(ν − 1)


Verwerfungsmethode: Beispiel

Das Maximum des Quotienten fX(z)gZ(z) ergibt sich fur z = ±1 und damit ist

c =fX(1)

gZ(1)=

2Γ(ν+12 )

Γ(ν2)(1 + 1

ν

)(ν+1)/2

√π

ν

−4 −2 0 2 4

0.0

0.1

0.2

0.3

0.4

x

f(x)


Squeezed Rejection Sampling

Die Zahl der Auswertungen von fX(z) kann mit einer untere Schranke s(z) ≤fX(z) verringert werden.

Seiz ∼ gZ und u ∼ U [0, 1] .

z wird akzeptiert, wenn

u ≤ s(z)

cgZ(x)

andernfalls wird uberpruft, ob u ≤ fX(z)cgZ(x).

Squeezed Rejection Sampling ist dann gunstig, wenn fX(x) aufwandig auszuwer-ten ist.


Adaptive Rejection Sampling

Automatische Erzeugung einer Hulle fur stetige, differenzierbare, log-konkaveDichtefunktion fX (Gilks and Wild, 1992):

• Wahl von Punkten x1 < · · · < xk

• Bestimmen vonHulle: Polygonzug der Tangenten von l(x) = log(fX)

in den Punkten x1, . . . , xk

unterer Schranke: Polygonzug, der die Punkte (x1, l(x1)), . . . , (xk, l(xk))

verbindet.

• Werte z, die im 2. Schritt akzeptiert wurden (d.h. nach Berechnung vonfX(x)) werden als neue Punkte xk+1, . . . hinzugefugt.


Simulationsbasierte Posteriori Inferenz

Mit i.i.d. Ziehungen ϑ(1), . . . ,ϑ(M) aus der Posteriori Verteilung p(ϑ|y) konnenKenngroßen der Posteriori-Verteilung approximiert werden, z.B.

• der Posteriori-Erwartungswert

E(g(ϑ)|y) =∫g(ϑ)p(ϑ|y)dϑ

durch den Mittelwert1

M

M∑

m=1

g(ϑ(m))

• Quantilen der Posteriori-Verteilung durch die entsprechenden Stichproben-quantile

• 100(1− α)%-HPD-Intervalle durch die kurzesten Intervalle, die 100(1− α)%der Stichprobe enthalten


Ziehen aus der Priori-Verteilung

Mit l(ϑ) = p(y|ϑ) ist der Posteriori-Erwartungswert von g(ϑ) gegeben als

E(g(ϑ)|y) =∫g(ϑ)p(ϑ|y)dϑ =

∫g(ϑ)p(y|ϑ)p(ϑ)dϑ∫p(y|ϑ)p(ϑ)dϑ =

E(g(ϑ)l(ϑ))

E(l(ϑ)),

Mit i.i.d.Ziehungen ϑ(1), . . . ,ϑ(M) aus der Priori-Verteilung p(ϑ) ist

E(g(ϑ)|y) ≈1M

∑Mm=1 g(ϑ

(m))p(y|ϑ(m))1M

∑Mm=1 p(y|ϑ(m))

.

Voraussetzung: eigentliche Priori-Verteilung

Monte-Carlo-Integration mit Ziehungen aus der Priori-Verteilung kann sehr inef-fizient sein !


Anwendung auf SFr Wechselkurs-Daten

Modell: yi i.i.d.∼ t4(0, σ2

)

Priori: σ2 ∼ G−1 (c0, C0) mit c0 = 1.5 und C0 = (c0 − 1)s2y(ν − 2)/ν = 0.1361

Monte Carlo Schatzer mit Ziehen aus der Priori-Verteilung: Bestimme

• fur σ2,(m),m = 1, . . . ,M

l(σ2,(m)) = log p(y|σ2,(m)) = −N2log σ2,(m) − 5/2

N∑

i=1

log

(1 +

y2i4σ2,(m)

)

und lmax = maxm l(σ2,(m))

• das Monte-Carlo-Integral

1

M

M∑

m=1

g(σ2,(m)) exp(l(σ2,(m))− lmax)


Anwendung auf SFr Wechselkurs-Daten• g(σ2) = 1: Schatzer fur die Normierungskonstante (marginale Likelihood):

log p(y|M) ≈ lmax + log

(

1

M

M∑

m=1

exp(l(σ2,(m)

) − lmax)

)

.

• Einen Schatzer fur den Posteriori-Erwartungswert erhalt man mit g(σ2) = σ2 aus

E(σ2|y) ≈

(

1

M

M∑

m=1

σ2,(m) exp(l(σ2,(m)) − lmax)

)

exp(lmax)

p(y|M)=

=1M

∑Mm=1 σ

2,(m) exp(l(σ2,(m)) − lmax)

1M

∑Mm=1 exp(l(σ

2,(m)) − lmax)

• Mit g(σ2) = (σ2)2 = σ4 kann ein Schatzer fur die Posteriori-Varianz bestimmt werden

Var(σ2|y) ≈

(

1

M

M∑

m=1

(σ2,(m)

)2exp(l(σ

2,(m)) − lmax)

)

exp(lmax)

p(y|M)− E(σ

2|y)2


SFr Wechselkurs-Daten: Resultate

1 2 3 4

0.285

0.29

0.295

0.3

0.305

0.31

0.315

0.32

0.325

Val

ues

Column Number1 2 3 4

0

2

4

6

8

10

12

x 10−4

Valu

es

Column Number

Abbildung 40: SFr Wechselkurs-Daten, Modell: i.i.d. t4(0, σ2

)mit

Priori σ2 ∼ G−1 (1.5, 0.1361)Schatzung von E(σ2|y) (links) und Var(σ2|y) (rechts) mit Trapezregel(1) bzw. Monte-Carlo Integration mit M = 100 (2), M = 500 (3) andM = 2000 (4) Ziehungen aus der Priori. Die Box-Plots zeigen die Variationfur je 100 Stichproben.


Importance Sampling

Ziehungen aus der Priori kommen aus Bereichen, wo p(ϑ) hoch ist. Eine bessereMonte-Carlo-Schatzung des Integrals ist moglich, wenn Ziehungen vorwiegendaus dem Bereich kommen, wo der Integrand groß ist =⇒ Importance Sampling

Importance Sampling basiert auf folgender Darstellung des Posteriori-Erwartungswertes von g(ϑ)

E(g(ϑ)|y) =∫g(ϑ)p(ϑ|y)

q(ϑ)q(ϑ)dϑ = Eq

(g(ϑ)p(ϑ|y)

q(ϑ)

). (41)


Importance Sampling

Mit M i.i.d. Ziehungen ϑ(1), . . . ,ϑ(M) aus der Verteilung mit Dichte q(ϑ)erhalt man den Schatzer

gISM =1

M

M∑

m=1

g(ϑ(m))p(ϑ(m)|y)q(ϑ(m))

=1

M

M∑

m=1

g(ϑ(m))w(ϑ(m))

Dabei ist

• q(ϑ) die Importance Dichte

• w(ϑ(m)) = p(ϑ(m)|y)

q(ϑ(m))sind die Importance-Gewichte.


Importance Sampling

Zur Berechnung des Importance Schatzers benotigt man die normierte Posteriori-Dichte p(ϑ|y). Die Normierungskonstante kann basierend auf

∫p?(ϑ|y)dϑ =

∫p?(ϑ|y)q(ϑ)

q(ϑ)dϑ = Eq

(p?(ϑ|y)q(ϑ)

)(42)

ebenfalls mit Importance Sampling geschatzt werden.

Die Varianz des Schatzers gISM ist

Var(gISM ) =1

MVarq

(g(ϑ)p(ϑ|y)

q(ϑ)

). (43)


Wahl der Importance-Dichte

q(ϑ) muss so gewahlt werden, dass

Varq

(g(ϑ)p(ϑ|y)

q(ϑ)

)=

∫ (g(ϑ)p(ϑ|y)

q(ϑ)− Ep(g(ϑ))

)2

q(ϑ)dϑ (44)

endlich ist. Eine hinreichende Bedingung dafur ist, dass die Funktion

∣∣∣∣g(ϑ)p(ϑ|y)

q(ϑ)

∣∣∣∣ (45)

nach oben beschrankt ist (Geweke, 1996).

Ist die Varianz endlich, so gilt

√M(gISM − Ep(g(ϑ))

)∼ N

(0,Varq

(g(ϑ)p(ϑ|y)

q(ϑ)

)). (46)


Wahl der Importance-Dichte

• Je kleiner die Varianz Var(gISM ), desto großer ist die Effizienz des Importance-Schatzers. Importance sampling kann mit gut gewahlter Importance Dichtesehr effizient sein.

• Mit q(ϑ) = g(ϑ)p(ϑ|y) hatte der Schatzer Varianz 0.

• Die Varianz wird im wesentlichen vom Verhalten des Quotienten (45) in denEnden der Importance Dichte q(ϑ) bestimmt. Wenn die Enden von q(ϑ)im Vergleich zur Posteriori p(ϑ|y) zu schnell abfallen, dann kann diesesVerhaltnis unbeschrankt sein (z.B. bei normaler Importance Dichte fur eineStudent-t-Posteriori)


Importance Sampling: Anwendung auf SFr Wechselkurs-Daten

• Berechnung von

h(σ2) = logp(y|σ2)p(σ2)

q(σ2)= log p(y|σ2) + log p(σ2)− log q(σ2)

fur Ziehungen σ2,(m),m = 1, . . . ,M aus der Importance-Dichte q(σ2), undbestimmen des Maximuxms hmax = maxm h(y|σ2,(m))

• Berechnung von

1

M

M∑

m=1

g(σ2,(m)) exp(h(σ2,(m))− hmax)

fur g(σ2) = 1, g(σ2) = σ2 und g(σ2) = σ4. Bestimmen der Schatzer analogzum Ziehen aus der Priori-Verteilung.


Importance Sampling: Anwendung auf SFr Wechselkurs-Daten

• Mit g(σ2) = 1 erhalt man die marginale Likelihood

log p(y|M) ≈ hmax + log

(

1

M

M∑

m=1

exp(h(σ2,(m)

) − hmax)

)

.

• Einen Schatzer fur den Posteriori-Erwartungswert erhalt man mit g(σ2) = σ2 aus

E(σ2|y) ≈

(

1

M

M∑

m=1

σ2,(m)

exp(h(σ2,(m)

) − hmax)

)

exp(hmax)

p(y|M)

• Mit g(σ2) = (σ2)2 = σ4 kann ein Schatzer fur die Posteriori-Varianz bestimmt werden

Var(σ2|y) ≈

(

1

M

M∑

m=1

(σ2,(m))2 exp(h(σ2,(m)) − hmax)

)

exp(hmax)

p(y|M)− E(σ2|y)2.


Importance Sampling basierend auf einem approximativen Modell

Modell: yi = µ+ εi, εi ∼ tν(0, σ2

)

Statt des gewunschten Modells Anpassung eines approximativen Modells, daseine konjugierte Analyse erlaubt:

yi = µ+ εi, εi ∼ N(0, σ2

ν

)

mitσ2ν =

ν

ν − 2σ2.

Transformation der Priori-Verteilung fur σ2 ergibt σ2ν ∼ G−1 (c0, C0) mit

c0 = 1.5, C0 = (c0 − 1)s2y.


Importance Sampling basierend auf einem approximativen Modell

Die Posteriori-Verteilung im approximativen Modell ist die inverse Gammavertei-lung σ2

ν|y ∼ G−1 (cN , CN) mit mit Parametern

cN = c0 +N/2, CN = C0 +N∑

i=1

y2i /2

Rucktransformation σ2 = (ν − 2)/νσ2ν ergibt die Importance-Dichte q(σ2):

q(σ2) = fIG(σ2; cN , CN(ν − 2)/ν).


Wahl der Importance Dichte

Mit Information uber Posteriori-Erwartungswert m und Varianz S kann dieImportance-Dichte basierend auf der Normal-Approximation als

q(σ2) = fN(σ2;m,S)

bzw. (robuster) als Dichte der Student-Verteilung mit kleinem Freiheitsgradgewahlt werden

q(σ2) = ftνq(σ2;m,S)

In beiden Fallen ist Stutzung notwendig, um negative Ziehungen zu vermeiden,d.h.

q(σ2) =1

cfN(σ2;m,S)I{σ2>0}, q(σ2) =

1

cftνq(σ

2;m,S)I{σ2>0},

Bestimmung von c analytisch oder als Anteil der positiven Ziehungen.



Eine Alternative ware die Normal- bzw.Student-tVerteilung fur den transformier-ten Parameter log σ2:

q(σ2) =1

σ2fN(log σ2;m,S)

bzw.

q(σ2) =1

σ2ftνq(log σ

2;m,S)

wobei m and S Approximationen fur Mittel und Varianz von θ = log σ2 sind.



Vergleich verschiedener Importance-Schatzer SFr Wechselkurs-Daten: Ver-gleich von Importance-Schatzern mit verschiedenen Importance-Dichten

Importance-Dichten:

• Dichte der Priori-Verteilung (in den Grafiken oben links): Ziehen aus derPriori-Verteilung entspricht Importance Sampling mit Priori-Verteilung alsImportance-Dichte

• Dichte der inverse Gamma-Verteilung aus dem approximativen Modell mitnormalverteilten Fehlern (oben rechts)

• Normal-Approximation (unten links)

• Student-t -Approximation (unten rechts)


Vergleich der Importance Dichte mit der Posteriori-Dichte

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.40

1

2

3

4

5

6x 10

−3

0.24 0.26 0.28 0.3 0.32 0.34 0.360

1

2

3

4

5

6

7

8x 10

−3

0.24 0.26 0.28 0.3 0.32 0.34 0.360

1

2

3

4

5

6

7x 10

−3

0.24 0.26 0.28 0.3 0.32 0.34 0.360

1

2

3

4

5

6x 10

−3

Abbildung 41: SFr Wechselkurs-Daten,Modell: yi i.i.d. ∼ t4

(0, σ2

); Priori: σ2 ∼ G−1 (1.5, 0.1361)

Vergleich verschiedener Importance-Dichten (strichliert) mit Posteriori-Dich-te (voll)


Vergleich der Importance Dichte mit der Posteriori-Dichte

0.24 0.26 0.28 0.3 0.32 0.34 0.360

5

10

15

20

25

30

0.24 0.26 0.28 0.3 0.32 0.34 0.360

2

4

6

8

10

12

14

16x 10

22

0.24 0.26 0.28 0.3 0.32 0.34 0.360

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

0.24 0.26 0.28 0.3 0.32 0.34 0.360

0.2

0.4

0.6

0.8

1

1.2

1.4


(0, σ2

); Priori: σ2 ∼ G−1 (1.5, 0.1361)

Quotient von Importance Dichte und Posteriori


Schatzer fur E(σ2|y)

1 2 3 4

0.285

0.29

0.295

0.3

0.305

0.31

0.315

0.32

0.325

Value

s


0.28

0.285

0.29

0.295

0.3

0.305

0.31

0.315

0.32

0.325

0.33

Value

s

Column Number

1 2 3 4

0.285

0.29

0.295

0.3

0.305

0.31

0.315

0.32

0.325

0.33

Value

s


0.285

0.29

0.295

0.3

0.305

0.31

0.315

0.32

0.325

0.33

Value

sColumn Number


(0, σ2

); Priori: σ2 ∼ G−1 (1.5, 0.1361)

Schatzer fur E(σ2|y): Trapezregel (1) bzw. M = 100 (2), M = 500 (3),M = 2000 (4) Ziehungen aus der Importance Dichte (100 Stichproben)


Schatzer fur Var(σ2|y)

1 2 3 40

1

2x 10

−4

Value

s


0

1

2x 10

−4

Value

s

Column Number

1 2 3 40

1

2x 10

−4

Value

s


0

1

2x 10

−4

Value

sColumn Number


(0, σ2

); Priori: σ2 ∼ G−1 (1.5, 0.1361)

Schatzer fur Var(σ2|y): Trapezregel (1) bzw. M = 100 (2), M = 500 (3),M = 2000 (4) Ziehungen aus der Importance Dichte (100 Stichproben)


SFr Wechselkurs-Daten: Ergebnisse

Bemerkungen:

• Fur die Normal- bzw. Student-Approximation wurden m und S mit 500Ziehungen aus der Priori bestimmt und νq = 10 gewahlt.

• Importance Schatzer mit der Dichte der inversen Gammaverteilung alsImportance-Dichte sind verzerrt (Importance Dichte und Posteriori-Verteilungstimmen zu wenig uberein).

• Normal- bzw. Student-Approximation ergeben bessere Schatzer.

• Die Student-Approximation hat im Vergleich zur Normal-Approximation ge-ringeren Stichprobenfehler.


Numerische Methoden der Bayes-Inferenz...

Documents

Transcript of Numerische Methoden der Bayes-Inferenz...