Bayesovská statistika M. Grendár i Abstrakt—Základy bayesovskaj štatistiky sú ilustrované na príklade hádzania mincou. I. ÚVOD BAYESIÁNCOM je každý, kto robí statistické úvahy pomocou Bayesovej vety. Bayesovská statistika je koncepčne velmi jednoduchá. Na príklade hádzania mincou sa pokúsime ilustrovať jej základné prvky. Posudzovanie modelu a bayesovské priemerovanie modelov budú predstavené prostredníctvom regresného modelu. Výpočtová stránka a niektoré ďalšie aspekty bayesovskej štatistiky sú spomenuté len okrajovo. Začnime s nebayesovskou analýzou výsledkov hádzania mincou. II. Nebayesovská štatistika: analýza hodov mincou Modelujme výsledok hodu mincou pomocou náhodnej premennej X ktorá nadobúda hodnoty z dvojprvkovej množiny X = {0,1} s pravdepodobnosťou P (X; 9) = 9x(l - 9f-x, kde parameter 9 je pravdepodobnosť úspechu P (X = 1). Rozdelenie P(X; 9) je bernoulliovské a náhodný výber X" = X\, X2, ■ ■ ■, Xn je známy aj ako bernoulliovská schéma. Nech sa v n-tici hodov mincou vyskytoval úspech n\ krát. Na základe tejto informácie chceme v rámci nebayesovskej štatistiky urobiť bodový a intervalový odhad hodnoty parametra 9, test hypotézy o 9 a predpoveď výsledku budúceho hodu. A. Odhad V nebayesovskej štatistike existuje veľké množstvo metód na odhadovanie parametrov: metóda najväčšej vierohodnosti, momentová metóda, zovšeobecnená momentová metóda, metóda empirickej vierohodnosti, robustné odhadovanie, atď. Najčastejšie používaná je asi metóda najväčšej vierohodnosti (angl., Maximum Likelihood, ML), v ktorej sa za odhad 9ml parametra 9 berie tá hodnota parametra, pri ktorej by daná realizácia výberu z danej parametrickej triedy rozdelení mala najväčšiu hodnotu vierohodnosti: 9ML = arg sup L(6>; x"), eee kde, v tomto prípade, vierohodnostná funkcia (čiže pravdepodobnosť realizácie výberu, chápaná ako funkcia parametra) L(9;x1) = ("J^^l - 6>)™~ni a parametrický priestor 0 = [0,1]. Ľahko sa zistí, že 9ml = —■ V prípade tohto jednoduchého zadania by asi všetky nebayesovské odhadova-cie metódy viedli na tento odhad. Katedra matematiky, FPV UMB, Tajovského 40, 974 01 Banská Bystrica. Ustav matematiky a informatiky, SAV a UMB, Banská Bystrica. Ustav merania SAV, Bratislava. Email: marian.grendar@savba.sk. Podporené VEGA grantom 1/0077/09. Dátum: 14. september 2009. Príklad: Nech n = 20, nx = 5. Potom 9Mh = ^ = ^ = 0.25. Teda, na základe 20-tich hodov mincou, v ktorých sa strana identifikovaná s X = 1 vyskytovala 5 krát, zoberieme ako odhad skutočnej, nám neznámej hodnoty parametra 9 hodnotu 0.25. Tento bodový odhad by teda naznačoval, že minca je nevyvážená. Ako prípravu na bayesovské úvahy poznamenajme, že mince zvyknú mať homogénne zloženie, bývajú vyvážené, takže by 0 a 1 mali padať s rovnakou pravdepodobnosťou. Otázkou je, ako zahrnúť túto mimodátovú informáciu do štatistických úvah, konkrétne do odhadovania parametra 91 Nebayesovská štatistika na to nemá odpoveď, o B. Inferencie Nebayesovský bodový odhad je v prípade bernoulliovskej schémy veľmi jednoduchý. To isté sa už nedá povedať o infe-renciách, teda o konfidenčných intervaloch a testoch. Tvorba konfidenčných intervalov pre 9 stále zamestnáva nebayesov-ských štatistikov, ako o tom svedčí napr. prehľadový článok [21]. Jeden z najobľúbenejších konfidenčných intervalov pre 9 je waldovský 95%-ný interval: §ml ± 1.96#ml(1 — Ôml)/™, ktorý je založený na Centrálnej limitnej vete (alebo, v tomto prípade, ekvivalentne, asymptotickej normalitě ML odhadu). Jeho pravdepodobnosť pokrytia je u výberov o veľkosti 20 iba 80%-ná1. Preto bolo navrhutých viacero korekcií waldovského intervalu. Iný populárny interval je Clopperov Pearsonov interval, zviazaný dualitou s rovnomenným testom. Príklad (pokr.): Clopperov Pearsonov test dáva pre naše dáta p-hodnotu 0.0414, teda v teste nulovej hypotézy H0 : 9 = 0.5 oproti alternatíve H\ : 9 7^ 0.5, nulovú hypotézu, na konvenčnej 5%-tnej hladine významnosti zamietame. A 95%-ný konfidenčný interval je (0.086,0.491). Asymptotický 95%-ný konfidenčný interval je (0.060, 0.440)2. Tieto nebayesovské inferencie vedú k záveru, že minca je nevyvážená, o Konfidenčný interval je jeden z mnohých fcefry-konceptov (angl., pre-data concepf) nebayesovskej štatistiky. Keby sme realizovali veľký (nekonečný) počet 20-tíc hodov peniazom, tak by daný 100(1 — a)%ný konfidenčný interval pokryl neznámu skutočnú hodnotu parametra v 100(1 — a)%-ách prípadov. Či v prípade konkrétnej realizácie 20-tich hodov mincou daný interval pokryl alebo nepokryl skutočnú hodnotu parametra sa ale nedozvieme: buď pokryl, alebo nepokryl. Nie vždy je toto keby-uvažovanie zmysluplné, nevraviac o tom, že v praxi je obyčajne dôležitejšie fceJ-uvažovanie (angl., post-data reasoning): keď máme túto konkrétnu realizáciu, aká je 'Ďalšou črtou tohto intervalu je, že ak 6>ml = 0, tak interval obsahuje, bez ohladu na n, jediný bod: nulu. 2Vypočítané v R. Kód: libray (Hmisc) ; binconf(5, 20, method = ' alľ ). R-kovský zdrojový kód ku všetkým výpočtom z tohto článku sa dá nájsť na www.savbb.sk/~grendar. 2 pravdepodobnosť, že skutočná hodnota parametra leží v tomto konkrétnom intervale? Táto otázka má zmysel v bayesovskej štatistike. Testovanie hypotéz je notoricky problematické, viď napr. [2], [24], [29], [9]. C. Predikcie V prípade náhodného výberu je predpovedanie veľmi jednoduché. Príklad (pokr): Pravdepodobnosť P(X = 1 | x"), že výsledok budúceho hodu mincou bude X = 1 je rovná 9. Aby sme ju odhadli, je prirodzené nahradiť neznáme 9 jeho odhadom 9ml, čo je v tomto prípade 0.25. Teda, podľa nášho odhadu padne v nasledujúcom hode X = 1 s pravdepodobnosťou 1/4. o III. Bayesovský rámec Rámec, v ktorom sa deje bayesovská štatistika je daný vetou Bayesa3. A. Bayesova veta Bayesovský štatistik musí mimodátovú informáciu o parametri 9 E 0 dáta-generujúceho rozdelenia pravdepodobnosti p(X" | 9) vyjadriť vo forme apriórneho rozdelenia pravdepodobnosti p(9). Prior, ako sa apriórne rozdelenie stručnejšie nazýva, vyjadruje statistikovu neistotu o hodnote parametra 9. Po tom ako je získaná realizácia x" výberu X", modifikuje bayesiánec svoju apriórnu informáciu pomocou Bayesovej vety p{9\xn1) = Jep(x^\9)p(9) de' (D a získa tak aposteriórne rozdelenie p(9 | x") parametra 9, pri daných dátach x". Posteriórne rozdelenie (alebo posterior) vyjadruje statistikovu neistotu o 9 po tom, ako boli pozorované dáta. Posterior samozrejme závisí aj od modelu - teda od apriórneho rozdelenia p(9) a dáta-generujúceho rozdelenia P(*ľ I 0). Niektoré bayesovské úvahy sa dajú robiť aj bez menovateľa v Bayesovej vete, ktorý sa zvykne nazývať evidencia (angl., evidence), alebo aj marginálna vierohodnosť, alebo prediktívna hustota. V takom prípade je zvykom písať Bayesovu vetu v stručnejšom tvare: p(6\x?) 0, 9 e [0,1]. Momenty: E9 = ^ľ_, a0 - . A mód(0) - OL~1 Var6> = (a+ßY(a+ß+l)- a+ß-2- Obr. 1. Beta(100, 100) prior. Skombinujeme . Beta(a = 100, (3 = 100) prior p{9) • s vierohodnosťou p(n\ \ 9) oc 9ni(l — 9)n~ni dát x™ v ktorých sa X = 1 vyskytovalo n\ krát, • pomocou Bayesovej vety, • a dostaneme posterior p(x" \ 9). Dá sa uvidieť, že posteriórne rozdelenie je Beta(a+ni = 105, /3+n — n\ = 115). Posterior ktorý sme obdržali patrí do tej istej triedy rozdelení ako prior. Prior môže mať akýkoľvek tvar. Prior, ktorý po bayesovskom skombinovaní s vierohodnostnou funkciou vedie na posterior patriaci do tej istej triedy rozdelení čo prior, sa nazýva konjugovaný (angl., conjugate prior); [9]. Pretože sú bayesovské výpočty s konjugovanými priormi ovela ľahšie, zvyknú sa takéto priory tiež nazývať pohodlnými priormi (angl., convenience priors). Mimo pohodlnosti, neexistujú vo všeobecnosti dôvody prečo by mala byť apriórna informácia vyjadrená konjugovaným priorom. Vierohodnostná funkcia, prior a posterior sú zobrazené4 na obr. 2. Maximum vierohodnosti sa dosahuje v 0.25 čo je 9ml-Prior je koncentrovaný okolo bodu 0.5. Dát je málo (n = 20), preto je aj posterior nimi len minimálne ovplyvnený, o B. Bodové charakteristiky posterioru Posteriórne rozdelenie obsahuje všetku dostupnú informáciu: sumarizuje model a dáta. 4library(LearnBayes); prior = c(100, 100) # params of Beta prior; data = c(5, 15) # 5 successes and 15 failures; triplot(prior, data) 3 Prior Likelihood Posterior 0.0 0.2 0.4 0.6 0.8 1.0 e Obr. 2. Triplot pre Beta(100, 100) prior. V prípade, že je ale nutné zhrnúť posterior do jediného 'bodu', dá sa tak urobiť v rámci bayesovskej teórie rozhodovania; [9]. Rozhodovanie sa deje za neurčitosti. Nesprávne rozhodnutie má za následok straty. Je nutné špecifikovať stratovú funkciu L(9,a), ktorá udáva závislosť straty od toho ako sa líši bayesovský odhad a E 0 od 9. Po tom ako bola špecifikovaná stratová funkcia je prirodzené hľadať také bodové zhrnutie a posterióru (tzv. bayesovský estimátor/odhad), ktoré minimalizuje priemernú posteriórnu stratu / L(e,a)P(e\x^)de. Je Niektoré základné stratové funkcie a príslušné bayesovské odhady: • posteriórna stredná hodnota, minimalizuje kvadratickú stratovú funkciu L(9, a) = (a — Q)2, • posteriórny medián, minimalizuje absolútnu stratovú funkciu L(9,a) = \a — 6\, • posteriórny mód, minimalizuje všetko-alebo-nič (angl., zero-oné) stratovú funkciu. Nie vždy je zmysluplné hľadať bodovú charakteristiku v rámci bayesovskej teórie rozhodovania. V mnohých prípadoch sa za bodový odhad berie posteriórna stredná hodnota. Ak má ale posterior viac než jeden mód, je zrejmé, že v tomto prípade posteriórna stredná hodnota nie je dobrou bodovou charakteristikou. Príklad (pokr.y. Nakoľko posterior je Beta(105,115), tak posteriórna stredná hodnota ie . n^5l._ = 0.47727. Medián j 105+115 posterióru je5 0.47720, a posteriórny mód je 0.47706. o C. Bayesovská robustnosť Do akej miery závisí posterior na voľbe prioru? Tejto otázke je v bayesovských výskumoch venovaná veľká pozornosť. Tvoria náplň toho, čo sa zvykne nazývať bayesovská robustnosť alebo analýza citlivosti (angl., sensitivity analysis), viď [17]. Navrhnuté boli viaceré techniky na zisťovanie citlivosti posteriórnych úvah na priore, ako aj kvantitatívne miery robustnosti, viď napr. [2]. 5qbeta(0.5, 105, 115) Príklad (pokr.y. Obmedzíme sa len na veľmi primitívnu formu analýzy citlivosti. Skúsime zobrať iný prior, ktorý by vyhovoval nášmu apriórnemu presvedčeniu P(0.44 < 0 < 0.56) = 0.9. Napríklad prior daný ako nasledovná zmes hustôt p{6) = 0.9Beta(500, 500) + 0.1Beta(l, 1) je taký. Aj posterior je potom zmes, a to p(0 \ x") = 0.9Beta(505, 515) + 0.1Beta(6,16). Na obrázku 3 je triplot6 prioru, vierohodnost-nej funkcie a posterióru. Prior Likelihood Posterior 0.0 0.2 0.4 0.6 0.8 1.0 e Obr. 3. Triplot pre zmesový prior. Posteriórna stredná hodnota je 0.472861, a teda je len o čosi menšia než hodnota 0.4772727, ktorú sme dostali pri Beta(100,100) priore, o D. Objektívne priory Nie vždy je k dispozícii taká jasná apriórna informácia ako v prípade hádzania mincou. Značná časť bayesovských výskumov je venovaná tomu, ako pretransformovať žiadne alebo veľmi slabé apriórne presvedčenie o skutočnej hodnote parametra do podoby apriórneho rozdelenia. Tento problém je taký zásadný, že rozdeľuje bayesiáncov na dve velké skupiny: subjektivistov a objektivistov. So značnou dávkou zjednodušenia sa dá povedať, že subjektivisti sú presvedčení, že nemá zmysel hľadať akési automatické (nesubjektívne) metódy na konštrukciu priorov. Naopak, objektivisti aktívne navrhujú metódy na prekladanie neznalosti (angl., ignorance) do podoby apriórneho rozdelenia. Jeden z prominentných predstaviteľov objektívnej bayesovskej štatistiky, James Berger, vraví: "V objektívnej bayesovskej analýze sa priory volia tak, aby predstavovali 'neutrálne' znalosti o neznámych". Väčšina bayesiáncov zastáva pragmatický postoj: v prípade modelu s mnohými parametrami sú pre nepodstatné parametre a parametre o ktorých hodnotách sa apriórne toho veľa nevie použité neinformatívne priory, a pre ostatné parametre sa subjektívne špecifikujú informatívne priory. 6x = seq(0, 1, 0 . 001) ; curve(0.9*dbeta(x,505, 515) + 0.1*dbeta(x, 6,6), col = ' reď ) # posterior; curve(0.9*dbeta(x,500, 500) + 0.1*dbeta(x, 1,1), col = 'darkgreen', add = TRUE) # prior; curve(dbeta(x, 6, 16), col = 'blue', add = TRUE) # like; legend (' toprighť , c ("Prior", "Likelihood", "Posterior"), col = c("darkgreen", "blue", "red")) 4 Existuje velké množstvo techník na konštrukciu objektívnych prioro v. Asi najznámejšou je technika navrhnutá Harol-dom Jeffreysom. Jeffreysov prior má tvar p(0) oc yjl(0), kde 1(0) je Fisherova informácia pre parameter 0. Takto skonštruovaný prior spĺňa požiadavku invariantnosti: pri transformácii 0 na g(0) sa Jeffreysova apriórna hustota transformuje tak, ako sa má hustota pravdepodobnosti transformovať, viď napr. [9]. Jefrreysovské priory sú zvyčajne neznormalizovateľné, nazývajú sa preto pseudo-priory (angl., improper priors). Vie sa, že pre mnohé jefrreysovské priory býva výsledný posterior už normalizovateľný, teda riadny (angl., proper prior). Neznormalizovateľnosť mnohých Jeffrey so vských priorov ale spôsobuje problémy pri bayesovskom porovnávaní hypotéz (pozri G). Iné veľmi populárne objektívne priory sú tzv. referenčné priory, viď [19]. V posledných rokoch sa rozbehlo štúdium a aplikácie tzv. slabo informatívnych priorov (angl., weakly informative priors). Príklad (pokr.y. Predpokladajme, že nemáme k dispozícii žiadnu apriórnu informáciu o minciach, netušíme, že mince sa razia strojovo zo značne homogénnych plechov kovu. Ináč povedané, apriórne nevieme o parametre 0 bernoulliovského rozdelenia nič. Vyjadrime našu neznalosť jeffreysovským pri-orom: p{0) oc O'1/2 (í - ô)"1/2, čo je Beta(0.5, 0.5) prior. Posterior je potom Beta(5.5,15.5). Z obrázku 4 je vidieť, že jeffreysovský prior naozaj vnáša do úvah minimum mimodá-tovej informácie. Obr. 4. Triplot pre Jeffreysov prior. Stredná hodnota posterioru je 0.262. Mód posterioru je 0.237. Obe hodnoty sú blízke k ML odhadu, o E. Intervalové charakteristiky posterioru Existuje viacero možných intervalových charakteristík posterioru. V prípade, že posterior je približne symetrický, je rozumnou charakteristikou jeho rozpätia 100(1 — a)%-ný posteriórny interval, mimo ktorého leží pod každým chvostom posterioru 100(a/2)% masy rozdelenia. Takýto bayesovský konfidenčný interval sa nazýva (rovnako)chvostový (angl., equal-tail). V opačnom prípade je lepšou intervalovou charakteristikou vrcholový posteriórny interval (angl., highest posterior density interval, HPD), ktorý pozostáva z takej oblasti s(x") paramet- rického priestoru 0, pre ktorú 1) P(0 G s(x") | X? = x") = 1 - a, 2) akŕ?aeS(x") afl^sK), potom p(0a | X? = x?) > p(0b | X? = x"). Samozrejme, možné sú aj iné konštrukcie intervalových charakteristík posterioru. Príklad (pokr): V prípade Beta(100,100) prioru a našich dát je posterior (t.j., Beta(105,115)) v podstate symetrický. 95%-ný chvostový posteriórny interval je (0.412, 0.543)7. Pri danom priore a dátach môžeme tvrdiť, že s 95%-nou pravdepodobnosťou leží neznáma hodnota parametra 0 v intervale (0.412,0.543). Lapiace, už v roku 1774 analyzoval bernoulliovskú schému pomocou rovnomerného prioru (ktorý je totožný s Beta(l, 1)). Posterior je v takom prípade Beta(ni + í, n — n\ + 1). Pre naše dáta by v takomto prípade 95%-ný posteriórny chvostový interval bol (0.132,0.437). o F. Testovanie v teórii rozhodovania Testovanie hypotéz je možné robiť v rámci bayesovskej teórie rozhodovania. Formuluje sa nulová hypotéza Hq : 0 E 0rj a alternatívna hypotéza H\ : 0 e 0$. Možné sú dve akcie: akcia a0 - prijatie H0, akcia a\ - prijatie H\. Je nutné špecifikovať stratovú funkciu, ktorá vyjadruje veľkosť straty v prípade chyby prvého a chyby druhého druhu, keď sa vykoná akcia a. Najbežnejšou stratovou funkciou je zovšeobecnená 0-1 strata: cn ak 0 e 6>§, 0 ináč, c: ak 0 e 6>0, 0 ináč. Zovšeobecnená je v tom, že strata cn v prípade chyby druhého druhu nemusí byť rovnaká ako strata c\ v prípade chyby prvého druhu. Bayesovský test spočíva v tom, že sa vyberie tá z hypotéz, ktorá aposteriórne spôsobí v priemere menšiu stratu. Posteriórna priemerná strata sa zvykne tiež nazývať riziko (angl., risk). Riziko má v prípade zovšeobecnenej 0-1 straty tvar: r(ar)K) = cn[l-P(íř0 true|x™)], r(ax | x") = ciP(H0 true | x"). Aby sme zamietli nulovú hypotézu musí byť teda r(ai \ x") < r(ao | x"). To je ekvivalentné tvrdeniu, že zamietni H0 : 0 e 60, ak P(0 e 60 | O < C" . ci + cn V opačnom prípade sa Hq prijíma. Ak sú straty u oboch chybných rozhodnutí rovnaké, tak test zamieta nulovú hypotézu ak je jej posteriórna pravdepodobnosť menšia než 1/2. Príklad (pokr.): Chceme testovať hypotézu, že 0 je z 0rj = (0.47, 0.53). Použijeme Beta(100,100) prior p(0). Prior určuje apriórnu pravdepodobnosť nulovej hypotézy p$ = j&o p(0) d0, 7qbeta(0.025, 105, 115); qbeta(0.975, 105, 115) 5 čo je 0.6048. Apriórne sme teda o čosi viac presvedčení o platnosti nulovej hypotézy, než o alternatíve. Predpokladajme že straty sú rovnaké u oboch chybných rozhodnutí. Nakoľko P (6 E 0o | x") = 0.5269, tak nulovú hypotézu nezamietame. Všimnime si, že aposteriórne, po tom ako sme obdržali dáta, naše presvedčenie o platnosti nulovej hypotézy, v porovnaní s apriórnym presvedčením, kleslo, o Za zmienku stojí, že ak je nulová hypotéza bodová, potom je bayesovské testovanie problematické, pretože si vyžaduje priradenie pravdepodobnostnej miery jednoprvkovej množine. Keď už bayesiánci musia testovať bodovú nulovú hypotézu, robia tak pomocou prioru obsahujúceho diracovskú mieru vo vyšetrovanom bode. G. Porovnávanie hypotéz Testovanie v rámci teórie rozhodovania je síce principiálne, ale zato veľakrát príliš zväzujúce. V prípadoch, keď nie je možné alebo žiadúce formulovať stratovú funkciu, má bayesiá-nec možnosť namiesto testovania, hypotézy Hq a H\ porovnať. Robí sa to pomocou veľmi prirodzeného nástroja, ktorý sa nazýva bayesovský faktor (angl., Bayes factof). Ak je prior riadny, potom je možné definovať apriórny pomer šancí (angl., prior odds ratio): Po = IeMd)dd Pi SeM^M' a aposteriórny pomer šancí (angl., posterior odds ratio): Pn,o = IeoP(0\xi)dd Bayesovský faktor v prospech Hq je gp = posterior odds = pn,o/pn,i 01 prior odds Po/pi Nakoľko nás skôr zaujíma či nulovú hypotézu zamietame, než jej potvrdenie, zvykne sa robiť s bayesovským faktorom BF10 proti H0. Už Jeffreys navrhol kalibráciu BF10 podľa ktorej sa určuje, aký silný je dôkaz o neplatnosti Hq. Pred nedávnom bola táto kalibrácia mierne modifikovaná v práci [27]. Uvedená je v Tabuľke 1. BFio dôkaz o neplatnosti Ho 1 až 3 3 až 20 20 až 150 > 150 takmer žiadny mierny silný veľmi silný Mnohí bayesovskí štatistici považujú hypotézy za neužitočnú formu štatistického usudzovania. Napr. v [3] sa hypotézy spomínajú len v súvislosti s kritikou nebayesovskej štatistiky. H. Predikcie Modelovanie sa často robí za účelom predpovedania, robenia inferencií o nových pozorovaniach x, teda, prediktívnych inferencií. Bayesiáncom na to slúži posteriórna prediktívna hustota (angl., posterior predictive density/distribution) p(x\xTi)= / p(x10,xri)p{eixi)de. (2) Tento dôležitý koncept sa nepoužíva len na predikcie ale aj na posúdenie (angl., validation) modelu (t.j., prioru a dáta-generujúceho rozdelenia). Podobne sa zavádza aj apriórna prediktívna hustota (angl., prior predictive distribution) p{x) = j p{o)p{x i e) de, nazývaná tiež marginálna hustota (angl., marginal density). Príklad (pokr.): Posteriórna prediktívna pravdepodonosť toho, že v budúcom výbere o veľkosti m bude y výskytov X = 1 je (viď [1]) (~\ ™\ (m\ B(a + y,í3 + m-y) kde y = 0,1,..., m a B(-,-) je beta funkcia; a, (3 sú parametre posterioru. Prior nech je Beta(100,100), a nech m = 20. Posteriórne prediktívne rozdelenie11 je na obr. 5. ypred Tabulka I Kalibrácia bayesovského faktora BFi0. Príklad (pokr): Prior nech je Beta(100,100); Potom BFi0 proti H0 : 6 e (0.47,0.53) je 1.37410, teda nie je takmer žiadny dôkaz toho, že by Hq bola neplatná, o su = pbeta(0.53, 100, 100); 100); pO = u-1 9up = pbeta(0.53, 105, 115), 115); pnO = up-lp 10pl = 1 - pO; pni = 1 - pnC pnl/pnO/(pl/pO) 1 = pbeta(0.47, 100, lp = pbeta(0.47, 105, ; BF10 = Obr. 5. Posteriórne prediktívne rozdelenie p(y \ a:™). 92%-ný posteriórny prediktívny interval pre y, pri daných dátach a priore je12 [6,13]. Teda, s 92%-nou pravdepodobnosťou bude v budúcich 20-tich hodoch mincou počet úspechov (X = 1) niekde medzi šesť až trinásť. V prípade, že by nás zaujímala pravdepodobnosť toho, že v budúcom jedinom hode nastane úspech, potom by sme 11 library (LearnBayes ) ; ab = c (105, posterior; m = 20; ys = 0:20; pred ys); plot(ys, pred, type = 'h') 12discint (cbind (ys , pred), 0.9) 115) # = pbetap(ab, 6 podľa (3) dostali p{X = 1 | z?) = = čo je priemer posterioru. V prípade rovnomerného prioru by sme dostali známe Laplaceovo pravidlo postupu (angl., rule of succession): p(X = 1 | x") = n^+2 ■ Pre Beta(100,100) prior je pravdepodobnosť, že po sekvencii 20-tich hodov, v ktorých došlo k úspechu 5 krát, nastane ďalší úspech, rovná 0.477 - priemer posterioru. Pripomeňme, že nebayesiánec by túto pravdepodobnosť odhadol ML odhadom, ktorý je 0.25. Pre rovnomerný prior je hľadaná pravdepodobnosť rovná 0.273. o Ako už bolo spomenuté, posteriórna prediktívna pravdepodobnosť sa používa aj na posteriórnu prediktívnu kontrolu modelu, na posúdenie, zhodnotenie modelu. Voľne povedané, ak je pozorovaný údaj v strede posteriórneho prediktívneho rozdelenia, potom je v súlade s fitom modelu. Ak ale nameraná hodnota leží na chvoste prediktívneho posterioru, potom ju model nevystihuje. Na formálne vyhodnotenie kvality modelu sa používa posteriórna prediktívna p-hodnota, viď [3]. Príklad (pokr.): Posteriórna prediktívna pravdepodonosťp(y < 5 |x™) toho, že v nasledovnej 20-tici hodov sa bude X = 1 vyskytovať 5-krát alebo menej je13 0.0392. To by naznačovalo, že náš model (t.j., prior a dáta-generujúce rozdelenie) nedostatočne vystihuje pozorované dáta, viď tiež obr. 5. Nakoľko dáta-generujúce rozdelenie je v našom prípade sotva spochybniteľné, mýliť sa môžme len v priore. Keďže máme ale dočinenia s malým výberom (n = 20) a apriórna informácia je dostatočne silná, nie je v tomto prípade dôvod meniť model. o Použitie prediktívneho posterioru na posúdenie modelu si zaslúži ešte jednu ilustráciu. Príklad: V R-kovej knižnici faraway sa nachádzajú dáta ozone obsahujúce merania hladiny ozónu ako aj deviatich ďalších premenných. Vyberieme z nich šesť (wind, humidity, temp, dpg, vis, doy), a s ich pomocou budeme chcieť v rámci regresného modelu modelovať hladinu ozónu 03. Po fáze budovania regresného modelu, skončíme napríklad pri modeli 03 ' ~ wind + humidity + poly(dpg,2) + poly(doy, 2) + poly(vis,2) + temp:humidity. Párový graf je na obr. 6. Na parametre modelu (3, a2 položíme jeffreysovský prior. Nájdeme posterior: (3 | cr2, y ~ n((3, a2Vfi) a a2 \ y ~ Inv-%2(n — k, s2), kde (3 = (X'X)'1 X'y, Vfi = (X'X)'1, k je počet prediktorov, n počet pozorovaní, s2 odhad variancie chýb. Hoci je v tomto prípade aj posteriórna prediktívna hustota PÍVi I Vi) vyjádřitelná analyticky (viď napr. [3]), stojí za zmienku, že predstava o prediktívnom posteriore I V\) budúcich dát y", sa dá získať aj bez rátania integrálu v (2). Stačí ak vieme generovať nové y z rozdelenia p(y,0\yi). Odhad hľadaného prediktívneho posterioru sa dá dostať napríklad jadrovým vyhladením vygenerovaných dát y. Nechajme si z posterioru vygenerovať 1000 nových dvojíc cr2), a pre každú následne vygenerujeme pri danej matici plánu X, n-ticu nových y. Jadrovo vyhladíme všetkých 1000 n-tíc a do výsledného odhadu prediktívneho posterioru ešte zakreslíme 13sum (pred [1:6]) 0 2(68 -0. 10 0. 05 -0. 05 0 . 10 03M 1 ! « n n ft n oo cSHpMb ™ 00 í _c i 0 -III 1_ 1 1 wind nantmamn n_ flMma "_ >¥m& _0_1 I hl ma_o_ MB 1MB f!_ . 1 0 ! ! -L B" S B BI o o t 8° o humid P Iii M, ! I I ! 1 1! D" Ii"' s 0 0 B I^ÄlV^Bo dpg-1 D 0 Ij s i 0 IS 11 0 0 1 if % doy.1 1 i h" Jj e o o llľ- so OOOCCCO 0 00 000 9 OS 0 OOOOCOOBO pjiJHSHft&BS -in 0 000000 00 ooooo 0 0 0 0 0 0 0 U OB 90 OtCOXOB Q co OS 1 Q1c0 00 KCtBB 0 u u 0osd9300s0 99 QS0 OS 0 EM CO 0 0 0 vis.1 —o- ■con000 o o nu soo o m&L 0 žžŽiggaiiifc5!! ! liijiiľ II TT Iii j e°s lil i i temp:hum 1. 0 2 . 0 3. 0 20 (0 60 80 -0. 10 0. 00 0 . 10 1000 (000 1000 h Obr. 6. Párový graf dát. pomocou zvislých čiar pozície nameraných y". Výsledok14 je na obr. 7. Z obrázka je zrejmé, že model v hrubých rysoch vystihuje dáta, ale zďaleka nie uspokojivo, nakoľko sa väčšina pozorovaných dát nachádza na pravom chvoste prediktívneho posterioru, a nie v strede. 12 3 4 ypred Obr. 7. Prediktívna aposteriórna hustota p(j/" | y") a pozície nameraných Vi- V inej forme je tá ista vada modelu viditeľná na obr. 8, zobrazujúcom graf závislosti priemerov z y-ov nasimulovaných 14 parsim = blinreg ( (oz$03) (1/3), as.matrix(Xra) f m = 1000); predy = blinregpred(as.matrix(Xm), parsim); plot(density(predy)) for (i in 1:330){abline(v = ( (oz$03) [i] ) ~ (1/3))} 7 z prediktívneho posteriori! oproti pozorovaným referenčnú, 45° čiaru15. Vi- sko aj kde p(x™ I M j) je integrovaná vierohodnosť (angl., integrated likelihood) (1/3) yobs = 03 Obr. 8. Graf priemerov z y-ov nasimulovaných z prediktívneho posterioru oproti pozorovaným j/". Z oboch obrázkov vyplýva, že by bolo potrebné modifikovať model tak, aby sa stred prediktívneho posterioru posunul doprava, o I. Priemerovanie modelov Aby sme mohli popísať ďalší užitočný bayesovský koncept - priemerovanie modelov, budeme musieť opustiť hádzanie mincou. Ideálna na to bude opäť regresia. Predpokladajme, že modelujeme regresným modelom náhodnú premennú y pomocou nejakej množiny prediktorov. Býva zvykom hľadať najlepší podmodel, teda podmnožinu množiny prediktorov, ktorá by v nejakom zmysle najlepšie vystihovala chovanie závislej premennej. Obyčajne sa kvalita posudzuje nejakým kritériom (napr. Akaikeho kritériom, viď [6] alebo [8]), ktoré zaručuje parsimóniu, teda rovnováhu medzi počtom prediktorov a tesnosťou fitu. V bayesovskej štatistike na tento účel slúži Bayesovské informačné kritérium (angl., Bayesian Information Criterion, BIC) (viď [6]), ktoré je aproximáciou bayesovského faktora. Snaha o výber 'naj' modelu je ale problematická, pretože výberom jediného modelu sa ignoruje neistota, ktorú máme o jednotlivých modeloch (angl., model uncertainty), a táto neistota veľmi často prevažuje všetky ostatné zdroje neistoty. Bayesovské priemerovanie modelov (angl., Bayesian Model Avergaing, BMA) poskytuje koherentný spôsob ako zobrať neistotu o modeli do úvahy. Vychádza sa pri tom z apriórneho (zvyčajne rovnomerného) rozdelenia p(-) na množine modelov16 M = {Mi,..., M^}. Prior sa následne bayesovský koriguje dátami X" a obdrží sa posterior p(Mj\x?) = p{x^\M3)p{M3) Y!l=1P{^\Ml)p{Mly (4) p(xl I M,) = / p(xl I 93,M3)p{63 I M,) d0j 15 plot ( (oz$03) " (1/3), colMeans(predy)); abline (0,1) 16Množina m môže byť vo všeobecnosti spojitá. (5) a 6j je parameter modelu M j. Posteriórne inferencie a predikcie sú teda založené na množine modelov, a nie na jedinom, akokoľvek optimálnom modeli. Nech je parameter, ktorý nás zaujíma. Posteriórna modelovo-spriemerovaná hustota p( | x") je fe p(4>\xi) = ^2p(4>\xi,mj)p(mj \ xi)- 3=1 No a jej stredná hodnota a variancia sú modelovo-spriemernené posteriórne charakteristiky parametra . Príklad: V R-kovej knižnici MAS S sa nachádzajú dáta UScrime obsahujúce údaje o 47 amerických štátoch, za rok 1960. Údaje boli zhromaždené za cieľom modelovania kriminality. Prediktorov je 15. Po logaritmickej transformácii prediktorov ako aj závislej premennej použijeme knižnicu BMA17 (viď [31]) na nájdenie modelovo-spriemerovaných po-steriórnych hustôt parametrov modelu18. Z nasledovného numerického zhrnutia sa môžeme dozvedieť, že 5 modelov s najvyššou posteriórnou pravdepodobnosťou má kumulatívnu pravdepodobnosť iba 0.3. Stĺpec p ! =0 obsahuje údaje, pre štyri z prediktorov, o tom, aká je (percentuálna) pravdepodobnosť, že daný prediktor je v modeli. V stĺpci EV sa nachádzajú stredné hodnoty modelovo-spriemernenej posteriórnej hustoty, a S D obsahuje směrodatné odchýlky, pre každú premennú. V ďalších stĺpcoch sú bayesovské odhady parametrov, v troch najlepších modeloch. 51 models were selected Best 5 models have cumulative post. probab. 0.3 p!=0 EV SD modeli model2 model3 M 97.5 1.4014 0.532 1.478 1.514 1.605 U2 82.7 0.2709 0.194 0.289 0.322 0.274 LF 3.4 0.0069 0.103 GDP 34.2 0.2006 0.357 . . 0.541 Na obr. 9 sú modelovo-spriemernené posteriórne hustoty zobrazené pre 4 vybrané prediktory. Veľkosť Diracovho impulzu v bode 0 udáva relatívny počet podmodelov (z celkového počtu 51 podmodelov s najvyššou posteriórnou pravdepodobnosťou), v ktorých sa daný prediktor nevyskytoval. Napríklad u prediktora M (percento mužov vo veku 14-24 rokov) je táto pravdepodobnosť (relatívna početnosť) veľmi malá. Krivka zobrazuje modelovo-spriemernenú posteriórnu hustotu, cez podmodely v ktorých sa daný prediktor nachádza. Predstava o tom, ktorý prediktor je ako významný sa dá rýchlo získať aj z BMA koberca (angl., image plot). Na obr. 10 sú 17V knižnici je implementovaná aproximácia integrovanej vierohodnosti pomocou BICu, čo značne znižuje komplexitu výpočtov. V knižnici BAS je na výpočet použité samplovanie bez navracania, z posteriórnej hustoty. 18library (MASS) ; library (BMA) ; data(UScrime); xcrime = UScrime[,-16]; xcrime = log(xcrime[,-2]); ycrime = log(UScrime[,16]); breg = bicreg(xcrime, ycrime); summary(breg, digits = 2); plot(breg, mfrow = c(2,2), include = c(1,10,5,11), include.inter = F); imageplot.bma(breg) 8 O 1 2 3 -0.20.00.20.40.60.81.0 -2 -1 0 1 2 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 Obr. 9. Modelovo-spriemernený posterior pre štyri parametre modelu. na x-ovej osi zoradené modely, zostupne, podľa posteriórnej pravdepodobnosti. Šírka binu vyjadruje veľkosť posteriórnej pravdepodobnosti daného modelu. Farebne je indikovaná po-steriórna pravdepodobnosť (v %), že prediktor (na osi y) je v modeli. Červená19 označuje hodnoty vyššie než 90%, modrá hodnoty20 v intervale (80, 90)%. 1 2 3 4 5 6 1 8 10 12 14 17 20 23 26 30 35 40 46 Model # Obr. 10. BMA koberec Pre porovnanie21, AIC (angl., Akaike Information Criterion) by do najlepšieho podmodelu vybralo prediktory M, Ed, Pol, M.F, Pop, NW, U2, GDP, Ineq, Prob, Time. Mnohé z nich sa v bayesovsky spriemerovanom modeli nenachádzajú. Naj-AIC model ani nepatrí medzi päť aposteriórne najpravdepodobnejších modelov, o Dôležitosť priemerovania modelov si vďaka bayesiáncom uvedomili aj nebayesiánci, viď [6]. /. A ďalšie 1) Inferencie - založené na (1), predikcie a posúdenie modelu - založené na (2), priemerovanie modelov - založené 19V tlačenej podobe tmavo sivá. 20V tlačenej podobe čierna. 21d = cbind (ycrime, xcrime); m = lm(ycrime ~ ., data = d); summary(m); stepAIC(m) na (4) a (5), ako aj ďalšie bayesovské operácie si vyžadujú výpočet integrálov. Bayesiánci sú spolu so štatistickými fyzikmi, od ktorých sa pred pár desaťročiami naučili integrovať mnohorozmerné integrály pomocou markovovského monte carlo simulovania (angl., Markov chain Monte Carlo, McMC), známi ako 'agresívni' výpočtári. Hoci sú základné McMC algoritmy (Metropolisov, Metropolisov-Hastingsov algoritmus, Gibbsov generátor (angl., Gibbs sampler)) pomerne ľahko popísateľné (viď napr. [9]), ich praktické zvládnutie si vyžaduje značnú dávku kumštu. 2) Bayesiánci sú majstri v modelovaní pomocou hierarchických modelov. S týmito modelmi úzko súvisí aj zaujímavá herézia, známa ako empirický bayes. Jednoduchým príkladom hierarchického bayesovského modelu je: Xí\9í ~ n(&i, a2), kde Xi (í = 1,2,..., n) sú nezávislé; a 6i ~ n(/i,r2), kde 6i (í = 1,2,..., n) sú zaměnitelné. Dá sa ukázať, že marginálne rozdelenie Xi je n(/x, a2 + t2). Ak je a2 známe, môžme teda parametre /x, t2 prioru (nazývajú sa tiež hyperparametrami) odhadnúť z dát. Vďaka marginálnemu rozdeleniu je tak možné určiť hyperparametre prioru z dát. V takomto empirickom bayesovskom prístupe sú teda dáta použité dva razy: na určenie hyperparametrov prioru a aj na jeho aktualizáciu (angl., updating). V hierarchickom bayesovskom modelovaní je podstatný predpoklad zameniteľnosti22 (angl., exchangeability) parametrov 6i. Parametre 6\, 62, ■ ■ ., 6]„ sú zameniteľné ak ich rozdelenie je invariantně na ich permutáciu. So zameniteľnosťou úzko súvisí aj zásadná de Finettiho veta. 3) Odhliadnuc od fundamentálnej odlišnosti bayesovského prístupu od všetkých nebayesovských prístupov, je možné chápať bayesovanie ako metódu, pomocou ktorej sa dajú (principiálnym spôsobom) regularizovať ML odhady. Napríklad známy hrebeňový odhad (angl., ridge estimator) parametrov gaussovského regresného modelu, používaný nebayesiánmi v prípade že matica plánu je takmer singulárna, sa dá bayesovsky obdržať pomocou gaussovského prioru, viď napr. [9]. 4) V posledných rokoch sa prudko rozvíja neparametrická bayesovská štatistika. Pravděpodobnostně modely parametrizované konečným počtom parametrov sú veľakrát príliš neflexibilné. V neparametrickom bayesovaní (viď napr. [5]) sa kladie prior na množinu všetkých rozdelení. 5) ABC (angl., Approximate Bayesian Computation) je zaujímavá technika, pomocou ktorej sa dajú robiť bayesovské úvahy v prípade, že vierohodnostná funkcia je numericky nezvládnuteľná. Na to, aby sa získala vzorka hodnôt z posterioru sa v ABC postupuje podľa tohto receptu: a) zvoľte pracovnú (angl., proposaľ) apriórnu hustotu p(6), b) vygenerujte kandidáta 6* z p(6), c) nasimulujte dáta Y™ z p(Y[l | 6*), d) vypočítajte vhodnú sumárnu štatistiku S(yi) pre nasimulované dáta, ako aj pre napozorované dáta S (x"), e) akceptujte 6* v prípade, že vhodná miera vzdialenosti spĺňa 6(S(x?),S(y?))