Bayesiánská analýza V. Nelineární regresní model Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 1 / 42 Obsah tématu 1 Nelineární regresní model 2 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Random Walk Chain M-H algoritmus Metropolis-within-Gibbs 3 Posteriorní predikční p-hodnota 4 Metoda Gelfanda-Deye 5 Predikce 6 Empirická ilustrace Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 2 / 42 Úvod Pátá kapitola z Koop (2003) resp. učebního textu. Doposud: yi = β1 + β2xi2 + . . . + βkxik + i . Cobb-Douglasova produkční funkce: y = α1xβ2 1 · . . . · xβk k . Logaritmování ln(yi ) = β1 + β2 ln(xi2) + . . . + βk ln(xik) + i . Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 3 / 42 CES funkce CES produkční funkce: yi =   k j=1 γjx γk+1 ij   1 γk+1 . Model: yi =   k j=1 γjx γk+1 ij   1 γk+1 + i . Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 4 / 42 Obecně Předpoklady: 1 ∼ N(0N, h−1 IN). 2 Všechny prvky matice X jsou pevná čísla (tj. nenáhodné proměnné) nebo náhodné veličiny nezávislé se všemi prvky vektoru + p(X|λ), kde λ je vektor parametrů, který neobsahuje žádný z ostatních parametrů modelu. Obecně pracujeme s modelem: y = f (X, γ) + . Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 5 / 42 Nelineární regresní model Obsah tématu 1 Nelineární regresní model 2 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Random Walk Chain M-H algoritmus Metropolis-within-Gibbs 3 Posteriorní predikční p-hodnota 4 Metoda Gelfanda-Deye 5 Predikce 6 Empirická ilustrace Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 6 / 42 Nelineární regresní model Věrohodnostní funkce Z normality náhodných složek: p(y|γ, h) = h N 2 (2π) N 2 exp − h 2 {y − f (X, γ)} {y − f (X, γ)} . Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 7 / 42 Nelineární regresní model Apriorní hustota Obecně p(γ, h). Neinformativní varianta: p(γ, h) = 1 h . Uniformní rozdělení pro γ a ln(h). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 8 / 42 Nelineární regresní model Posteriorní hustota Součin věrohodnostní funkce a apriorní hustoty: p(γ, h|y) ∝ p(γ, h) h N 2 (2π) N 2 exp − h 2 {y − f (X, γ)} {y − f (X, γ)} . Marginální hustota pro γ (při neinformativní apriorní hustotě): p(γ|y) ∝ {y − f (X, γ)} {y − f (X, γ)} −N 2 . Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 9 / 42 Metropolis-Hastings algoritmus Obsah tématu 1 Nelineární regresní model 2 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Random Walk Chain M-H algoritmus Metropolis-within-Gibbs 3 Posteriorní predikční p-hodnota 4 Metoda Gelfanda-Deye 5 Predikce 6 Empirická ilustrace Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 10 / 42 Metropolis-Hastings algoritmus Obecný algoritmus 1 Zvolíme počáteční hodnotu, θ(0). 2 Vygenerujeme kandidátský výběr θ∗ ze zvolené kandidátské hustoty q(θ(s−1); θ). 3 Spočítáme akceptační pravděpodobnost (acceptance probability), α(θ(s−1), θ∗). 4 Přiřadíme θ(s) = θ∗ s pravděpodobností α(θ(s−1), θ∗) a θ(s) = θ(s−1) s pravděpodobností 1 − α(θ(s−1), θ∗). 5 Opakujeme Krok 1, 2 a 3 celkem S krát. 6 Spočítáme průměr S výběrů g(θ(1)), . . . , g(θ(S)). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 11 / 42 Metropolis-Hastings algoritmus Akceptační pravděpodobnost 1 Princip: procházet hlavně oblast vysoké hustoty (ale nejen ji) + tendence vracet se do této oblasti. 2 Více zamítat kandidáty z oblasti nízké hustoty. 3 Podoba akceptační pravděpodobnosti: α(θ(s−1) , θ∗ ) = min p(θ = θ∗|y)q(θ∗; θ = θ(s−1)) p(θ = θ(s−1)|y)q(θ(s−1); θ = θ∗) , 1 4 Kandidátská hustota: q (obecně nemusí být symetrická). 5 Potřeba dobré volby kandidátské hustoty! Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 12 / 42 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Obsah tématu 1 Nelineární regresní model 2 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Random Walk Chain M-H algoritmus Metropolis-within-Gibbs 3 Posteriorní predikční p-hodnota 4 Metoda Gelfanda-Deye 5 Predikce 6 Empirická ilustrace Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 13 / 42 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Princip Kandidátské hustoty nezávislé na výběrech. q(θ(s−1); θ) = q∗(θ). Užitečný v případech, kdy existuje vhodná aproximace posteriorní hustoty (kandidátská hustota). Zjednodušení kandidátské hustoty: α(θ(s−1) , θ∗ ) = min p(θ = θ∗|y)q∗(θ = θ(s−1)) p(θ = θ(s−1)|y)q∗(θ = θ∗) , 1 . Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 14 / 42 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Souvislost s importance sampling Pokud definujeme váhy: w(θA ) = p(θ = θA|y) q∗(θ = θA) . q(θ(s−1); θ) = q∗(θ). Akceptační pravděpodobnost: α(θ(s−1) , θ∗ ) = min w(θ∗) w(θ(s−1)) , 1 . Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 15 / 42 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Volba kandidátské hustoty Není obecné pravidlo. Obvykle využití „klasické“ metody maximální věrohodnosti (maximal likelihood). Klasická ekonometrie: ML estimátor asymptoticky normální s kovarianční maticí var(θML) = I(θ)−1 . Informační matice: I(θ) = −E ∂2 ln(p(y|θ)) ∂θ∂θ . Záporný inverzní Hessián: var(θML). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 16 / 42 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Volba kandidátské hustoty – pokračování Pokud velikost vzorku dostatečně velká + relativně neinformativní apriorní hustota ⇒ posteriorní hustota aproximativně normální se střední hodnotou θML a kovarianční maticí aproximativně var(θML). Někdy maximalizace posteriorní hustoty → θmax a kovarianční maticí aproximativně var(θmax) (lepší v případě informativního prioru). q∗(θ) = fN(θ|θML, var(θML)) dobrá × obvyklejší q∗(θ) = ft(θ|θML, var(θML), ν) (důležitost konců hustoty; nejméně tak tlusté jako posteriorní hustota). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 17 / 42 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Další otázky Konvergenční diagnostiky. Problémy: multimodální rozdělení nebo posterior jen v omezené oblasti (např. gama rozdělení). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 18 / 42 Metropolis-Hastings algoritmus Random Walk Chain M-H algoritmus Obsah tématu 1 Nelineární regresní model 2 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Random Walk Chain M-H algoritmus Metropolis-within-Gibbs 3 Posteriorní predikční p-hodnota 4 Metoda Gelfanda-Deye 5 Predikce 6 Empirická ilustrace Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 19 / 42 Metropolis-Hastings algoritmus Random Walk Chain M-H algoritmus Princip Pokud nenalezneme dobrou aproximaci posteriorní hustoty. Kandidátské výběry dle: θ∗ = θ(s−1) + z. z: increment random variable. θ∗ a θ(s − 1) vstupují do vztahu symetricky ⇒ q(θ∗; θ = θ(s−1)) = q(θ(s−1); θ = θ∗) (Metropolis algoritmus). Akceptační pravděpodobnost: α(θ(s−1) , θ∗ ) = min p(θ = θ∗|y) p(θ = θ(s−1)|y) , 1 . Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 20 / 42 Metropolis-Hastings algoritmus Random Walk Chain M-H algoritmus Volba z Obvykle vícerozměrné normální rozdělení. θ(s−1) je střední hodnota. Potřeba kovarianční matice Σ! q(θ(s−1) ; θ) = fN(θ|θ(s−1) , Σ). Problém s akceptační pravděpodobností (v průměru). Není obecné pravidlo, z praxe pro dimenzi blížící se k nekonečnu, 0.23. Jiné pravidlo: zhruba 0.5. Použití konvergenčních diagnostik! Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 21 / 42 Metropolis-Hastings algoritmus Random Walk Chain M-H algoritmus Volba Σ Snadné experimenty pro skalární Σ. Jinak obvyklé Σ = cΩ, kde c skalár a Ω odhad posteriorní kovarianční matice → experimenty s různými c. Nalezení Ω jako odhadu var(θ|y): První způsob: začít s Σ = cIp a najít c pro které neextrémní akceptační pravděpodobnost (0.000001 nebo 0.99999) → hrubý odhad Ω a volba Σ = cΩ a hledání nové hodnoty c (nalézání lepších Ω a následně Σ). Alternativně: Ω = var(ΩML). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 22 / 42 Metropolis-Hastings algoritmus Metropolis-within-Gibbs Obsah tématu 1 Nelineární regresní model 2 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Random Walk Chain M-H algoritmus Metropolis-within-Gibbs 3 Posteriorní predikční p-hodnota 4 Metoda Gelfanda-Deye 5 Predikce 6 Empirická ilustrace Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 23 / 42 Metropolis-Hastings algoritmus Metropolis-within-Gibbs Princip Kombinace Gibbsova vzorkovače a M-H algoritmu. V nelineárním regresním modelu volba neinformativní nebo nezávislé gama apriorní hustoty pro h implikuje p(h|y, γ) jako gama hustotu. p(γ|y, h) není známá hustota ⇒ M-H algoritmus. Následně p(h|y, γ) jako gama rozdělení → platné výběry γ(s) a h(s) pro s = 1, . . . , S. Obecně: systém plně podmíněných posteriorních hustota, kdy z některých výběry snadné, z jiných M-H algoritmem. Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 24 / 42 Posteriorní predikční p-hodnota Obsah tématu 1 Nelineární regresní model 2 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Random Walk Chain M-H algoritmus Metropolis-within-Gibbs 3 Posteriorní predikční p-hodnota 4 Metoda Gelfanda-Deye 5 Predikce 6 Empirická ilustrace Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 25 / 42 Posteriorní predikční p-hodnota Princip Analýza kvality modelu, alternativa k posteriornímu podílu šancí. Rozlišení mezi aktuálně pozorovanými daty y a pozorovatelnými daty y◦ generovanými modelem → y◦ náhodný vektor rozměru N × 1 s funkcí hustoty pravděpodobnosti p(y◦|θ) . Funkce g(·)→ p(g(y◦)|y) jako souhrn veškeré informace z modelu o g(y◦) po shlédnutí dat. Snadno g(y) pro pozorovaná data → pokud g(y) je z odlehlého konce p(g(y◦)|y), potom model nemůže dobře vysvětlovat g(y), tedy g(y) není tím typem charakteristiky dat, které lze přijatelně generovat modelem. Formální získání pravděpodobnosti okrajových oblastí způsobem podobným p-hodnotě z klasické statistiky. Posteriorní predikční p-hodnota je pravděpodobnost modelu generovat datovou sadu, která bude mít extrémnější vlastnosti než ta, kterou skutečně pozorujeme (jednostranná nebo oboustranná p-hodnota). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 26 / 42 Posteriorní predikční p-hodnota Výpočet p(g(y◦)|y) možno spočítat použitím simulačních metod: p(g(y◦ )|y) = p(g(y◦ )|θ, y)p(θ|y)dθ = p(g(y◦ )|θ)p(θ|y)dθ. Díky podmíněnosti vektorem parametrů θ, aktuální data nepřinášejí žádnou dodatečnou informaci o y◦. Posteriorní simulátor poskytuje výběry z p(θ|y) a jsme tak schopni simulovat výběry z p(g(y◦)|θ) pro daný posteriorní výběr parametrů θ. Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 27 / 42 Posteriorní predikční p-hodnota Využití Posteriorní predikční p-hodnotu lze využít dvěma způsoby: 1 Jako měřítko souladu modelu s daty, tedy jaká je pravděpodobnost, že data byly generována dle tohoto modelu. 2 Pro porovnání různých modelů, tedy jestliže jeden model poskytuje posteriorní predikční p-hodnoty výrazně nižší než druhý model, potom je to důkaz v neprospěch tohoto druhého modelu. Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 28 / 42 Posteriorní predikční p-hodnota Příklad pro test normality – úvod Příklad z úvodu (pro i = 1, . . . , N): y◦ i = f (Xi , γ) + i . Z předpokladů na náhodnou složku: p(y◦ |γ, h) = fN(y◦ |f (X, γ), h−1 IN). Zjednodušení pro neinformativní apriorní hustotou (h lze vyintegrovat): p(y◦ |γ) = ft(y◦ |f (X, γ), s2 IN, N), kde s2 = [y − f (X, γ)] [y − f (X, γ)] N . Nalezení příslušného percentilu, který vytváří g(y) v rámci hustoty p(g(y◦)|y). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 29 / 42 Posteriorní predikční p-hodnota Příklad pro test normality – definice Možnosti volby g(·): analýza reziduí z hlediska vlastností (splnění předpokladů). V bayesovském kontextu jsou chyby i pro i = 1, . . . , N dány jako i = yi − f (X, γ) Předpoklad normality = nulovost Skew = √ N N i=1 3 i N i=1 2 i 3 2 Kurt = N N i=1 4 i N i=1 2 i 2 − 3. Nepozorujeme i → na základě posteriorního simulátoru: E[Skew|y] = E    √ N N i=1[yi − f (Xi , γ)]3 N i=1[yi − f (Xi , γ)]2 3 2 y    . Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 30 / 42 Posteriorní predikční p-hodnota Příklad – praktický postup 1 Vezmeme výběr γ(s) využitím posteriorního simulátoru. 2 Vygenerujeme reprezentativní datovou sadu y◦ z p(y◦|γ(s)) užitím odpovídajících vztahů. 3 Definujeme (s) i = yi − f (Xi , γ(s)) pro i = 1, . . . , N a vyhodnotíme hodnotu šikmosti v tomto bodě, abychom získali Skew(s). 4 Definujeme ◦(s) i = y ◦(s) i − f (Xi , γ(s)) pro i = 1, . . . , N a vyhodnotíme hodnotu špičatosti v tomto bodě, abychom získali Skew◦(s). 5 Opakujeme Krok 1, 2, 3 a 4 celkem S krát. 6 Spočítáme průměr S výběrů Skew(1), . . . , Skew(S) pro odhad E[Skew|y]. 7 Spočítáme podíl S vzorků Skew◦(1), . . . , Skew◦(S), které jsou menší než odhad E[Skew|y]. Odhad p-hodnoty, pokud výsledek menší než 0.5 (jinak jedna mínus). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 31 / 42 Metoda Gelfanda-Deye Obsah tématu 1 Nelineární regresní model 2 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Random Walk Chain M-H algoritmus Metropolis-within-Gibbs 3 Posteriorní predikční p-hodnota 4 Metoda Gelfanda-Deye 5 Predikce 6 Empirická ilustrace Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 32 / 42 Metoda Gelfanda-Deye Princip Pro porovnání nevnořených modelů, či vnořených modelů, kdy nelze snadno využít Savage-Dickey density ratio. Obecnější metoda vyhodnocení posteriorního podílu šancí. Porovnání různých funkčních vztahů f (·). Princip: inverzi marginální věrohodnosti pro model Mi , který závisí na vektoru parametrů θ, lze zapsat jako E[g(θ)|y, Mi ] pro specifickou volbu g(·). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 33 / 42 Metoda Gelfanda-Deye Metoda Gelfanda a Deye Nechť p(θ|Mi ) označuje apriorní hustotu, p(y|θ, Mi ) věrohodnostní funkci a p(θ|y, Mi ) posteriorní hustotu pro model Mi definovaný v oblasti Θ. Pokud f (θ) je funkce hustoty pravděpodobnosti definovaná na oblasti obsahující Θ, potom E f (θ) p(θ|Mi )p(y|θ, Mi ) y, Mi = 1 p(y|Mi ) Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 34 / 42 Metoda Gelfanda-Deye Důkaz E f (θ) p(θ|Mi )p(y|θ, Mi ) y, Mi = f (θ) p(θ|Mi )p(y|θ, Mi ) p(θ|y, Mi )dθ = f (θ) p(θ|Mi )p(y|θ, Mi ) p(θ|Mi )p(y|θ, Mi ) p(y|Mi ) dθ = 1 p(y|Mi ) f (θ)dθ = 1 p(y|Mi ) . Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 35 / 42 Metoda Gelfanda-Deye Problémy Na první pohled lze pro jakoukoliv p.d.f f (θ) nastavit: g(θ) = f (θ) p(θ|Mi )p(y|θ, Mi ) . Následné využití výsledků posteriorního simulátoru k odhadu E[g(θ)|y, Mi ]. Pro úspěch potřeba obezřetné volby f (θ). Geweke (1999): f (θ) p(θ|Mi )p(y|θ,Mi ) musí být shora omezena, tedy musí nabývat konečných hodnot pro jakoukoliv volbu θ. Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 36 / 42 Metoda Gelfanda-Deye Geweke (1999) – úvod Strategie definování f (θ) jako funkce normální hustoty pravděpodobnosti s ohraničenými okraji. Důvod: obtížné ověřit, zda-li f (θ) p(θ|Mi )p(y|θ,Mi ) je konečné na okrajích hustoty normálního rozdělení → odříznutím je f (θ) pro tyto potenciálně problematické oblasti nulová. Formálně: θ a Σ jsou odhady E(θ|y, Mi ) a var(θ|y, Mi ) z posteriorní simulace. Pro určitou oblast pravděpodobnosti p ∈ (0, 1) označuje Θ oblast definičního oboru funkce f (θ): Θ = {θ : (θ − θ) Σ−1 (θ − θ) ≤ χ2 1−p(k)}. χ2 1−p(k) je (1 − p) procentní kvantil rozdělení chí-kvadrát s k stupni volnosti (k je počet prvků vektoru θ). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 37 / 42 Metoda Gelfanda-Deye Geweke (1999) – dokončení Geweke doporučuje f (θ) jako funkci vícerozměrné normální hustoty pravděpodobnosti omezené do oblasti Θ: f (θ) = 1 (1 − p)(2π) k 2 |Σ|−1 2 exp − 1 2 (θ − θ) Σ−1 (θ − θ) 1(θ ∈ Θ). 1() je indikační funkce + očekává se, že nejlepších výsledků dosahujem při volbě nízkých hodnot p (např. p = 0.01) (při odhadu marginální věrohodnosti zahrnujeme mnohem více vzorků). Dodatečný náklad zkoušení různých hodnot p je velmi malý. Standarní způsob spočítání NSE – využití pro vyhodnocení přesnosti odhadu marginální věrohodnosti (např. balík BACC). Metoda pro jakýkoliv model: nutný posteriorní simulátor a známé p(θ|Mi ) a p(y|θ, Mi ) (netriviální, někdy známe jen jádrové hustoty). Gewekova implementace jen v případě Θ ∈ Θ (pokud ne, nabízeny další postupy). Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 38 / 42 Predikce Obsah tématu 1 Nelineární regresní model 2 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Random Walk Chain M-H algoritmus Metropolis-within-Gibbs 3 Posteriorní predikční p-hodnota 4 Metoda Gelfanda-Deye 5 Predikce 6 Empirická ilustrace Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 39 / 42 Predikce Princip Využijeme výběry z p(γ|y), které nám poskytuje M-H algoritmus. Tvorba vzorků z podmíněné predikční hustoty p(y∗|y, γ(s)). Lze ověřit: p(y∗ |y, γ) = ft(y∗ |f (X∗ , γ), s2 IT , N). Výběry z y∗ podmíněné vektorem γ lze tedy získat velmi snadno. Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 40 / 42 Empirická ilustrace Obsah tématu 1 Nelineární regresní model 2 Metropolis-Hastings algoritmus Independence Chain M-H algoritmus Random Walk Chain M-H algoritmus Metropolis-within-Gibbs 3 Posteriorní predikční p-hodnota 4 Metoda Gelfanda-Deye 5 Predikce 6 Empirická ilustrace Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 41 / 42 Empirická ilustrace Příklad CES produkční funkce. Bayesiánská analýza (BAAN) V. Nelineární regresní model Podzim 2018 42 / 42