Bayesiánská analýza
III. Normální lineární regresní model s přirozeně konjugovanou apriorní
hustotou (více vysvětlujících proměnných)
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 1 / 41
Obsah tématu
1 Lineární regresní model
2 Věrohodnostní funkce a apriorní hustota
3 Posteriorní hustota
4 Porovnání modelů
5 Intervaly nejvyšší posteriorní hustoty
6 Predikce
7 Monte Carlo integrace
8 Empirická ilustrace
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 2 / 41
Úvod
Koop (2003) – kapitola 3
Normální lineární regresní model s více vysvětlujícími proměnnými.
Použití maticové algebry.
Přirozeně konjugovaný prior ⇒ analytické výsledky (podobně jako u
modelu s jedinou vysvětlující proměnnou).
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 3 / 41
Lineární regresní model
Obsah tématu
1 Lineární regresní model
2 Věrohodnostní funkce a apriorní hustota
3 Posteriorní hustota
4 Porovnání modelů
5 Intervaly nejvyšší posteriorní hustoty
6 Predikce
7 Monte Carlo integrace
8 Empirická ilustrace
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 4 / 41
Lineární regresní model
LRM v maticovém vyjádření I
k vysvětlujících proměnných xi1, . . . , xik pro i = 1, . . . , N a model:
yi = β1 + β2xi2 + . . . + βkxik + i .
Úrovňová konstanta: xi1 = 1.
Vektory N × 1 a k × 1:
y =







y1
y2
·
·
yN







=







1
2
·
·
N







β =







β1
β2
·
·
βk







Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 5 / 41
Lineární regresní model
LRM v maticovém vyjádření II
Matice vysvětlujících proměnných rozměru N × k
X =







1 x12 · · x1k
1 x22 · · x2k
· · ·· ·
· · ·· ·
1 xN2 · · xNk







Regresní model:
y = Xβ + .
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 6 / 41
Lineární regresní model
Předpoklady
Předpoklady o a X determinují formu věrohodnostní funkce.
Klasické předpoklady (později uvolněny):
1 ∼ N(0N, h−1
IN), kde h = σ−2
;
2 všechny prvky X jsou nenáhodné veličiny.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 7 / 41
Věrohodnostní funkce a apriorní hustota
Obsah tématu
1 Lineární regresní model
2 Věrohodnostní funkce a apriorní hustota
3 Posteriorní hustota
4 Porovnání modelů
5 Intervaly nejvyšší posteriorní hustoty
6 Predikce
7 Monte Carlo integrace
8 Empirická ilustrace
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 8 / 41
Věrohodnostní funkce a apriorní hustota
Věrohodnostní funkce I
Věrohodnostní funkce s využitím OLS odhadů:
ν = N − k,
β = (X X)−1
X y,
s2
=
(y − Xβ) (y − Xβ)
ν
.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 9 / 41
Věrohodnostní funkce a apriorní hustota
Věrohodnostní funkce II
Věrohodnostní funkce:
p(y|β, h) =
1
(2π)
N
2
h
1
2 exp −
h
2
(β − β) X X(β − β)
× h
ν
2 exp −
hν
2s−2
.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 10 / 41
Věrohodnostní funkce a apriorní hustota
Apriorní hustota
Normální-gama rozdělení, kdy β|h z vícerozměrného normálního
rozdělení, h je stále z gama rozdělení:
β|h ∼ N(β, h−1
V ),
h ∼ G(s−2
, ν.)
Přirozeně konjugovaný prior pro β a h:
β, h ∼ NG(β, V , s−2
, ν).
β: k-rozměrný vektor apriorních středních hodnot pro k regresních
koeficientů, β1, . . . , βk.
V : je pozitivně definitní apriorní část kovarianční matice rozměru
k × k.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 11 / 41
Posteriorní hustota
Obsah tématu
1 Lineární regresní model
2 Věrohodnostní funkce a apriorní hustota
3 Posteriorní hustota
4 Porovnání modelů
5 Intervaly nejvyšší posteriorní hustoty
6 Predikce
7 Monte Carlo integrace
8 Empirická ilustrace
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 12 / 41
Posteriorní hustota
Posterior I
Posteriorní hustota – proporcionální součinu věrohodnostní funkce a
apriorní hustoty
Obdobné NLRM s jedinou vysvětlující proměnnou (použití matic).
Posteriorní hustota parametrů:
β, h|y ∼ NG(β, V , s−2
, ν).
Přičemž platí:
V = (V −1
+ X X)−1
,
β = V (V −1
β + X Xβ),
ν = ν + N.
s−2 je implicitně definováno následovně:
νs2
= νs2
+ νs2
+ (β − β) [V + (X X)−1
]−1
(β − β).
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 13 / 41
Posteriorní hustota
Vlastnosti posteriorní hustoty I
Marginální posteriorní hustota pro β má vícerozměrné t-rozdělení:
β|y ∼ t(β, s2
V
2
, ν).
Z definice t-rozdělení:
E(β|y) = β,
var(β|y) =
νs2
ν − 2
V .
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 14 / 41
Posteriorní hustota
Vlastnosti posteriorní hustoty II
Platí intuice z NLRM s jednou vysvětlující proměnnou (jen vektory
resp. matice).
β je vektor.
Matice (X X)−1: podobná role jako skalár 1
x2
i
v jednoduchém LRM.
V : matice rozměru k × k.
Posteriorní střední hodnota parametru β, β: maticově vážený průměr
apriorní a datové informace.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 15 / 41
Posteriorní hustota
Neinformativní prior I
Nastavení ν = 0 a V −1
„blížící“ se hodnotou k nule.
Obvykle V −1
= cIk, kde c je skalár blížící se nule.
Nepravá apriorní hustota:
p(β, h) ∝
1
h
.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 16 / 41
Posteriorní hustota
Neinformativní prior II
Získáme OLS odhady:
β, h|y ∼ NG(β, V , s−2
, ν),
kde
V = (X X)−1
,
β = β,
ν = N,
νs2
= νs2
.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 17 / 41
Porovnání modelů
Obsah tématu
1 Lineární regresní model
2 Věrohodnostní funkce a apriorní hustota
3 Posteriorní hustota
4 Porovnání modelů
5 Intervaly nejvyšší posteriorní hustoty
6 Predikce
7 Monte Carlo integrace
8 Empirická ilustrace
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 18 / 41
Porovnání modelů
Omezení ve tvaru nerovnosti
Předpokládáme omezení ve tvaru nerovnosti:
Rβ ≥ r.
R je známá matice rozměru J × k a r je známý J-rozměrný vektor.
Porovnáváme dva modely:
M1 : Rβ ≥ r,
M2 : Rβ ≥ r.
Značení v definici M2 znamená, že jedno nebo více z J omezení v M1
není splněno.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 19 / 41
Porovnání modelů
Posteriorní podíl šancí I
Snadný výpočet posteriorního podílu šancí, použití neinformativních
priorů není problém:
PO12 =
p(M1|y)
p(M2|y)
=
p(Rβ ≥ r|y)
p(Rβ ≥ r|y)
.
Posteriorní hustota pravděpodobnosti pro β je z vícerozměrného
t-rozdělení ⇒ p(Rβ|y) je z vícerozměrného t-rozdělení.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 20 / 41
Porovnání modelů
Omezení ve tvaru rovnosti
Případ 1: M1 obsahující omezení Rβ = r a M2, u kterého toto
omezení neplatí (vnořené modely – nested models).
Případ 2: M1 : y = X1β(1) + 1 vzhledem k M2 : y = X2β(2) + 2, kde
X1 a X2 jsou matice obsahující i různé vysvětlující proměnné
(nevnořené modely – non-nested).
Oba případy v definicí modelů pro j = 1, 2:
Mj : yj = Xjβ(j) + j.
Porovnání nevnořených modelů zahrnuje y1 = y2
Porovnání vnořených modelů definuje M2 jako neomezenou regresi a
M1 obsahuje omezení Rβ = r (vyžaduje předefinování vysvětlované a
vysvětlujících proměnných).
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 21 / 41
Porovnání modelů
Posteriorní podíl šancí II
Posteriorní podíl šancí – využití marginálních věrohodností.
Marginální věrohodnost:
p(yj|Mj) = cj
|V j|
|V j|
1
2
(νjs2
j )−
νj
2 ,
pro j = 1, 2, kde
cj =
Γ
νj
2 (νjs2
j )
νj
2
Γ
νj
2 π
N
2
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 22 / 41
Porovnání modelů
Posteriorní podíl šancí III
Posteriorní podíl šancí porovnávájící M1 s M2:
PO12 =
c1
|V 1|
|V 1|
1
2
(ν1s2
1)−
ν1
2 p(M1)
c2
|V 2|
|V 2|
1
2
(ν2s2
2)−
ν2
2 p(M2)
.
Posteriorní podíl šancí závisí na apriorním podílu šancí, zvýhodňuje
soulad modelu s daty, koherenci mezi apriorní a datovou informací a
šetrnost pokud jde o počet vysvětlujících proměnných.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 23 / 41
Porovnání modelů
Porovnání modelů – neinformativní priory
Neformální pravidlo: použití neinformativních priorů pro společné
parametry a čistě informativní priory pro všechny ostatní parametry.
Nastavením ν1 = ν2 = 0: posteriorní podíl šancí má rozumnou
interpretaci (soulad modelu s daty, koherenci mezi apriorní a datovou
informací, atd.)
Rozumné užití neinformativního prioru pro přesnost chyby.
Užití neinformativního prioru pro β(j) – vážné problémy v případě
k1 = k2.
Neinformativní prior pro β(j) založený na V −1
j = cIkj a c → 0.
|V j| = 1
ckj
se vzájemně nevyruší.
Jesliže k1 < k2, PO12 je nekonečno, v případě k2 > k1, PO12 se blíží k
nule bez ohledu na data.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 24 / 41
Intervaly nejvyšší posteriorní hustoty
Obsah tématu
1 Lineární regresní model
2 Věrohodnostní funkce a apriorní hustota
3 Posteriorní hustota
4 Porovnání modelů
5 Intervaly nejvyšší posteriorní hustoty
6 Predikce
7 Monte Carlo integrace
8 Empirická ilustrace
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 25 / 41
Intervaly nejvyšší posteriorní hustoty
Přípustné intervaly
Volná analogie ke konfidenčním intervalům (testování hypotéz).
Definice pro jediný regresní koeficient βj.
95% přípustný interval pro βj je interval [a, b] takový že:
p(a ≤ βj ≤ b|y) = 0.95.
Existuje nekonečně mnoho přípustných intervalů.
Např. βj|y ∼ N(0, 1) ⇒ 95% přípustné intervaly [−1.96, 1.96],
[−1.75, 2.33], [−1.64, ∞] atd.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 26 / 41
Intervaly nejvyšší posteriorní hustoty
Highest Posterior Density Intervals
95% HPDI je 95% přípustný interval, který má nejmenší rozsah oproti
jekémukoliv jinému 95% přípustnému intervalu.
Např. [−1.96, 1.96] je nejkratší přípustný interval.
HPDI existuje vždy, když existuje posteriorní hustota
pravděpodobnosti.
Lze jej využít i při neinformativním prioru.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 27 / 41
Predikce
Obsah tématu
1 Lineární regresní model
2 Věrohodnostní funkce a apriorní hustota
3 Posteriorní hustota
4 Porovnání modelů
5 Intervaly nejvyšší posteriorní hustoty
6 Predikce
7 Monte Carlo integrace
8 Empirická ilustrace
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 28 / 41
Predikce
Prediční hustota
Chceme předpovědět:
y∗
= X∗
β + ∗
.
Předpověď je založena na:
p(y∗
|y) = p(y∗
|y, β, h)p(β, h|y)dβdh.
Lze ukázat:
y∗
|y ∼ t(X∗
β, s2
{IT + X∗
V X∗
}, ν).
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 29 / 41
Monte Carlo integrace
Obsah tématu
1 Lineární regresní model
2 Věrohodnostní funkce a apriorní hustota
3 Posteriorní hustota
4 Porovnání modelů
5 Intervaly nejvyšší posteriorní hustoty
6 Predikce
7 Monte Carlo integrace
8 Empirická ilustrace
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 30 / 41
Monte Carlo integrace
Monte Carlo integrace pro NLRM
Porovnání modelů, predikci a posteriorní analýzu týkající se β lze
provést analyticky (netřeba posteriorních simulátorů).
Theorem (Monte Carlo integrace)
Nechť β(s) pro s = 1, . . . , S je náhodný výběr (vzorek) z p(β|y) a g(·) je
funkce a definujme
gS =
1
S
S
s=1
g(β(s)
),
potom gS konverguje k E[g(β)|y] pro S jdoucí k nekonečnu.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 31 / 41
Monte Carlo integrace
Algoritmus
1 Vygenerujeme náhodný výběr β(s) z posteriorní hustoty pro β
použitím generátoru náhodných čísel pro vícerozměrné t-rozdělení.
2 Spočítáme g(β(s)) a zachováme tento výsledek.
3 Opakujeme předchozí kroky S-krát.
4 Spočítáme průměr S náhodných výběrů g(β(1)), . . . , g(β(S)).
Těmito kroky získáme odhad E[g(β)|y] pro jakoukoliv funkci parametrů,
která nás zajímá.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 32 / 41
Monte Carlo integrace
Numerická standardní chyba
Monte Carlo integrace: aproximace E[g(β)|y] (volbou S řídíme
velikost chyby aproximace).
Číselné měřítko chyby aproximace – využitím centrální limitní věty:
√
S {gS − E[g(β)|y]} → N(0, σ2
g ),
pro S jdoucí k nekonečnu, přičemž σ2
g = var[g(β)|y].
Aproximativní 95% konfidenční interval pro E[g(β)|y]
gS − 1.96
σg
√
S
, gS + 1.96
σg
√
S
.
Alternativně lze využít i numerickou standardní chybu (NSE)
σg
√
S
(obsahuje stejnou informaci v kompaktnější podobě).
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 33 / 41
Empirická ilustrace
Obsah tématu
1 Lineární regresní model
2 Věrohodnostní funkce a apriorní hustota
3 Posteriorní hustota
4 Porovnání modelů
5 Intervaly nejvyšší posteriorní hustoty
6 Predikce
7 Monte Carlo integrace
8 Empirická ilustrace
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 34 / 41
Empirická ilustrace
Reálná data
Reálná data o prodejích domů ve Windsdoru, Kanada v roce 1987,
N = 546.
Proměnné:
yi = prodejní cena i-tého domu (v Kanadských dolarech),
xi2 = rozloha i-tého domu (ve čtverečních stopách),
xi3 = počet ložnic v i-tém domě,
xi4 = počet koupelen v i-tém domu,
xi5 = počet pater v i-tém domu.
Hlavní soubory chapter03.m + skript chapter03_neinf.m + řada
dalších podpůrných funkcí.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 35 / 41
Empirická ilustrace
Reálná data
Předpokládámě vliv rozlohy mezi 0 a 20 dolary za stopu čtvereční ⇒
var(β2) = 25; pro ostatní koeficienty volíme var(β3) = 25002 a
var(β4) = var(β5) = 50002.
var(β) =
νs2
ν − 2
V .
„Téměř“ neinformativní prior a informativní prior s s−2 = 4.0 × 10−8,
ν = 5 a
β =







0.0
10
5000
10000
10000







V =







2.40 0 0 0 0
0 6.0 × 10−7 0 0 0
0 0 0.15 0 0
0 0 0 0.60 0
0 0 0 0 0.60







.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 36 / 41
Empirická ilustrace
Apriorní a posteriorní střední hodnoty pro β (směrodatné
odchylky v závorkách)
Prior Posterior
Informativní (Inf. prior) (Neinf. prior)
β1
0
(10000)
−4035.05
(3530.16)
−4009.55
(3590.11)
β2
10
(5)
5.43
(0.37)
5.43
(0.37)
β3
5000
(2500)
2886.81
(1184.93)
2824.61
(1210.43)
β4
10000
(5000)
16965.24
(1708.02)
17105.18
(1728.18)
β5
10000
(5000)
7641.23
(997.02)
7634.90
(1004.34)
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 37 / 41
Empirická ilustrace
Apriorní a posteriorní střední hodnoty pro h (směrodatné
odchylky v závorkách)
Prior Posterior
Informativní (Inf. prior) (Neinf. prior)
Stř. hodnota 4 × 10−8 3.05 × 10−9 3.03 × 10−9
Sm. odchylka 2.53 × 10−8 1.84 × 10−10 1.83 × 10−10
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 38 / 41
Empirická ilustrace
Porovnání modelů zahrnující parametr β
Post. podíl šancí
p(βj > 0|y) 95% HPDI 99% HPDI pro βj = 0
Informativní apriorní hustota
β1 0.13 [-10969.27,2899.17] [-13159.75,5089.64] 4.14
β2 1.00 [4.71,6.15] [4.49,6.38] 0.00
β3 0.99 [559.29,5214.34] [-175.96,5949.58] 0.39
β4 1.00 [13610.20,20320.27] [12550.37,21380.10] 0.00
β5 1.00 [5682.82,9599.65] [5064.16,10218.30] 0.00
Neinformativní apriorní hustota∗
β1 0.13 [-11061.65,3042.54] [-13289.45,5270.34] —
β2 1.00 [4.71,6.15] [4.48,6.38] —
β3 0.99 [446.96,5202.27] [-304.15,5953.38] —
β4 1.00 [13710.50,20499.85] [12638.10,21572.25] —
β5 1.00 [5662.07,9607.73] [5038.84,10230.96] —
* Jako neinformativní apriorní hustota byla použita apriorní hustota blížící se čistě
neinformativní hustotě. Podíl šancí tak bylo možno spočítat, ale není zde uváděn.
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 39 / 41
Empirická ilustrace
Predikční hustota pro dům o rozloze 5000 čtverečních stop.
0 2 4 6 8 10 12 14
x 10
4
0
0.5
1
1.5
2
2.5
x 10
−5
Predikovana cena domu, y*
Hustotapravdepodobnosti
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 40 / 41
Empirická ilustrace
Posteriorní výsledky pro parametr β2 spočítané
alternativním způsobem
Směrodatná
Stř. hodnota odchylka NSE
Analyticky 5.43 0.37 —
Počet replikací
S = 10 5.44 0.27 0.085
S = 100 5.43 0.38 0.038
S = 1000 5.44 0.36 0.011
S = 10000 5.44 0.36 0.004
S = 100000 5.43 0.37 0.001
Bayesiánská analýza (BAAN) III. NLRM s NCP (více) Podzim 2018 41 / 41