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 Masarykova univerzita jaroslav hajekcenter Ivana Horová, Jan Koláček, Kamila Vopatová 1/68 Obsah O Úvod @ 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 © Simulovaná data © Reálná data O Reference Ivana Horová, Jan Koláček, Kamila Vopatová 2/68 > Úvod Obsah O Úvod © 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 © Simulovaná data 0 Reálná data © Reference 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 c/-rozměrný náhodný výběr X-|, ..., Xn z rozdělení s hustotou f definujeme jádrový odhad hustoty ?(x, H) = ±J2 K»(* ~ X') = lW/2 É K(H-V2(x ~ X,)), /'=1 /'=1 • K je c/-rozměrná jádrová funkce, pro kterou platí fRd K(x) dx = 1 • H je matice vyhlazovacích parametrů z množiny T symetrických pozitivně definitních matic typu d x d . x = (x1,...,xc/)rGMQř Ivana Horová, Jan Koláček, Kamila Vopatová 5/68 > Úvod Data - Jádro - Hustota Pro daná data Ivana Horová, Jan Koláček, Kamila Vopatová > Úvod Data - Jádro - Hustota vybereme jádro, např. Epanečnikovo x x x x x 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 optimalityje 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,..., h%), S : třída nejjednodušších matic typu H = /t2 • ld Pro třídy V a S lze psát odhad hustoty v jednodušším tvaru n /Ji xd-Xd V f(x,H) S[ ?(X'H) = J^E^((x-x')/^) ;=1 Ivana Horová, Jan Koláček, Kamila Vopatová 12/68 > Metody pro odhad H - Metody založené na MISE Obsah O Úvod @ 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 © Simulovaná data © Reálná data O Reference Ivana Horová, Jan Koláček, Kamila Vopatová 13/68 > Metody pro odhad H - Metody založené na MISE MISE Kvalitu odhadu f můžeme měřit např. pomocí střední integrální kvadratické chyby (MISE) [Wand & Jones, 1995] MISEf(-, H) f(x,H) - f(x)) dx varf(x,H)dx + / (bias/^x,H)) dx, kde varf(x,H)dx=- / K£(x-y)r(y)dy ■ KH(*-y)f(y) J (bias?(x,H))2dx = J(| /<„(x-y)f(y,H)dy-?(x,H))2dx Ivana Horová, Jan Koláček, Kamila Vopatová 14/68 Metody pro odhad H - Metody založené na MISE LSCV Nevychýlená metoda křížového ověřování odhaduje integrální kvadratickou chybu (ISE) [Sain et al, 1994] ISEf{-, H) = í (f{x, H) - f{x)Y dx = j[f{x, H)]2dx-2 j ?(x, H) ■ f(x) dx + j f2{x) dx LSCV(H) = J [?(x, H)]2 üx--nJ2 f_i(Xh H), kde ?_,-(x, H) = £jLi Kh{x - Xj) LSCV(H) = n'2 £ (KH * KH - 2K„)(X,- - Xy) + 2^1 KH(0) Ivana Horová, Jan Koláček, Kamila Vopatová 15/68 > Metody pro odhad H - Metody založené na AMISE 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 0 Reálná data 0 Reference 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 I ( f{x, H) - f{x)Y dx var f(x, H) dx + / (bias f(x, H)) dx, AMISE f (■ ,H) = -n[j K2H{x- y)r(y) dy AlVar + 1(1 K„(x-y)r(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 = (32{K)ld, kde Í32{K) = f xfK(x) dx nezávisí na / a / oo. • Všechny prvky matice 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) = AlVar + AlBias2 Vychýlení (Bias) lze s užitím vícerozměrné Taylorovy věty rozepsat (KH * f)(x) - f(x) = I KH(x - y)f(y) dy - f(x) = j K(z)f(x - H1/2z) dz - f(x) = f(x) + \f52{K)\r[HV2{x)] + o(trH) - f(x) Bias ~ l(32{K)\r[HV2{x)] Ivana Horová, Jan Koláček, Kamila Vopatová 19/68 Metody pro odhad H - Metody založené na AMISE AMISE Podobně můžeme rozepsat i rozptyl (Var) 1 var f(x, H) = -[J K^(x - y)ŕ(y) dy - (J KH(x - y)r(y) = 1|H|-V2 J K2(z)f(x-HV2z)óz = J[-\H\-^2V{K)f{x) + o{n-^\H\-^2) Var- l\H\-VzV(K)f(x) Tedy AMISE(H) = ^\H\-^2V{K) + \ß\{K) j\ŕ[HV2{x)] dx Ivana Horová, Jan Koláček, Kamila Vopatová 20/68 Metody pro odhad H - Metody založené na AMISE AMISE Dále ještě upravíme tvar vychýlení j tr2 [HVf(x)] dx = (vech H)7^ vech H. Operace vech (z angl. vector half): A={cd) => VQChA Pro d = 2 je matice vi/^ tvaru V>4,0 2^3,1 V>2,2 2-03,1 4-02,2 2^1)3 02,2 20-, )3 0o,4 kde [dk+ef(x)f. ^ = 7 afa^ř(x)dx' Ivana Horová, Jan Koláček, Kamila Vopatová 21/68 > Metody pro odhad H - Metody založené na AMISE AMISE AMISE(H) = -\H\-V2V(K) + l/3|(/<)(vech H)r^vechH 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(/7f, ...,H^)= diag(h2) a AMISE můžeme psát ve tvaru AMISE(H) = n^{K)h + }/?22(/<)(h2)r^(h2) kde vi/v obsahuje prvky V2e/+2e;> je jednotkový vektor s 1 na /-té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ßUK)h* j [£ d2 f {x) /=1 dx? dx /(D2) Jen v tomto případě lze vyjádřit optimální hodnotu vyhlazovacího parametru h 'AMISE V(K) nß2{K)l{D*) 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: AMISE(H) = + l/?2(K)2(>TÍV4,o + 2^1.V>2,2 + ^0,4), kde • V{K) = ///<2(x1,x2)dx1dx2 < 00 f32{K) = j j xfK{x,, x2) dx,dx2, / = 1,2 ' ^k,e = Il(^f)f(^,x2)dx,dx2, k,£ = 0,2,4, k + t = 4 Ivana Horová, Jan Koláček, Kamila Vopatová 24/68 Metody pro odhad H - Metody založené na AMISE Lemma Označme Hamise = arg m\nHeT,AMISE(H). Pak pro vyhlazovací matici H aMise s prvky '1 amise '2 amise 1/2 ,1/2x 1/2 ,1/2x 1/6 1/6 platí var f(x, Hamise) d*i dx2 = 2 / / (bias f(x, Hamise) ) dx: dx2 Ivana Horová, Jan Koláček, Kamila Vopatová 25/68 Metody pro odhad H - Metody založené na AMISE AMISE Odhad AMISE ÁMISE(H) = [[var?(x,H)dx1dx2 + f f (bías?(x,H))2 óx:óx2, kde Jj var?(x, H)óx:óx2 = ^JJ\H\-^2V{K)f{x, H)óx:óx2, bíis?(x, H))2 6x,6x2 = ff^ff KH(x - y)?(y, H)dyidy2 -?(x,H))2dx1dx2. Ivana Horová, Jan Koláček, Kamila Vopatová 26/68 > Metody pro odhad H - Metody založené na AMISE ÁMÍSE - optimální H Úpravou předešlých rovnic získáme //var?(x;H)dx,dx2 = ^. JI (bías?(x; H))2 óx:óx2 ^ n n = -fjžYl ^2(Kh *KH*Kh*Kh i=1 y=1 - 2KH *KH*KH + KH* KH)(Xi - Xj), kde * značí operaci konvoluce. Ivana Horová, Jan Koláček, Kamila Vopatová 27/68 > Metody pro odhad H - Metody založené na AMISE ÁMÍSE - optimální H Označme n g{h,,hz)=YtKH*KH*KH*KH-2KH*KH*KH + KH*KH){Xi-Xj) Myšlenka této metody je založena na Lemmatu, tedy hledáme taková fy a h2, aby platilo IVar = 2 • IBias2 ^ = 2-^(fy,fy) /ify h2 n nV{K) = 2h,h2g{h^h2) Ivana Horová, Jan Koláček, Kamila Vopatová 28/68 Metody pro odhad H - Metody založené na AMISE ÁMÍSE - optimální H Ivana Horová, Jan Koláček, Kamila Vopatová 29/68 Metody pro odhad H - Metody založené na AMISE Optimální H -> M1 Jak najít další vztah mezi fy a fy? O Scottovo pravidlo: fy = ůjn-V6, (/' = 1,2) [Scott, 1992] fy = C fy , C = ^ m í 2fyfyg(fy,fy) = A7l/(/<) | fy = Cfy K odhadu a, můžeme použít výběrovou směrodatnou odchylku ^TT^É^-*'-)2' / = 1'2' nebo robustnější odhad n- 1 7=1 a/'/0fí--1^49-' ' = 1'2- Doporučuje se použít menší z odhadů: min(a,, číjor), / = 1,2. Ivana Horová, Jan Koláček, Kamila Vopatová 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 Optimální H -> M2 © Podle vztahu pro optimálními hodnoty fy a fy (vzhledem k AMISE) platí fy f VHP V/4 O ^^04 = ^1^40, / ífd2f\\V n-2\-(d2K» d2K»\(Y Ví V>40 9^ dx^V&^W-X,) f/ 9x2 ^ J) SX,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 h^ýoA = ^1^40 plyne í 2h,h2g(fh,h2) = nV(K) M2 ^-2E^@«§)(X,-X;.) = ( =^-2E^@«^)(x,-x,) 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 r dk+íf(yč\ ^k/ = J ^^•f(X)dX' M = 0'2'4' k + £ = 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 BCV O ^w > ■'(x'H) dx" n J 'H(x) n n *«=»-TE(S*'í")(Xi-x') i=1 y=1 ° 1 a 2 /0*+V(x)\ 1 f9*4ř(Xi,H) 9x^9x1 / n ^ 9x^9x1 /=1 y=1 aX1 aX2 kde ?_/(x, A7) = (n — 1 )~1 £jLi,#/ _ [Duon9 & Hazelton, 2004] Ivana Horová, Jan Koláček, Kamila Vopatová > Metody pro odhad H - Metody založené na AMISE Plug-in metoda = F/a*+žf(xh 1 A dk+rf(xh G) dX+'KaiXi-Xj) i=1 y=1 kde předpokládáme G = g2l Qamse 2dk+lK(0) dxfdxt 4+k+i [Wand & Jones, 1994] Ivana Horová, Jan Koláček, Kamila Vopatová 37/68 > Jiné metody pro odhad H Obsah O Úvod @ 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 © Simulovaná data © Reálná data O Reference Ivana Horová, Jan Koláček, Kamila Vopatová 38/68 > Jiné metody pro odhad H PLCV Metoda křížového ověřování pomocí pseudověrohodnostní funkce [Cao, 1993] n L(H) = l\f_i(Xi,H) max £(H) = In L(H) n n = -n\n(n - 1) - g Indet(H) + Eln [E K("~1/2(X/ - X,)) i=1 y=1 Ivana Horová, Jan Koláček, Kamila Vopatová 39/68 Jiné metody pro odhad H Metoda referenční hustoty Za předpokladu, že f je c/-rozměrná normální hustota, pak optimální hodnota matice vyhlazovacích parametrů je 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 Hr ((d + 8)^/*V(K)\^ , \ i6A?r(^)(d + 2) ) pro d = 2 H.HT = (<*s*vm)i.t V 96/1 J přičemž Z značí odhad kovarianční matice. Pro H = diag{/71,A72} [Sain et al, 1994] Ivana Horová, Jan Koláček, Kamila Vopatová 41/68 > Simulovaná data Obsah O Úvod @ 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 © Simulovaná data © Reálná data O Reference 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{x,, x2) = A(1 _ X12)(1 - xf), X1, x2 g [-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 I n= 100, R= 100 Ivana Horová, Jan Koláček, Kamila Vopatová 44/68 > Simulovaná data Ivana Horová, Jan Koláček, Kamila Vopatová 45/68 > Simulovaná data Normální hustota I - rekonstrukce LSCV (o) M1 (x) M2 ( ) > Simulovaná data Normální hustota II n= 100, R= 100 X ~0.25A/2(0,0; 1,1,0) + 0.75A/2(4,3;4,3, 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 (+) Ivana Horová, Jan Koláček, Kamila Vopatová 49/68 > Simulovaná data Normální hustota III n = 100, R = 100 Ivana Horová, Jan Koláček, Kamila Vopatová 50/68 > Simulovaná data Normální hustota III - H 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 ~r(5) • ř(5) Ivana Horová, Jan Koláček, Kamila Vopatová 53/68 > Simulovaná data Ivana Horová, Jan Koláček, Kamila Vopatová 54/68 > Simulovaná data Studentovo rozdělení - rekonstrukce LSCV (o)__M1 (x)__M2 (+) Ivana Horová, Jan Koláček, Kamila Vopatová 55/68 > Simulovaná data Weibullovo rozdělení n = 60,R= 100 X ~IMb(2,2) • Wb{2,3) Ivana Horová, Jan Koláček, Kamila Vopatová 56/68 > Simulovaná data Weibullovo rozdělení - H O ^LSCV Ivana Horová, Jan Koláček, Kamila Vopatová 57/68 > Simulovaná data Weibullovo rozdělení - rekonstrukce LSCV (o)__M1 (x)__M2 (+) ^^^^^^^^^^^^^ 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 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) f(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 Ú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 0 Reference 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, 500 - • - cholesterol [mg/100 ml], 400 - • X2 - triglycerid [mg/100 ml]. 300 - gl-1-1-1-1-1-1-1-1 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) th =42.31, h2 = 31.86 M2( ) th = 14.99, h2 = 25.58 0 .' 50 100 160 200 250 300 350 400 460 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 Ivana Horová, Jan Koláček, Kamila Vopatová 64/68 > Reference Obsah O Úvod @ 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 © Simulovaná data © Reálná data O Reference Ivana Horová, Jan Koláček, Kamila Vopatová 65/68 > Reference Reference t> Cao R., Cuevas A., Gonzalez Manteiga W. A comparative study of several smoothing methods on density estimation Comput. Statist. Data Anal. 17, 1994. t> Duong T., Hazelton M.L. Convergence rates for unconstrained bandwidth matrix selectors in multivariate kernel density estimation. J. Multivariate Anal. 93, 2005. t> Duong T., Hazelton M.L. Cross-validation Bandwidth matrices for Multivariate Kernel Density Estimation. Scand. J. Statist. 31, 2004. t> Horovä I. era/. Bandwidth choice for kernel density estimates. Proceedings IASC. Yokohama: IASC, 2008. t> Sain S.R., Baggerly K.A. and Scott D.W. Cross-Validation of Multivariate Densities. J. Amer. Statist. Assoc. 89, 807-817, 1994. t> Scott D.W. Multivariate Density Estimation: Theory, Practice, and Visualization New York: John Wiley & Sons, 1992. t> Wand M.R, Jones M.C. Multivariate Plug-in Bandwidth Selection. Comput. Statist. 9, 1994. t> 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ě. o + o o 'b o o °" 5 o ŕ h° ° o "ä9 i° o, 0 a c ^ + o -tg o o o + ° oA co ° o c ° 0 H A ° 3 o o o o o o c + C3 O A O O o A 3 O c -k) O A °A O o 3. C Ivana Horová, Jan Koláček, Kamila Vopatová 67/68 Pouštní rostliny - dead Ivana Horová, Jan Koláček, Kamila Vopatová 68/68