AndrzejTorój Wykład8:MCMC.PróbnikGibbsa.Numerycznaocenagęstościbrzegowej EkonometriaBayesowska

(1)

Ekonometria Bayesowska

Wykład 8: MCMC. Próbnik Gibbsa. Numeryczna ocena gęstości brzegowej

Andrzej Torój

Instytut Ekonometrii – Zakład Ekonometrii Stosowanej

(2)

Próbnik Gibbsa Marginal likelihood

Plan wykładu

1 Próbnik Gibbsa

2 Wiarygodność brzegowa modelu – ocena numeryczna

Andrzej Torój Instytut Ekonometrii – Zakład Ekonometrii Stosowanej

(8) Ekonometria Bayesowska 2 / 12

(3)

Plan prezentacji

1 Próbnik Gibbsa

(4)

Próbnik Gibbsa

Próbnik Gibbsa – idea

Czy potrafimy losować z tych rozkładów?

f (x , y ) ∝ e⁻

(xy −µ)2

2σ2 (x 6= 0, y 6= 0)

f (x|y) ∝ e⁻

(xy−µ)2

2σ2 = e

−(^{x −}^µy)²

2(^σy)²

f (y|x) ∝ e⁻

(xy −µ)2

2σ2 = e

−(^{y −}^µ_x)²

2(^σ_x)²

(5)

Próbnik Gibbsa – idea

f (x , y ) ∝ e⁻

(xy −µ)2

2σ2 (x 6= 0, y 6= 0)

f (x|y) ∝ e⁻

(xy−µ)2

2σ2 = e

−(^{x −}^µy)²

2(^σy)²

f (y|x) ∝ e⁻

(xy −µ)2

2σ2 = e

−(^{y −}^µ_x)²

2(^σ_x)²

(6)

Próbnik Gibbsa

Próbnik Gibbsa – idea

f (x , y ) ∝ e⁻

(xy −µ)2

2σ2 (x 6= 0, y 6= 0)

f (x|y) ∝ e⁻

(xy−µ)2

2σ2 = e

−(^{x −}^µy)²

2(^σy)²

f (y|x) ∝ e⁻

(xy −µ)2

2σ2 = e

−(^{y −}^µ_x)²

2(^σ_x)²

(7)

Próbnik Gibbsa – idea

f (x , y ) ∝ e⁻

(xy −µ)2

2σ2 (x 6= 0, y 6= 0)

f (x|y) ∝ e⁻

(xy−µ)2

2σ2 = e

−(^{x −}^µy)²

2(^σy)²

f (y|x) ∝ e⁻

(xy −µ)2

2σ2 = e

−(^{y −}^µ_x)²

2(^σ_x)²

(8)

Próbnik Gibbsa

Próbnik Gibbsa – przypadek dwuwymiarowy

Rozważmy wektor parametrów θ = (θ₁, θ2) o nieznanej gęstości a posteriori p (θ₁, θ2|y ).

1 Wybieramy startową wartość θ₂⁽⁰⁾.

2 Losujemy θ⁽¹⁾₁ z rozkładu warunkowego p

θ1|θ⁽⁰⁾₂ , y

.

3 Losujemy θ⁽¹⁾₂ z rozkładu warunkowego p

θ₂|θ⁽¹⁾₁ , y .

4 Losujemy θ⁽²⁾₁ z rozkładu warunkowego p

θ1|θ⁽¹⁾₂ , y .

5 Losujemy θ⁽²⁾₂ z rozkładu warunkowego p

θ2|θ⁽²⁾₁ , y .

6 Powtarzamy kroki 4 i 5 naprzemiennie S razy, każdorazowo warunkując wynikiem poprzedniego losowania.

7 Otrzymujemy θ⁽¹⁾, θ⁽²⁾, ..., θ^(S).

(9)

Próbnik Gibbsa – przypadek dwuwymiarowy

θ1|θ⁽⁰⁾₂ , y

.

θ₂|θ⁽¹⁾₁ , y .

θ1|θ⁽¹⁾₂ , y .

θ2|θ⁽²⁾₁ , y .

(10)

Próbnik Gibbsa

Próbnik Gibbsa – przypadek dwuwymiarowy

θ1|θ⁽⁰⁾₂ , y

.

θ₂|θ⁽¹⁾₁ , y .

θ1|θ⁽¹⁾₂ , y .

θ2|θ⁽²⁾₁ , y .

(11)

Próbnik Gibbsa – przypadek dwuwymiarowy

θ1|θ⁽⁰⁾₂ , y

.

θ₂|θ⁽¹⁾₁ , y .

θ1|θ⁽¹⁾₂ , y .

θ2|θ⁽²⁾₁ , y .

(12)

Próbnik Gibbsa

Próbnik Gibbsa – przypadek dwuwymiarowy

θ1|θ⁽⁰⁾₂ , y

.

θ₂|θ⁽¹⁾₁ , y .

θ1|θ⁽¹⁾₂ , y .

θ2|θ⁽²⁾₁ , y .

(13)

Próbnik Gibbsa – przypadek dwuwymiarowy

θ1|θ⁽⁰⁾₂ , y

.

θ₂|θ⁽¹⁾₁ , y .

θ1|θ⁽¹⁾₂ , y .

θ2|θ⁽²⁾₁ , y .

(14)

Próbnik Gibbsa

Próbnik Gibbsa – przypadek dwuwymiarowy

θ1|θ⁽⁰⁾₂ , y

.

θ₂|θ⁽¹⁾₁ , y .

θ1|θ⁽¹⁾₂ , y .

θ2|θ⁽²⁾₁ , y .

(15)

Próbnik Gibbsa – uzasadnienie

p (θ1, θ2|y ) = p (θ₁|θ₂, y ) p (θ2|y )

Losowanie (θ₁, θ2) z rozkładu łącznego można zastąpić

losowaniem θ₁ z rozkładu warunkowego (względem θ₂) oraz θ₂ z rozkładu brzegowego.

Nie znamy jednak rozkładu brzegowego! Nasz wybór θ₂⁽⁰⁾ nie jest więc losowaniem.

Nie wpływa to jednak na rozkład, o ile liczba losowań S jest odpowiednio długa (często odrzucamy S₀ pierwszych losowań jako tzw. burn-in i zostawiamy S₁ = S − S0 pozostałych.

(16)

Próbnik Gibbsa

Próbnik Gibbsa – uzasadnienie

p (θ1, θ2|y ) = p (θ₁|θ₂, y ) p (θ2|y )

(17)

Próbnik Gibbsa – uzasadnienie

p (θ1, θ2|y ) = p (θ₁|θ₂, y ) p (θ2|y )

(18)

Próbnik Gibbsa

Próbnik Gibbsa – uzasadnienie

p (θ1, θ2|y ) = p (θ₁|θ₂, y ) p (θ2|y )

(19)

Próbnik Gibbsa – przypadek ogólny

Rozważmy wektor parametrów θ = (θ1, θ₂, ..., θ_K) o gęstości a posteriori p (θ|y ).

1 Wybieramy wektor wartości startowychθ⁽⁰⁾.

θ1|θ⁽⁰⁾₂ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ₂|θ⁽¹⁾₁ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

4 Losujemy θ⁽¹⁾₃ z rozkładu warunkowego p

θ3|θ⁽¹⁾₁ ,θ⁽¹⁾₂ ,θ₄⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

5 Kontynuujemy aż do kroku K, czyli uzyskania całego wektoraθ⁽¹⁾.

6 Powtarzamy kroki 1-5 zθ⁽¹⁾ jako wektorem wartości startowych.

7 Powtarzamy tę sekwencję S razy.

(20)

Próbnik Gibbsa

Próbnik Gibbsa – przypadek ogólny

θ1|θ⁽⁰⁾₂ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ₂|θ⁽¹⁾₁ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ3|θ⁽¹⁾₁ ,θ⁽¹⁾₂ ,θ₄⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

(21)

Próbnik Gibbsa – przypadek ogólny

θ1|θ⁽⁰⁾₂ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ₂|θ⁽¹⁾₁ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ3|θ⁽¹⁾₁ ,θ⁽¹⁾₂ ,θ₄⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

(22)

Próbnik Gibbsa

Próbnik Gibbsa – przypadek ogólny

θ1|θ⁽⁰⁾₂ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ₂|θ⁽¹⁾₁ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ3|θ⁽¹⁾₁ ,θ⁽¹⁾₂ ,θ₄⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

(23)

Próbnik Gibbsa – przypadek ogólny

θ1|θ⁽⁰⁾₂ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ₂|θ⁽¹⁾₁ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ3|θ⁽¹⁾₁ ,θ⁽¹⁾₂ ,θ₄⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

(24)

Próbnik Gibbsa

Próbnik Gibbsa – przypadek ogólny

θ1|θ⁽⁰⁾₂ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ₂|θ⁽¹⁾₁ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ3|θ⁽¹⁾₁ ,θ⁽¹⁾₂ ,θ₄⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

(25)

Próbnik Gibbsa – przypadek ogólny

θ1|θ⁽⁰⁾₂ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ₂|θ⁽¹⁾₁ ,θ₃⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

θ3|θ⁽¹⁾₁ ,θ⁽¹⁾₂ ,θ₄⁽⁰⁾, ...,θ⁽⁰⁾_K , y .

(26)

Plan prezentacji

1 Próbnik Gibbsa

(27)

Wiarygodność brzegowa modelu

Metody klasy MCMC, m.in. próbnik Gibbsa (jak również próbniki zaimplementowane w rstan) pozwalają na symulację rozkładu a posteriori (tzn. wielokrotne losowanie z niego).

Wyzwaniem pozostaje wyznaczenie wiarygodności brzegowej (ang. marginal likelihood) modelu, czyli całki licznika gęstości a posteriori ze względu na wszystkie parametry (lub, równoważnie, mianownika gęstości a posteriori).

Ta stała skalująca jest używana m.in. do wyznaczenia czynnika Bayesa dla pary modeli.

Czynnik Bayesa

BF^{(i ,j )}= P y |M^{(i )} P y |M^{(j )}

P y |M^{(j )}

=

˙

Θ

P

y |θ^{(j )}, M^{(j )} P

θ^{(j )}|M^{(j )} d θ^{(j )}

¯

(28)

Metody

Estymator średniej harmonicznej

Wzór Bayesa: P (θ|y ) =P(y |θ)P(θ)

P(y ) → ^P(θ)_{P(y )} = ^{P(θ|y )}_{P(y |θ)}

P (y ) =

1 P(y )

−1

=

¯

Θ

P(θ)d θ P(y )

!−1

=

¯

Θ

P(θ) P(y )

d θ

−1

=

¯

Θ

P(θ|y ) P(y |θ)

d θ

−1

Wyrażenie w nawiasie kwadratowym można interpretować jako wartość oczekiwaną _{P(y |θ)}¹ ze względu na gęstość a posteriori P (θ|y ). Dysponując próbkami z gęstości a posteriori

θ⁽¹⁾, θ⁽²⁾, ..., θ^(S¹⁾, możemy wyznaczyć dla każdej z nich wartość f.

wiarygodności, a wyniki uśrednić (Newton & Raftery, 1994):

P (y ) '

"

1 S₁

S1

X

s=1

1 P y |θ^(s)

#⁻¹

(29)

Krytyka i inne metody

Radford Neal (2008), "Worst Monte Carlo method ever": The bad news is that the number of points required for this

estimator to get close to the right answer will often be greater than the number of atoms in the observable universe.

Estymator średniej harmonicznej miewa w wielu przykładach nieskończoną wariancję. (Polecam replikację przykładów R.

Neala.)

Inne podejścia:

metodaGelfanda-Deya (1994)– zob. Koop (2003), roz.

5.7

losowanie z funkcji ważności

zob. też przeglądy Gronau i in. (2017, 2018)

(30)

Metody

Bridge sampling

Meng i Wong (1996). Pakiet bridgesampling w R (dobrze współpracuje z rstan, ale jest uniwersalny).

P (y )=

´P(y |θ)P(θ)h(θ)g (θ)d θ

´ P(y |θ)P(θ)

P(y ) h(θ)g (θ)d θ =

´P(y |θ)P(θ)h(θ)g (θ)d θ´ h(θ)g (θ)P(θ|y )d θ

Licznik to wartość oczekiwana wyrażenia P (y |θ) P (θ) h (θ) ze względu na funkcję generującą kandydatów g (θ) – konceptualny odpowiednik funkcji ważności. Częsty wybór to gęstość (wielowymiarowa) normalna z wartością oczekiwaną i wariancją jak w rozkładzie a posteriori. Losujemy N razy z tej funkcji, i dla wszystkich wyników θ⁽ⁿ⁾ uśredniamy

P

y |θ⁽ⁿ⁾ P

θ⁽ⁿ⁾ h

θ⁽ⁿ⁾ .

Mianownik to wartość oczekiwana wyrażenia h (θ) g (θ) ze względu na gęstość a posteriori. Po przeprowadzeniu symulacji a posteriori metodą MCMC, uśredniamy h

θ^(s) g

θ^(s) . h (θ): ang. bridge function. Odpowiada za to, by próbki z ogonów nie zdominowały zachowania estymatora, i w ten sposób stabilizuje jego wariancję. Meng i Wong (1996)

proponują iteracyjną metodę wyznaczania h (θ) minimalizującą MSE estymatora.