Kapitola 5 Jádrové odhady dvourozměrných hustot Výstupy z výukové jednotky Student • bude znát součinová a sférická dvourozměrná jádra pro odhady dvourozměrných hustot. • porozumí principu vyhlazování dvourozměrných hustot. • pochopí nejjednodušší metody pro volbu prvků diagonální vyhlazovací matice. • zvládne použití příslušného toolboxu v Matlabu pro simulační studii i pro zpracování reálných dat. 1 Motivace V této kapitole se budeme zabývat rozšířením jádrových odhadů pro jednorozměrné hustoty na odhad vícerozměrných hustot. Ovšem ve vícerozměrném případě nevystačíme s jedním vyhlazovacím parametrem, ale je třeba specifikovat matici vyhlazovacích parametrů. Tato matice řídí jak hladkost, tak i orientaci vícerozměrného vyhlazení. Budeme se zabývat jádrovým odhadem, který je přímým rozšířením jednorozměrného odhadu (3.1) v kapitole 3, a zaměříme se zejména na odhad dvourozměrné hustoty. 3,-.---g-.-1 3 g I-,-,-,-1 _gl-,-,-1— 2 1 0 1 2-2-1 0 1 Obrázek 5.1: Náhodný výběr a jeho jádrový odhad 60 Poznámka 1.1. Jádrové odhady dvourozměrných hustot se obvykle znázorňují pomocí vrstevnic, které umožňují snazší náhled na odhadnutou funkci. 2 Základní typy odhadů Podobně jako u odhadů hustoty můžeme použít histogram, ale ten má zmíněné nevýhody - jde o schodovitou funkci a je citlivý na volbu počtu a šířky třídicích obdélníků - viz obrázek 5.2. Příklad 2.1. Mějme dán datový soubor o velikosti n — 100 generovaný ze směsi tří normálních hustot1. N(0, -1; 1/3,1/3, 0), N(0, 2; 1,1, 0) a N(0,4; 1/3,1/3, 0) f(x v) = lJLe-l(*2+(y+tr)+l±_ e-H-2+fe-2)2)+IAe-!(-2+fe-4)2) Jy 3 2tt 3 2tt 3 2tt (Data jsou v tabulce 6.4.) Z obrázku 5.2 je patrné, že histogram nepostihuje charakteristické rysy hustoty pravděpodobnosti dat. (a) 7 x 7 obdélníků (b) 5 x 10 obdélníků Obrázek 5.2: Histogramy s různými počty třídicích obdélníků Předpokládejme, že máme k dispozici náhodný výběr ([Xi, Yi],..., [Xn, Yn]) z dvourozměrného spojitého rozdělení s hustotou f(x,y). Jádrový odhad hustoty / v bodě [x,y] G M2 je definovaný vztahem i n i n f(x, y; h) = - Kn(x -XllV- Y), = ^ ]T ^(h^ -XllV- y,)T) (5.1) n i=i n' ' i=i přičemž h je matice vyhlazovacích parametrů a K je dvourozměrné jádro. Jádro K je dvourozměrná funkce, kterou můžeme získat pomocí jednorozměrného symetrického jádra K\ (K\ e £02)- Existují dva typy těchto jader: • součinové jádro Kp(x, y) — K\(x) ■ K\(y), • sféricky symetrické jádro Ks(x, y) — c^1Kí (\Jx2 + y2^, ck — J J K1 (\Jx2 + y2) dx dy. Příklad2.2. Epanečnikovo jádro, které je v jednorozměrném případě tvaru K(x) — x2)I[_i ^ (x), má následující dvourozměrné varianty Kp(x,y) = ^(1 - x2)(l - y2) pro -l, která zahrnuje diagonální matice, • třída T', která obsahuje tzv. plné matice. Rozdíly mezi jednotlivými maticemi jsou patrné z tabulky 5.1, kde jsou zobrazeny vrstevnice součinového Epanečnikova jádra v závislosti na třídě matic. Tabulka 5.1: Třídy vyhlazovacích matic S V T Budeme se zabývat jádrovými odhady s diagonální vyhlazovací maticí. Jádrový odhad s maticí třídy S dává ve všech směrech stejnou míru vyhlazení, což neponechává příliš mnoho prostoru pro zachycení variabilty dat. Na druhou stranu při použití matice třídy T je potřeba odhadnout větší počet parametrů, což znamená vyšší výpočetní náročnost. Konstrukce jádrového odhadu je analogická konstrukci jednorozměrného odhadu. Tedy v každém bodě [Xi, Yi] sestrojíme jádro Ku a odhad v bodě [x, y] je průměr n hodnot jader v tomto bodě - viz obrázek 5.4. 3 Statistické vlastnosti jádrových odhadů hustoty Stejně jako u jádrových odhadů jednorozměrných hustot můžeme kvalitu jádrového odhadu hustoty popsat lokálně pomocí střední kvadratické chyby: MSE f(x, y; H) = - * /)(*, V) ~ {Ku * f?{x, y)) + ((Ku * f)(x, y) - f(x, y)f, var bias 62 nebo globálně pomocí střední integrální kvadratické chyby MISE/(x,y;H) = J J MSE f (x, y; H) dx dy. Poznámka 3.1. Podobně jako v jednorozměrném případě se definuje konvoluce funkcí dvou proměnných. Nechf jsou dány funkce f ag, pro které platí f f f2 (x, y) dx dy < oo a f f g2 (x, y) dx dy < oo. Konvoluci f * g definujeme vztahem (f*g)(x,y)= // f{t,u)g{x -t,y -u)átáu. Optimální vyhlazovací matice minimalizuje MISE. Je zřejmé, že tyto optimální hodnoty vyhlazovacích parametrů není možné z MISE přímo vyjádřit. Stejně jako u odhadu jednorozměrných hustot se budeme zabývat asymptotickou střední integrální kvadratickou chybou AMISE. Věta 3.1. Předpokládejme, že funkce f, jádro K a matice vyhlazovacích parametrů2 H — diag(/if, h\) splňují následující předpoklady. (i) Nectí H — H„ je posloupnost matic takových, že (n \ H |) ~1 a prvky matice H konvergují k nule pro n —> oo. (íí) Dále nectí všechny druhé parciální derivace funkce / jsou ohraničené, spojité a íntegrovatelné se čtvercem. (ííí) Jádro K splňuje xK(x,y)dxdy — J J yK(x, y) dx dy — 0 x2K(x, y) dxdy — jjy2 K (x, y) dxdy — Í32(K). Pak platí kde MISE(H) = AMISE(H) + o(hj + h22) + o((/ii/i2n)"1), AMISE(H) = AMISE/(■,H) = + -^2{K)(h\V{fxx) + 2h2h2V(fxfy) + h42V(fyy)), (5.2) přičemž označení je ve shodě s předchozími kapitolami, tj. V(g) — JJ g2(x, y) dx dy. 'hl o 2Užíváme zde zkráceného zápisu: diag(fti, /12) : 0 h22 63 Důkaz věty o tvaru AMISE je založen na Taylorově rozvoji funkce f(x, y) a lze jej nalézt např. v knize [14]. Hodnoty parametrů h\, h^, pro které AMISE(/ii, h^) nabývá minimální hodnoty, jsou dány vztahy: h = (_V^(fyy)V(K) \1/6 'OPt \nft(K)VV*(fxx)[V(fxfy) + y/V(fxx)V(fyy)\ J ' ( V3/4(fxx)V(K) h2'°Pt = \nPl{K)VV^fyy)[V{fxfy) + y/VUxx)VUyy)\ ) ' (5'4) Z těchto vztahů plyne, že množina přípustných vyhlazovacích parametrů je tvaru ain-1/6 < hi < birT1^, a^n-1/6 < hi < 1/r6 pro vhodné konstanty 0 < a\ < b\ < oo, 0 < a,2 < &2 < oo. 4 Volba jádra Podobně jako u odhadu jednorozměrné hustoty není volba jádra podstatná. Je vhodné zvolit součinový tvar optimálního jádra. Tím zajistíme jistou hladkost výsledného odhadu a navíc výpočty s využitím součinových jader jsou jednodušší. Poznámka 4.1. V literatuře se také využívá Gaussovo jádro K(x, y) — (2tt)^1 e-^' +v ^2, které se zdá být výhodnějším při studiu asymptotických vlastností jádrového odhadu. Na druhou stranu má nevýhodu, že jeho nosičem je celá reálná osa, což způsobuje „nedokonalost" při odhadech hustot s omezeným definičním oborem. 5 Volba vyhlazovacího parametru 5.1 Metoda referenční hustoty Předpokládejme, že náhodný výběr ([Xi, Y\],..., [Xn,Yn]) pochází z normálního rozdělení s hustotou f(x,y) = --e "i -i Z7T(7l(72 a jádro K je dvourozměrnou standardizovanou normální hustotou, tj. K(x,y) = — e ^ 2 . Pak podle metody referenční hustoty lze získat tyto odhady vyhlazovacích parametrů /ij.ref = CTjn-1/6, 2 = 1,2. Tento vztah je také znám jako Scottovo pravidlo ([11]). Příklad 5.1. Pro simulovaná data z příkladu 2.1 vychází matice vyhlazovacích parametrů podle metody referenční hustoty s využitím součinového Epanečnikova jádra takto: /0,3595 0 Href — \ 0 0,9972^ Na obrázku 5.5 je vykreslen odhad hustoty s touto maticí a porovnání odhadu se skutečnou hustotou. Poznámka 5.1. V toolboxu, který doplňuje tato skripta, se uvádí druhá mocnina matice vyhlazovacích parametrů, tedy H2. Z obrázku 5.5 je patrné, že jádrový odhad je podhlazený, zejména ve směru osy x. Je to způsobeno odhadem směrodatné odchylky, která je základem metody referenční hustoty. 64 -2 0 2 -2 0 2 (a) Skutečná hustota (b) Data (c) Odhad s Href (d) Porovnání odhadu (—) se skutečnou hustotou (—) Obrázek 5.5: Jádrový odhad dvourozměrné hustoty - referenční hustota 5.2 Metoda křížového ověřování Metoda křížového ověřování je založena na odhadu hustoty v bodě [Xi, li] s vynecháním tohoto pozorování. Funkci metody křížového ověřování CV můžeme zapsat ve tvaru CV(H)= //(/(x,y,H))2dxdy--^/_í(Xí,yí,H). n ■ i=l kde 1 " f-i(Xi, li, H) = —- Ka(Xi - Xj,Yí - Yj) i=i Někdy se metoda C V nazývá nevychýlená metoda křížového ověřování (unbiased cross-validation), důvodem je jednoduchý vztah mezi CV a MISE, který uvádí následující věta. Věta 5.1. Funkce CV je nevychýleným odhadem MISE, tj. platí ECV(H) = MISE(H) - J J f (x, y) dxdy. Důkaz. Vypočtěme střední hodnotu CV: £CV(H) = £ (i f2(x,y,H)dxdy --Y.Ef^X^n) J J n i=l = E J J f2(x,y,n) dxdy -2EKU(X1 - X2,Y1 -Y2) = E /2(^,y,H)dxdy - 2 //// Ku(x - u, y - v)f(x, y)f(u, v) dxdy du dv 65 a úpravou MISE MISE(H) = E J J (f(x,y,H)-f(x,y)) dxdy = eJJ f2(x,y,H)dxdy-2E J f {x,y,K) f {x,y) dxdy + J J f {x,y) dxdy = E íl f2 (x, y,H) dxdy - 2 jjjj KH(x - u, y - v) f (x, y) f (u, v) dxdy dudv f (x, y) dx dy. Porovnaním upravených výrazů dostaneme tvrzení. □ Optimální matici vyhlazovacích parametrů vzhledem k metodě C V označíme H^v/ tj. Hcv — arg min CV(H). Hex> Příklad 5.2. Použijeme-li součinové Epanečnikovo jádro, pak pro simulovaná data z příkladu 2.1 dostanem matici vyhlazovacích parametrů určenou podle metody křížového ověřování v následujícím tvaru: '1,2055 0 0 0,3783, Na obrázku 5.6 je vykreslen odhad hustoty s touto maticí. H cv -202 (a) Skutečná hustota x** * «*„*»«-* •3<« * x x x * x x »x*»' -2 0 2 (b) Data (c) Odhad s Hc (d) Porovnání odhadu (—) se skutečnou hustotou (—) Obrázek 5.6: Jádrový odhad dvourozměrné hustoty - metoda křížového ověřování Odhad hustoty na obrázku 5.6 se jeví podhlazený, zejména ve směru osy y. Metoda křížového ověřování nedává dobré výsledky pro odhad hustoty vícerozměrných dat, částečně je to způsobeno problémy při minimalizaci funkce křížového ověřování. 6 Aplikace na reálná data Tento datový soubor pochází ze studie koncentrace lipidů v krevní plazmě, která vyšla v časopise Circulation v roce 1980. Výběrový soubor, který jsme převzali z [11] a s nímž zde pracujeme, obsahuje měření množství cholesterolu a triglyceridu v krevní plazmě u 320 pacientů, kteří si stěžovali na bolest v hrudníku. Data jsou shrnuta v tabulkách 6.9 a 6.10. 66 Obrázek 5.7: Vrstevnicové grafy odhadnutých hustot pro koncentraci lipidů - na ose x je vyneseno množství cholesterolu (v miligramech na 100 ml plazmy) a na ose y množství triglyceridu v krevní plazmě (mg/100 ml) Matice vyhlazovacích parametrů určené podle metody referenční hustoty a metody křížového ověřování jsou následující: (16,45 0 \ /42,39 0 \ \ 0 38,94y y 0 29,88y Na obrázku 5.7 jsou znázorněna data a vrstevnice jádrového odhadu s Epanečnikovým součinovým jádrem. 67 Shrnutí Odhad dvourozměrné hustoty pravděpodobnosti f(x, y) v bodě [x, y] je tvaru 1 n n' ' i=i Dva typy jader: • součinové jádro: Kp(x,y) — K\(x) ■ K\(y), sféricky symetrické jádro: Ks{x, y) — ck xK\{^Jx2 + y2), Ck — jj K\{^J x2 + y2) dx dy. Asymptotická střední integrální kvadratická chyba dvourozměrného jádrového odhadu AMISE(H) = + -^l{K)(h\V{fxx) + 2hjhlV(fxfy) + h22V{fyy)). Optimální vyhlazovací parametry vzhledem k AMISE \ 1/6 _V3/\fyy)V(K) 'OPt \nft(K)W*(fxx) [V(fXfy) + y/V(fxx)V(fyy)\ h\ opt — \np2(K)v-^(;xx)[V(;x;y) + ^v(Ux)v(hv)\ ) 1/6 _ /_V3/\fxx)V(K)_ l2'°Pt " \n^(K)V^(fyy)[V(fxfy) + ^V{fxx)V{fyy)\ Metody pro odhad optimálních hodnot matice vyhlazovacích parametrů H — diag(/ii, h2) • metoda referenční hustoty • metoda křížového ověřování HCv = argminCV(H), CV(H) = í í (f(x, y, H))2 dxdy - - V f^Xi, Yit H). i=i Dodatky a cvičení 1. Určete součinové a sféricky symetrické dvourozměrné jádro odvozené z kvartického jádra K{x) = ^{i-x2y. 2. Odvoďte vztahy (5.3) a (5.4) pro optimální vyhlazovací parametry. 3. Aplikujte metodu referenční hustoty a metodu křížového ověřování na simulovaná i reálná data. 68