> Volba vyhlazovací matice pro jádrové odhady vícerozměrných hustot Ivana Horová, Jan Koláček, Kamila Vopatová Oddělení aplikované matematiky ■ ■ ■ Přírodovědecká fakulta Hfl Masarykova univerzita jaroslav hájek center Ivana Horová, Jan Koláček, Kamila Vopatová 1/68 Obsah O Úvod 0 Metody pro odhad H - Metody založené na MISE O Metody pro odhad H - Metody založené na AMISE O Jiné metody pro odhad H 0 Simulovaná data © Reálná data O Reference Ivana Horová, Jan Koláček, Kamila Vopatová 2/68 > Úvod Obsah O Úvod ^^K^^^k WL H m III ■ ■ . \/ f Ivana Horová, Jan Koláček, Kamila Vopatová 3/68 > Úvod Motivace Ivana Horová, Jan Koláček, Kamila Vopatová 4/68 > Úvod Jádrový odhad vícerozměrné hustoty Pro d-rozměrný náhodný výběr X1, ..., Xn z rozdělení s hustotou f definujeme jádrový odhad hustoty ?(x, H) = 1 JT KH(x - X,) = 1|H|-V2 J2 K(H~y\x - X/)), /=1 /=1 • /C je d-rozměrná jádrová funkce, pro kterou platí JRd K(x)áx = 1 • H je matice vyhlazovacích parametrů z množiny J?7 symetrických pozitivně definitních matic typu d x d • x = (x1,...,xd)7GlRaf Ivana Horová, Jan Koláček, Kamila Vopatová 5/68 > Úvod Data - Jádro - Hustota Pro daná data Ivana Horová, Jan Koláček, Kamila Vopatová 6/68 > Úvod Data - Jádro - Hustota vybereme jádro, např. Epanečnikovo (mm 1 1 1 1 1 1 - X - X X X X X X X X 1 1 1 1 X 1 1 1 1 1 Ivana Horová, Jan Koláček, Kamila Vopatová 7/68 > Úvod Data - Jádro - Hustota vyčíslíme v každém bodě Ivana Horová, Jan Koláček, Kamila Vopatová 8/68 > Úvod Data - Jádro - Hustota a získáme rekonstruovanou hustotu. Ivana Horová, Jan Koláček, Kamila Vopatová 9/68 > Úvod Jádrový odhad vícerozměrné hustoty Typ jádra - z jednorozměrného jádra lze vytvořit dvě různá vícerozměrná jádra: • součinové jádro • sféricky symetrické jádro Ivana Horová, Jan Koláček, Kamila Vopatová 10/68 > Úvod Výběr mnohorozměrného jádra Vybrat si součinové jádro Kp nebo sféricky symetrické Ks? Měřítkem optimality je funkcionál Cd [Wand & Jones, 1995] Cd{K) = (V(K)4/32(K)2dy/{d+4) - min Řešením této úlohy jsou optimální jádra. Mezi součinovými jádry nejlépe vychází Epanečnikovo součinové jádro. Ivana Horová, Jan Koláček, Kamila Vopatová 11/68 > Úvod Matice vyhlazovacích parametrů H • nejdůležitější složka při jádrovém vyhlazování • má vliv na orientaci jádra a jeho „šířku" • 3 třídy vyhlazovacích matic T : obecná třída symetrických pozitivně definitních matic s \d(d + 1) nezávislými prvky, V : diagonální matice tvaru H = diag(/?f,..., h2d), S : třída nejednodušších matic typu H = h2 ld Pro třídy V a S lze psát odhad hustoty v jednodušším tvaru S: ř(x'W> = 7Sfí>((x-xi)/") /=1 Ivana Horová, Jan Koláček, Kamila Vopatová 12/68 > Metody pro odhad H - Metody založené na MISE Obsah O Metody pro odhad H - Metody založené na MISE Ivana Horová, Jan Koláček, Kamila Vopatová 13/68 > Metody pro odhad H - Metody založené na MISE MISE Kvalitu odhadu f můžeme merit napr. pomocí strední integrální kvadratické chyby (MISE) [Wand & Jones, 1995] MISEf(- ,H) = E I (7(x, H) - f(x)j dx = ívar?(x,H)dx+ f fbias?(x,H)]2 dx. kde 1 var f (x, H) dx = - Kg(x-y)/(y)dy- ( KH(x-y)f(y) 2i bias?(x,H))2dx= í (/" KH(x-y)/(y,H)dy-?(x,H))2dx, Ivana Horová, Jan Koláček, Kamila Vopatová 14/68 > Metody pro odhad H - Metody založené na MISE Nevychýlená metoda křížového ověřování odhaduje integrální kvadratickou chybu (ISE) [Sain et al, 1994] ISEf(-,H)= í (?(x, tf)-/(x))2dx = /[?(x,H)]2dx-2 f ?(x, H) ■ f(x) dx + J ř2(x)dx n LSCV(H) = í[?(x, H)f dx - ^ ?_/(X/, H) J ;'=1 kde L,(x, H) = T^r) EyLi KH(x - Xy) LSC= n-2 £ (ACH * KH - 2KH){X, - Xy) + 2n~1 KH(0) Ivana Horová, Jan Koláček, Kamila Vopatová 15/68 > Metody pro odhad H - Metody založené na AMISE Obsah O Metody pro odhad H - Metody založené na AMISE Ivana Horová, Jan Koláček, Kamila Vopatová 16/68 > Metody pro odhad H - Metody založené na AMISE Od MISE k AMISE MISEf{- ,H) = E J (?(x, H) - r(x))2 dx = / var?(x,H)dx + f (bias?(x, H))2 dx. AMISE f (■, H) = 1 [ / Kg(x - y)/(y) dy -v-' AlVar + l{ I K„(x-y)/(y,H)dy-?(x,H))2dx AlBias2 Ivana Horová, Jan Koláček, Kamila Vopatová 17/68 > Metody pro odhad H - Metody založené na AMISE AMISE Předpoklady pro přechod od MISE k asymptotickému tvaru střední integrální kvadratické chyby (AMISE): • K je ohraničená jádrová funkce s kompaktním nosičem, pro niž platí J K(x) dx = 1, J xK(x) dx = 0, J xxTK(x) dx = 02{K)ld, kde @2(K) = / xfK(x) dx nezávisí na / a ld je jednotková matice řádu d. • h = hn\e posloupnost vyhlazovacích matic takových, že n-i |/-/|-V2 a všechny prvky h se blíží k nule pro n —► oo. • Všechny prvky matice Vf druhých derivací hustoty f jsou po částech spojité a integrovatelné se čtvercem. Ivana Horová, Jan Koláček, Kamila Vopatová 18/68 > Metody pro odhad H - Metody založené na AMISE AMISE Asymptotický tvar střední integrální kvadratické chyby [Wand & Jones, 1995] AMISE(H) = AI Var + AlBias2 Vychýlení (Bias) lze s užitím vícerozměrné Taylorovy věty rozepsat (K„ * 0(x) - m = J KH(x - y)ř(y) dy - f(x) = í K{z)f(x - H1/2z) dz - /(x) = f(x) + l/3z(K)tr[HVf(x)] + o(trH) - f(x) Bias ~ \p2{K)Vc[HV%x)\ Ivana Horová, Jan Koláček, Kamila Vopatová 19/68 > Metody pro odhad H - Metody založené na AMISE Podobne můžeme rozepsat i rozptyl (Var) var?(x, H) = \ [J K2 (x - y)f(y)dy - (| KH(x-y)f(y) 2i l|H|-1/2 ľK2(z)f(x-H1/2z)dz 1 n H\-112 V{K)f(x) + o(A7"11 H\-112) Var ~ ±\H\-V2V(K)f(x) Tedy AMISE(H) = -\H\-^2V{K) + \ß\{K) ítr2[HD2(x)] dx Ivana Horová, Jan Koláček, Kamila Vopatová 20/68 > Metody pro odhad H - Metody založené na AMISE Dále ještě upravíme tvar vychýlení Jtr2 [HVf(x)] dx = (vech H)TWTvech H. Operace vech (z angl. vector half): Pro d = 2 je matice tvaru ^4,0 2^3;1 ^2,2 2^3,1 4^2,2 2^1,3 ^2,2 2^1 |3 ^o,4 kde [dk+ef(x),f , . ^=y^fř(x)dx Ivana Horová, Jan Koláček, Kamila Vopatová 21/68 > Metody pro odhad H - Metody založené na AMISE AMISE AMISE (H) = ^-\H\-^2V(K) + \ßl(K)(\/ec\\ H)TtyTvechH I přes tato zjednodušení nelze vyjádřit optimální H vzhledem k AMISE a musí se počítat numericky. [Wand & Jones, 1995] V : Je-li matice H diagonální, pak lze psát H = diag(A?f, ..., /72) = diag(h2) a AMISE můžeme psát ve tvaru AMISE(H) = n^{K)hd + l/322(/C)(h2)^^(h2) kde i|/£> obsahuje prvky V;2e/+2ey5 je jednotkový vektor s 1 na Mém místě. Ivana Horová, Jan Koláček, Kamila Vopatová 22/68 > Metody pro odhad H - Metody založené na AMISE AMISE S : V případě H = h2 • ld dostaneme AMISE(H) = ^ + lß22(K)h4 J d E ;=1 d2f(x)H2 dx Jen v tomto případe lze vyjádřit optimální hodnotu vyhlazovacího parametru h hAMISE d V(K) nß2(K)l(D2)l 1/(cř+4) Ivana Horová, Jan Koláček, Kamila Vopatová 23/68 > Metody pro odhad H - Metody založené na AMISE AMISE-2D Zaměříme se na odhady dvourozměrné hustoty s diagonální vyhlazovací maticí H. Pak asymptotický tvar střední integrální kvadratické chyby (AMISE) lze jednoduše vyjádřit takto: kde • V{K) = J J K\xA, x2) dxi dx2 < oo • ^k,e = ll(^0l)f(^,x2)dx^x2, M = 0,2,4, k + £ = 4 Ivana Horová, Jan Koláček, Kamila Vopatová 24/68 > Metody pro odhad H - Metody založené na AMISE Lemma Označme Hamise = arg minW€£>/WW/SE(/-/). Pak pro vyhlazovací matici HAMISE s prvky platí h 1 amise 4a v(K) 1 1/6 ß2(Kf'^(^2,2 + %A^4,o)n 1/2-1/2- h 2AMISE <7o V(K) 1 1/6 var f (x, HAMISE) d*i dx2 = 2 / / (bias f (x, HAMISE)) dxi dx2 Ivana Horová, Jan Koláček, Kamila Vopatová 25/68 > Metody pro odhad H - Metody založené na AMISE Odhad AMISE ÁMISE(H) = ffvar?(x, H)dx1dx2 + ff fbías?(x, H))2 dx1dx2, kde v%?(x,H)d*ld^ 'bías?(x, H)f áx,áx2 = f f( f f KH(x - y)?(y, H) 6y,dy2 - f (x, H))2 dxi dx2 Ivana Horová, Jan Koláček, Kamila Vopatová 26/68 > Metody pro odhad H - Metody založené na AMISE Úpravou předešlých rovnic získáme __ V(K) varř(x;H)dx1dx2 = 2 bias f(x; H)) dx-|dx2 ^ n n H * K/-i * K/-i * K/-i /=1 y=i - 2KH * K h * Ku + Ku * KH)(Xi - Xy), kde * značí operaci konvoluce. Ivana Horová, Jan Koláček, Kamila Vopatová 27/68 > Metody pro odhad H - Metody založené na AMISE ÁÍVliSE - optimální H Označme n 0(/7i ,h2)=Y,KH*KH*KH*KH-2KH *KH*KH + KH* KH)(X, - Xy) Myšlenka této metody je založena na Lemmatu, tedy hledáme taková /Ji a h2, aby platilo IVar = 2 ■ IBias2 V(K) n 1 . , nh-\ h2 n nV(K) = 2kh2g(huh2) Ivana Horová, Jan Koláček, Kamila Vopatová 28/68 > Metody pro odhad H - Metody založené na AMISE ÁMISE - optimální H Ivana Horová, Jan Koláček, Kamila Vopatová 29/68 > Metody pro odhad H - Metody založené na AMISE Jak najít další vztah mezi h-\ a /72? O Scottovo pravidlo: h, = ovrr1/6, (/' = 1,2) [Scott, 1992] h2 = C/7!, C = 11 ✓v M1 2h,h2g(h,1h2) = nV(K) h2 = C/7i K odhadu a, můžeme použít výběrovou směrodatnou odchylku /=i nebo robustnejší odhad aiJQR Xj[3n/4] ~~ ^/[n/4] 1.349 7 = 1,2. Doporučuje se použít menší z odhadů: min(a/, S7,/cv?), / = 1,2. Ivana Horová, Jan Koláček, Kamila Vopatová 30/68 > Metody pro odhad H - Metody založené na AMISE Řešení metody M1 Ivana Horová, Jan Koláček, Kamila Vopatová 31/68 > Metody pro odhad H - Metody založené na AMISE O Podle vztahu pro optimálními hodnoty fy a h2 (vzhledem k AMISE) platí Tr=\ — ) "2^04 = «1^40 '2 _ / ^40 N '1 "i ^04, d2f^ 2 ) dx d2f^ 2 ) dx n d*KH d*KH 2 n , ^2 i^.. ^2 Ivana Horová, Jan Koláček, Kamila Vopatová 32/68 > Metody pro odhad H - Metody založené na AMISE Optimální H -> M2 Pak ze vztahu hfeýoA = ri\i>Ao plyne /,y=1 * ^ /,y=1 1 1 2hih2g(huh2) = nV(K) ^2E&m(^*^)(x,-x,) H (X/ - Xy) Ivana Horová, Jan Koláček, Kamila Vopatová 33/68 > Metody pro odhad H - Metody založené na AMISE Řešení metody M2 Ivana Horová, Jan Koláček, Kamila Vopatová 34/68 > Metody pro odhad H - Metody založené na AMISE BCV a Plug-in metoda Cílem těchto metod je odhadnout funkci ^ = /^^"ř(X)dX' M = 0<2'4' k + i = 4 • Vychýlená metoda křížového ověřování (BCV) - 2 typy, • plug-in metoda. Ivana Horová, Jan Koláček, Kamila Vopatová 35/68 > Metody pro odhad H - Metody založené na AMISE -. „ • f(x, H)dx — dx(dx£ v ' n dk+eKH{x) ~dxjdxf KH(x)dx n n -2 i=i j= 2 * kh) (X 2 1 y=1 1 * -Xy) ' V dxkdxí, ) ~ n ^ dxfdxl ;'=1 -1 ö*+'/CH(X/ - Xy) /=1 y=i dxfdx. £ 2 kde /L,(x, H) = {n - 1 )"1 E"=i >yy/ Kh(X, - Xy). [Duong & Hazelton, 2004] Ivana Horová, Jan Koláček, Kamila Vopatová 36/68 > Metody pro odhad H - Metody založené na AMISE Plug-in metoda fdk+if(x)\ 1 A 9*+ Jiné metody pro odhad H Obsah O Jiné metody pro odhad H Ivana Horová, Jan Koláček, Kamila Vopatová 38/68 > Jiné metody pro odhad H PLCV Metoda krížového overovaní pomocí pseudoverohodnostní funkce [Cao, 1993] n L(H) = JI f-i(Xh H) —> max ;=1 £(H) = In L(H) n n _n|n(n - 1) - £ In det(H) + Eln ÍE K(AT1 /2(X,- - X,)) /=1 y'=i Ivana Horová, Jan Koláček, Kamila Vopatová 39/68 > Jiné metody pro odhad H Metoda referenční hustoty Za předpokladu, že f je d-rozměrná normální hustota, pak optimální hodnota matice vyhlazovacích parametrů je Hw, = (_l_)^>žn-2/(-+4> Z je odhad kovarianční matice. [Wand & Jones, 1995] Ivana Horová, Jan Koláček, Kamila Vopatová 40/68 > Jiné metody pro odhad H Metoda maximálního vyhlazení Hms je řešením rovnice 1 16/?r(^)(d + 2) 1 pro d = 2 1 96n ) přičemž Z značí odhad kovariancní matice. Pro H = diagjfy, h2} 1 ft/ 1 96n J / = 1'2 [Sain et al, 1994] Ivana Horová, Jan Koláček, Kamila Vopatová 41/68 > Simulovaná data Obsah 0 Simulovaná data Ivana Horová, Jan Koláček, Kamila Vopatová 42/68 > Simulovaná data Simulace Aplikace na simulovaných datech - porovnání metod M1 a M2 s metodou LSCV. • Epanečnikovo součinové jádro K{xux2) = ^(1 -^2)(1 -xf), xux2e [-1,1], • n - počet pozorování v náhodném výběru, • R - počet opakování. Ivana Horová, Jan Koláček, Kamila Vopatová 43/68 > Simulovaná data Normální hustota n= 100, R= 100 X~Ak(0,0; 1,1,0) Ivana Horová, Jan Koláček, Kamila Vopatová 44/68 > Simulovaná data Normální hustota I - H Ivana Horová, Jan Koláček, Kamila Vopatová 45/68 > Simulovaná data Normální hustota I - rekonstrukce LSCV(o) M1 (x) M2(+) -1-1-1-1-1 i-1-1-1-1-1-1 i-1-1-1- Ivana Horová, Jan Koláček, Kamila Vopatová 46/68 > Simulovaná data Normální hustota II n= 100, R= 100 X~0.25Afe(0,0; 1,1,0) + 0.75A/2(4,3;4,3,0) Ivana Horová, Jan Koláček, Kamila Vopatová 47/68 > Simulovaná data Ivana Horová, Jan Koláček, Kamila Vopatová 48/68 > Simulovaná data Normální hustota II - rekonstrukce LSCV(o) M1 (x) M2(+) -1-1-1-1-1 i-1-1-1-1-1-1-1 i-1-1-1-1— Ivana Horová, Jan Koláček, Kamila Vopatová 49/68 > Simulovaná data Normální hustota III n= 100, R= 100 v 3AÍ , , n 9 49 63 x X —-A/oí-1 0- — - -) 7 21 ' '25'100'250; 3A(/, 2 9 49 + 7A/2(1'7!:25'1ÖÖ'0) 1 Al 2 9 49 oX + 7A/2(1'-7!:25'TÖO'0) Ivana Horová, Jan Koláček, Kamila Vopatová 50/68 > Simulovaná data Ivana Horová, Jan Koláček, Kamila Vopatová 51/68 > Simulovaná data Normální hustota III - rekonstrukce LSCV(o) M1 (x) M2(+) Ivana Horová, Jan Koláček, Kamila Vopatová 52/68 > Simulovaná data Studentovo rozdělení n = 60, R= 100 X ~ř(5) • ř(5) Ivana Horová, Jan Koláček, Kamila Vopatová 53/68 > Simulovaná data Studentovo rozdělení - H Ivana Horová, Jan Koláček, Kamila Vopatová 54/68 > Simulovaná data Studentovo rozdělení - rekonstrukce LSCV(o) M1 (x) M2(+) -1-1-1-1-1 i-x-1-1-1-1-1-1 i-x-1-1-1- Ivana Horová, Jan Koláček, Kamila Vopatová 55/68 > Simulovaná data Weibullovo rozdělení n = 60, R= 100 X - Wb(2,2) • Wb{2,3) Ivana Horová, Jan Koláček, Kamila Vopatová 56/68 > Simulovaná data Weibullovo rozdělení - H Ivana Horová, Jan Koláček, Kamila Vopatová 57/68 > Simulovaná data Weibullovo rozdělení - rekonstrukce LSCV(o) M1 (x) M2(+) -1-1-1-1-1-1 i-1-1-1-1-1-1-1-1 i-1-1-1-1- Ivana Horová, Jan Koláček, Kamila Vopatová 58/68 > Simulovaná data Průměr relativních chyb AMISE Relativní chyba AMISE AMISE(Hmetoda) - AMISE(HAMISE)\ AMISE(HAMISE) data 1 1 LSCV M1 M2 Norm 1 0.5046 0.4486 0.0841 Norm II 0.5392 0.7227 0.7200 Norm III 0.4548 0.9232 0.6392 Student 0.5407 0.4052 0.3063 Weibull 0.8424 0.2987 0.1619 Ivana Horová, Jan Koláček, Kamila Vopatová 59/68 Simulovaná data Průměry lAE Integrální absolutní chyba (IAE) IAE(H) = / ?(x, H) - f(x) dx data LSCV M1 M2 Norm 1 0.3213 0.3016 0.2687 Norm II 0.3619 0.3468 0.3313 Norm III 0.4041 0.3828 0.3700 Student 0.4098 0.3378 0.3313 Weibull 0.3895 0.2942 0.2920 Ivana Horová, Jan Koláček, Kamila Vopatová 60/68 > Reálná data Obsah O Reálná data Ivana Horová, Jan Koláček, Kamila Vopatová 61/68 > Reálná data Lipidy v krvi Koncentrace lipidů v krevní plazmě pacientů. 320 pacientů, u kterých bylo zjištěno zúžení cév, X1 - cholesterol [mg/100 ml] X2 - triglycerid [mg/100 ml]. 600 500 400 300 200 100 n i r * * ♦ ****** ** ** ♦+ *** V .>* y*\ J_L 0 50 100 150 200 250 300 350 400 450 Ivana Horová, Jan Koláček, Kamila Vopatová 62/68 > Reálná data Lipidy v krvi - rekonstrukce LSCV (o) M2 (+) fy = 42.31, h2 = 31.86 fy = 14.99, h2 = 25.58 50 1GG 150 200 250 300 350 400 450 50 100 150 200 250 300 350 400 450 Ivana Horová, Jan Koláček, Kamila Vopatová 63/68 > Reálná data Lipidy v krvi - rekonstrukce LSCV M2 áf Ivana Horová, Jan Koláček, Kamila Vopatová 64/68 > Reference Obsah O Reference Ivana Horová, Jan Koláček, Kamila Vopatová 65/68 > Reference Reference > Cao R., Cuevas A., Gonzalez Manteiga W. A comparative study of several smoothing methods on density estimation Comput. Statist. Data Anal. 17, 1994. > Duong T., Hazelton M.L. Convergence rates for unconstrained bandwidth matrix selectors in multivariate kernel density estimation. J. Multivariate Anal. 93, 2005. > Duong T., Hazelton M.L. Cross-validation Bandwidth matrices for Multivariate Kernel Density Estimation. Scand. J. Statist. 31, 2004. > Horovä I. et ai Bandwidth choice for kernel density estimates. Proceedings IASC. Yokohama: IASC, 2008. > Sain S.R., Baggerly K.A. and Scott D.W. Cross-Validation of Multivariate Densities. J. Amer. Statist. Assoc. 89, 807-817, 1994. > Scott D.W. Multivariate Density Estimation: Theory, Practice, and Visualization New York: John Wiley & Sons, 1992. > Wand M.R, Jones M.C. Multivariate Plug-in Bandwidth Selection. Comput. Statist. 9, 1994. > Wand, M.P and Jones, M.C. Kernel Smoothing. London: Chapman and Hall, 1995. Ivana Horová, Jan Koláček, Kamila Vopatová 66/68 > Pouštní rostliny Měření dlouhověkosti u rostlin vybrané lokality Mohavské pouště, ^7 h A e ° o 8 QA° o i-sp-!-o—r ® O O + o A o <& o o 8 o o o o 0 o o o + o o o o o o A o ° o o o co + 8 o + o ° 0 o 00 0 °° A o o + o oA co 00 o o o + Oo° & (X- A O O o ^5 o Oj o * o o o o 'AA ° o + Q o A o o o o o o o o + O + o + o o o o o A o o o o " °S o A co + o o o o o a o o + Q3 OA O ^o o o o o + o o o o o o o o o o o o o + o o o o C£A HO o o o o + o J_r~\_I_y_L Ho O A oA o o o o o o. o _I_L Ivana Horová, Jan Koláček, Kamila Vopatová 67/68 > Pouštní rostliny - dead LSCV (o) M2 (+) fy = 6.01, h2 = 7.38 fy = 7.36, h2 = 7.96 Ivana Horová, Jan Koláček, Kamila Vopatová 68/68