Kapitola 3 Jádrové odhady hustoty Výstupy z výukové jednotky Student • bude znát tvar jádrových odhadů hustoty pravděpodobnosti. • bude schopen analyzovat statistické vlastnosti těchto odhadů. • se seznámí s metodami pro volbu vyhlazovacího parametru. • porozumí automatické proceduře pro simultánní volbu parametrů vyhlazování. • zvládne použití toolboxu v Matlabu a dokáže pro daný soubor dat zkonstruovat jádrový odhad hustoty a jejích derivací. 1 Motivace Hustota pravděpodobnosti je základním pojmem ve statistice [1, 2]. Odhadem hustoty rozumíme rekonstrukci hustoty z množiny naměřených dat. Tato rekonstrukce může poskytnout důležité informace o dané množině dat. Předpokládejme, že máme k dispozici nezávislé náhodné proměnné Xi,..., Xn, které mají tutéž spojitou hustotu /. Můžeme předpokládat, že neznámá hustota patří do třídy hustot, které závisejí na nějakých parametrech. Pro odhad hledané hustoty je tedy třeba odhadnout tyto parametry. Tento přístup se nazývá parametrický. My se zaměříme na neparametrický přístup, ve kterém se předpokládá pouze jistá hladkost odhadované hustoty (tj. dostatečný počet spojitých derivací). 2 Základní typy neparametrických odhadů Nejstarším neparametrickým odhadem hustoty je histogram [12, 11, 14]. Histogram zobrazuje relativní četnosti třídicích intervalů jako plochy obdélníků sestrojených nad těmito intervaly. Pak definujeme odhad hustoty četnosti f(x, h) — —(počet Xi ve stejném intervalu jako x), kde h značí šířku třídicích intervalů (obvykle se volí stejná šířka pro všechny intervaly). Nevýhody histogramu: • Histogram je citlivý na počet tříd a jejich šířku. 28 • Histogram je schodovitá funkce, ale přitom předpokládáme, že neznámá hustota je spojitá. Příklad 2.1. Mějme dán datový soubor generovaných ze směsi dvou normálních hustot N(0; 0,25) a N(2; 0,25), který má rozsah n = 100. 1 2 1 0-2)2 f (x) = 0,3^== e~fe +0,7 _ e ~ (Data jsou v tabulce 6.2.) Na obr. 3.1 je patrné, že histogram nevystihuje hustotu pravděpodobnosti dat. -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 (a) 7 intervalů o šířce h = 0,55 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 (b) 13 intervalů o šířce h = 0,29 Obrázek 3.1: Histogramy s různými počty třídicích intervalů Výše uvedené problémy lze odstranit použitím jádrových odhadů. Jádrový odhad hustoty / v bodě x G M je definovaný vztahem (Rosenblatt (1956), Parzen (1962)) 1 n Y 1 n (3.1) K e Sofc a. h je vyhlazovací parametr nebo také šířka vyhlazovacího okna. Jádrový odhad hustoty závisí na třech parametrech: jádře, které hraje roli váhové funkce, vyhlazovacím parametru, který řídí hladkost odhadu, a na řádu jádra, který odpovídá předpokládanému počtu derivací neznámé hustoty. Popíšeme konstrukci jádrového odhadu. V každém bodě Xi sestrojíme jádro Kh a odhad v bodě x je průměr n hodnot jader v tomto bodě - viz obrázek 3.2(a). Nyní uvedeme ještě vztah pro jádrový odhad v-té derivace hustoty. Budeme předpokládat, žeO < v < k — 2 a k, v jsou stejné parity. Pak 1 n Y (3.2) Konstrukce jádrového odhadu derivace je stejná jako konstrukce odhadu hustoty - obr. 3.2(b) - pouze místo jádra K e Sq/c používáme jádro třídy Svk- 3 Statistické vlastnosti jádrových odhadů hustoty Stejně jako u jádrových odhadů regresní funkce lze kvalitu jádrového odhadu hustoty popsat lokálně pomocí střední kvadratické chyby. 29 0.8 (a) Hustota -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 (b) Derivace hustoty Obrázek 3.2: Konstrukce jádrového odhadu hustoty a její derivace Věta 3.1. MSE f(x,h) = E(f(x,h) - f(x)f = -((^ * f)(x) (Kh * /)2(x)) + ((Kh * /)(*) - /(x)V Důkaz. Spočítejme střední hodnotu odhadu f(x,h) 1 " /" ä) = E- Kh(x - Xi) = ~X)= Kh(x - y)f(y) dy = * f)(x). i=0 Vychýlení (bias) pak bude mít tvar bias/(x, h) — Ef(x, h) — f(x) — (Kh * f)(x) — f(x). Dále upravíme vztah pro rozptyl 1 " 1 " 1 var f(x, h) — var — N Kh(x — Xi) — — var N Kh(x — Xi) — — var Kh(x — X) i=l i=l = -EKlix -X)- - ÍEKh(x - X))2 n n 1 ľ Kl{x-y)f{y)dy--({Kh*f){x))2 n .1 n 1 ((K2*f)(x)-(Kh*f)2(x)). □ Důsledek. mise = I (/"(Kč * /)(a.) da. _ ľ(Kh * /)2(a.) da.^ + /*^ t /)(a;) _ }{x)y dx Podobně jako u odhadu regresní funkce můžeme použít globální pohled na kvalitu odhadu, a to pomocí střední integrální kvadratické chyby (mise) a jejího asymptotického tvaru (amise). Věta 3.2. Nectí funkce f má spojité derivace až do řádu k0 (tj. / e Ck°) pro 0 < k < k0, K e Sok a j(/(fe) (x)) dx < oo, dále předpokládejme h —>• 0 a nh —>• oo pro n —>• oo. PflA: pZarí MISE/(■,/!) = MISE(/i) = + ^h2k(32(K)V(f^) + o(h2k + (nh)-1), kdeV(f^) = J(/(fe)(x))2da;. 30 Důkaz. Nejprve vypočteme střední hodnotu Ef{x,h) = (Kh*f)(x) = j Kh(x-y)f(y)dy = J l-K(^-) f {y) dy K(z)f(x — h z) dz dále použijeme Taylorův rozvoj: f(x — hz) — f(x) — f'(x)hz + • • • + ^~fcV f^k\x)hkzk + o(h K(z) [f(x) - f'(x)hz + ■■■+ (-^f(-k\x)hkzk + o(hk)] dz = m + {~l)kf^){x)hkh{K) + o(hk). Tedy vychýlení odhadu je tvaru Ef(x, h) - f(x) = {-1)kfk^{x)hk(3k(K) + o(hk). Odtud plyne, že Ef(x, h) = f(x) + o(l). Nyní dokážeme vztah pro rozptyl. Víme, že a dále počítáme :f(x,h) = ^{{Kl * f)(x) - (Kh * f)2(x)) 1 '^(V9/(y)^-^/(ic)+o(1))2 =/(z)+o(l) \ í K\z){f{x) + o{l))dz--{f{x)+o{l)Y nh / n Tedy ^ / if2(z)dz + o((n/l)-1). nh /(*) /"^^ , „,,„m-^ , í(-l)kf(k)(x)uk, MSE/(x, /i) = / A"2(*)dz + o((n/0_1) + i-; v ' hk/3k(K) + o(hk) nh \ k\ a pak MISE/(-,/i) = / MSE/(x, h) dx ™ + -^ti*Pl{K)V• 0, nh —>• oo pro n —>• oo, paA: / je konzistentním odhadem /, í/. E1/ fl var / —>• 0. 31 Stejně jako u odhadu regresní funkce má význam asymptotická integrální střední kvadratická chyba AMISE/(■, h): MISE/(-, /i) = AMISE/(■, h) + o(h2k + (nh)-1), kde AMISE je tvaru AMISE/(-, h) = AMISE(/0 = + J^h2kfá(K)V(fM). (3.3) V dalších částech textu budeme využívat označení jednotlivých částí chyby AMISE, která je součtem asymptotického tvaru integrálu rozptylu AIV (asymptotic integrated variance) a asymptotického tvaru integrálu druhé mocniny vychýlení AISB (asymptotic integrated squared bias): AlV(h) = W- AISB(/i) = j^h2k(32(K)V(f^), tedy AMISE (h) = AIV(h) + AISB(/i). Jr^2/(2fc+l)aJ?2fc+1 _ _ Užitím vztahů T(K) = (l/(^)fe/3fc(^))2/(2fe+1) a = pro £ľ G S0k lze AMISE zapsat ve tvaru AMISEW=T(/0(^ + ^p). 0.4) Odtud je zřejmé, že vyhlazovací parametr, pro nějž AMISE nabývá minimální hodnoty, je dán vztahem > 2fc+l _ ^0fc + 1(^Q2 °Pt'°'fe 2knV(f^)' tj. Äopt.o.fc = 0(n-1/(2fe+D). Vypočtěme hodnotu AMISE při dosazení optimálního parametru /iOpt,0,fe: AMISE(= '^VW1' (2fc(fc^/(2fc+1) ■ (3-5) tj. AMISE(/loptj0,fc) = 0(n-2fe/(2fe+D). I v tomto případě, podobně jako u odhadhu regresní funkce, platí vztah mezi AIV(h) a AISB(/i): AIV(/iopt,o,fc) = 2fcAISB(/iopt,0,fc)- (3.6) Nyní uvedeme zajímavou vlastnost vyhlazovacího parametru. Poznámka 3.1. Nechť Ä" e So2- Pak optimální hodnota vyhlazovacího parametru je h5 - 02 opt'0'2 nv(/") ■ Počítejme derivace AMISE (3.3) í^ = Ä+ewä(Jf)v(n. Řešením rovnice d3 AMISE / dh3 — Oje V(K) č,5 02 nP2{K)V{f") nV(f") opt'0'2' tj. hopt02 také realizuje minimum d2AMISE / dh2. 32 Obecně lze ukázat, že d2AMISE(/(-,/i)) dh2 , _ 2k-2 N — 0{n 2k+1 / ^—hopt,0,k a to znamená, že pro jádra vyšších řádů je minimum AMISE plošší a tedy volba h blízká optimální hodnotě /i0pt,o,fe nevede k velkému růstu AMISE (obr. 3.3). Obrázek 3.3: AMISE pro jádra vyšších řádů s vyznačenými minimálními hodnotami pro jádra řádu 2,4, 6 Vztah pro optimální hodnotu vyhlazovacího parametru poskytuje informaci, že asymptoticky je /i ~ n-1/(2fe+1). Ale vztah má pouze teoretický charakter, protože optimální parametr závisí na neznámé hustotě /. Poznámka 3.2. Z předchozích úvah je zřejmé, že množina přípustných hodnot vyhlazovacích parametrů je dána vztahem H„ = [afen-Wi);&fen-i/(2fe+i)]; kde cifc, bk jsou konstanty, 0 < ak < bk < oo. 3.1 Odhad derivace hustoty Pojednáme nyní stručně o statistických vlastnostech jádrových odhadů derivace hustoty. Připomeňme, že jádrový odhad derivace hustoty je dán vztahem (3.2), tj. 1 n Y i=l Předpokládejme nyní, že platí 0 < v < k — 2, linin^oo h — 0, linin^oo nh + — oo, / e C 0 (k < ko) a V(f^) — f (f(k\x))2 dx < oo. Pak lze ukázat, že asymptotická střední kvadratická chyba AMISE /»(■, h) je tvaru AMISE/»(,/*) = + ^/^-^(tfM^/W). Důkaz je založen na použití vhodného Taylorova rozvoje hustoty /, podobně jako u důkazu tvaru AMISE u odhadu hustoty. Optimální hodnota vyhlazovacího parametru je dána vztahem fe+1 5lk+\2v + l)(klf V(KM) °P^k 2n{k-v)V{f^) ' vk 33 Tento vzorec umožňuje výpočet optimálního vyhlazovacího parametru pro pomocí /iOpt,0,fe a hopt,i,k- Předpokládejme nejdříve, že v a k jsou sudá čísla. Pak hopt,,,k _ /(2^+l)AA1/(2fe+1) 6: uk 'lopt,0,k \ k — V ) Sok Pro ľ a k lichá platí hopt,,,k _ ({2v+\){k-\)\1/{2k+1) Ô; hopt,i,k V -v) J ôik Speciálně pro v — 2, k — 4 dostáváme velmi užitečný vztah -,1/9^24 i lopt,2,4 — J-l přičemž KoptM*) = ^(*2 - !)(7a;2 - 3); ^04 = 2,0165, 105 ^(2)(x) = KMA{x) = —(1 - x2)(5x2 - 1), Ô24 = 1,3925. (3.7) (3.8) hopt,2,4 — 101/'972^/iopt,0,4, (3-9) 004 4 Volba jádra Volba jádra není z asymptotického hlediska podstatná, jak je zřejmé z faktu (3.5). Je vhodné zvolit optimální jádro, které minimalizuje fukcionál T(K), neboť tato jádra jsou spojitá na M a hladkost jádra také „zdědí" odhadovaná hustota. 5 Volba vyhlazovacího parametru Volba vyhlazovacího parametru pro jádrový odhad hustoty je, stejně jako u regrese, zásadním problémem. I když tomuto problému byla a je věnována značná pozornost, doposud neexistuje univerzální přístup k řešení tohoto problému. Nejjednodušší metoda je „okometrická". Je účelné „nakreslit" několik křivek s různými vyhlazovacími parametry dříve než užijeme nějakou automatickou proceduru. Je třeba zdůraznit, že z hlediska analýzy všechny volby vyhlazovacího parametru vedou k užitečnému odhadu hustoty. Velká šířka okna charakterizuje globální strukturu hustoty a naopak malá šířka odhaluje lokální strukturu, která může nebo nemusí být přítomná v přesné hustotě. Tuto myšlenku ilustruje obrázek 3.4, na němž jsou zobrazeny odhady pro simulovaná data (n — 100) s hodnotami vyhlazovacího parametru z intervalu [0,05,1]. Jednotlivé odhady hustoty příslušející těmto hodnotám jsou znázorněny tenkými čarami. Silná křivka znázorňuje odhad s optimální hodnotou h — 0,4217. Třída těchto odhadů ukazuje široký rozsah vyhlazení od pod-hlazení až k přehlazení. V dalších úvahách bude užitečná následující definice: Definice 5.1. Nechť h je odhad /iOpt,0,fe- Řekneme, že h konverguje k /iOpt,0,fe s relativní rychlostí ■nTa, jestliže hopt,0,k 34 Obrázek 3.4: Volba vyhlazovacího parametru 5.1 Metoda referenční hustoty Nejčastěji se pro odhad neznámé veličiny V(f^) (viz rovnice (3.3)) používá parametrické třídy hustot. Jednou z možností je použít standardní normální hustotu / s rozptylem a2, tj. předpokládáme, že 1 m — V tomto případě je odhad optimálního vyhlazovacího parametru tvaru pro K e Sofc %ef = i,,/ ) S0k_1 je standardní normální kvantilová funkce a číslo X^n/^, respektive X[„/4], je horní, respektive dolní výběrový kvartil. Je vhodné volit mm{dsD,&iQii}-Poznámka 5.1. Pokud za jádro K zvolíme Gaussovo jádro (k — 2) pak dostaneme jednoduchý vztah / 4 \ V5 ZlREF = ( ^~ CT- (3-13) \3n< Příklad 5.1. Použijme odhad vyhlazovacího parametru pro data z příkladu 2.1 metodou referenční hustoty. Pro Epanečnikovo jádro, které je řádu k — 2, se vztah (3.10) zjednoduší na tvar /o f=\ 1/5 "ref — 3~ i 02 Dále odhadneme směrodatnou odchylku: cťsd — 1,0325, viqr — 1,3344, tedy a — 1,0325. Po dosazení počtu prvků n — 100 a parametru S02 — 1,7188 získáme hodnotu vyhlazovacího parametru pro odhad hustoty /iref — 0,9639. Na obrázku 3.5 je vykreslen odhad hustoty s tímto parametrem. 35 Obrázek 3.5: Odhad hustoty s /iref — 0,9639, odhad (—), původní funkce (—) 5.2 Metoda maximálního vyhlazení Princip maximálního vyhlazení (maximal smoothing) - ms (nebo přehlazení) znamená, že vybereme největší stupeň přehlazení kompatibilní s odhadovanou hustotou. Získáme tak horní hranici pro odhad optimální šířky vyhlazovacího okna. Tato hodnota pak může sloužit jako počáteční aproximace pro některé z dalších metod. Princip spočívá v tom, že hledáme hustotu, pro kterou V(f^) nabývá minimální hodnoty, a tedy vztah pro /iOpt,0,fe nabývá maximální hodnoty. Věta 5.1 (Terrell 1990). Mezí všemi hustotami f s nosičem [—1,1] má hustota rozdělení Beta(k+2, k + 2) ř (2fc+3)! Cl x2)k+1 \x\<í gk(x) = \ (fe+l)!22fe+3^ 1 1 1 - ' 10 jinak, nejmenší hodnotu integrálu (/(fe) (x))2 dx. Lze ukázat, že platí L al = I-i x29k{x) dx = 2^5. 2. Pro r > 0, f (rg(k\rx))2 dx — r2k+1 J (g^ (x))2 dx pro každou hustotu, pro kterou integrál existuje. 3. Jestliže hustota g má rozptyl a2, pak hustota ^-g(^-x} má rozptyl a2. (Podrobněji např. [5].) Jestliže / je neznámá hustota s rozptylem a2 a gi~ je hustota rozdělení Beta(/j + 2, k + 2), pro kterou je V(g^) minimální, pak lopt,0,k Hodnotu cr lze odhadnout pomocí dříve uvedených vztahů a Ofe = 2fc1|_5. Hodnotu V(g^) lze vypočítat pomocí speciálních ortogonálních polynomů [5]: v(Q^) ŕ (Q^(x))2dx _(2fc+3)!(2fc+2)!_ V(9k f-J^x dx~ 2^(2k+l)(2k + 5)(k + l)p- (3.14) Použijeme-li poslední vyjádření (3.14), dostaneme horní hranici pro vyhlazovací parametry hopt,o,k < hMS = dn-^^+^bk, 36 pncemz ^/——/22k+2V(K)(2k + l)(2fc + 5)(fc + l)2(fc!) ^(Ä-)(2A: + 3)!(2fe + 2)! Tabulka 3.1: Hodnoty bk pro optimální jádro Koptto,k <= Sofc 10 2,5324 3,3175 3,9003 4,3949 4,8349 Příklad 5.2. Určeme hodnotu Iims Pro odhad hustoty s jádrem řádu k — 2, tj. K e Sofc-Podle věty 5.1 je hustota g% tvaru 35 92{x) = —(1 -x2)3, a; e [-1,1], a dále z vlastností funkce #2 a ze vztahu (3.14) plyne 1 f „ 2 x g2(x) dx — — a / dx = 35. Pak pro Kopt^j, hus = 1/E 1/5fV(K) 243 xV5 /322(^) 35 Příklad 5.3. Pro data z příkladu 2.1 bude vyhlazovací parametr určený metodou maximálního vyhlazení s Epanečnikovým jádrem (k — 2, V(K) — 3/5, ^(K) = 1/5) roven ^ _1/5/ 3/5 243 \ VŠ ,0409, ,1/25 35 protože a — 1,0325 a n — 100. Výsledný odhad je vidět na obr. 3.6. 0.6i-1-1-1-1-1— Obrázek 3.6: Odhad hustoty s hyis — 1,0409, odhad (—), původní funkce (—) Poznámka 5.2. Hodnota hyis může sloužit jako horní hranice pro množinu vyhlazovacích parametrů volených podle jiné metody, např. metody křížového ověřování. Tedy Hn — [hg, /ims]/ kde hi je nejmenší vzdálenost mezi po sobě jdoucími body Xi, i — 1,..., n. 37 5.3 Metoda křížového ověřování Metoda křížového ověřování patří mezi nejužívanější metody pro odhad vyhlazovacího parametru. Myšlenka této metody je založena na minimalizaci MISE, jak je zřejmé z následující úvahy: MISE/(-,/0=£ I (f(x,h)-f(x))2dx = E / f2(x,h)dx -2E / f(x,h)f(x)dx+ / f(x)dx MISE f (x, h) -j f(x)dx = E\^j f(x,h)dx-2 j f(x,h)f(x)dx Definujme funkci krížového overovaní ľ 2 " CV(/i) = J f (x, h)dx--J2 f-ÁXi, h), (3.15) i=l kde f-i (Xi, h) je odhad v bodě Xi bez použití tohoto bodu. Věta 5.2. Platí ECV(h) = WSEf(;h)- J f(x)dx t). CV(h) je nevychýleným odhadem El f2(x,h)dx-2 / f(x,h)f(x)dx Důkaz. Střední hodnota prvního členu rovnice (3.15) je zřejmá, potřebujeme spočítat střední hodnotu druhého členu tohoto vyjádření. 1 n 1 n 1 n E- f-AXi, h) = E-J2-7 E Kh(xi - xj) n '-^ n '-^ n — 1 Odtud Í — 1 Í — 1 j — 1 . n n E-7T E E Kh(xi - xj) = EK^X, - X2) n(n — 1) '-^'-^ v ' i=i j=i Kh(x -y) f (y)J (x) dxdy = E J f (x, h) f (x) dx. E f (x,h) ECV(h)=E J f2(x,h)dx-2E J f(x,h)f(x)dx. □ Odhad /ioptío,fc je dán vztahem/icv — argmin^g^ C V (h). Odtud plyne, že C V (h) + j f2 (x) dx je pro každé h nevychýleným odhadem MISE(/i). Protože j f2 (x) dx nezávisí na h, minimalizace EGV(h) odpovídá minimalizaci MISE. Jestliže předpokládáme, že minCV(/i) ~ min E GV(h), dostaneme dobrou aproximaci optimální hodnoty h. Poznámka 5.3. Předpokládejme, že k — 2. Pak vychýlení odhadu může být velké, jestliže /" nabývá velkých hodnot, tj. křivost hustoty je velká. Při vyhlazování se tato objevuje ve „vrcholech", kde je vychýlení záporné, nebo v „údolích", kde je vychýlení kladné. Odhad má tendenci „vyhladit" tyto jevy, jak je patrné z obrázku 3.8. Příklad 5.4. Jádrový odhad hustoty dat z příkladu 2.1 je zobrazen na obr. 3.9. Pro rekonstrukci bylo použito Epanečnikovo jádro a vyhlazovací parametr určený metodou křížového ověřování hcv = 0,5628. 38 o h 1 Obrázek 3.7: Porovnání minima MISE (červenou) a minima funkce křížového ověřování CV (modrou) pro simulovaná data Obrázek 3.8: Zahlazení vrcholů a údolí při odhadu hustoty směsi normálních rozdělení -2-1 0 1 2 3 4 Obrázek 3.9: Odhad hustoty s hcv — 0,5628, odhad (—), původní funkce (—) 39 5.4 Iterační metoda Připomeňme nyní vztah (3.6): Přepíšeme tuto rovnici AIY(/iopt,o,fc) = 2A; AISB(/iopt,o,fc). nh \k\y (3.16) Problém nalézt hoptaM, pro které AMISE /(-, h) nabývá minimální hodnoty, je tedy ekvivalentní řešení této rovnice. Zde se ovšem vyskytuje stejný problém - neznáme hodnotu V(f^), a proto budeme počítat s odhady rozptylu a vychýlení. Tyto odhady uvažujeme ve tvaru var/(x, h) = ^ J K2 (y) f (x - hy) dy bias/(x, h) — bias/(x, h) — (K * /)(.r, h) — f (x, h) = f(x-hy,h)K(y)dy-f(x,h). Odtud plyne AISB/(-,/i)= / ( I f(x-hy,h)K(y)dy-f(x,h)j dx výraz lze upravit pomocí konvolucí a dostaneme = — V (K * K * K * K - 2K * K * K + K * K) 1 1 ' J dx i J = l 1 " 4XA Xt — X; Funkce A(z) — (K * K * K * K — 2K * K * K + K * K)(z) má tyto vlastnosti zJ'A(2)dr = 0, j = 0,l,...,2A:-l, 2^ fcA(z)dz= (^)/^, 1 . /r AISB/(-, /i) je vychýleným odhadem AISB, a proto budeme uvažovat _ i n AISB/(-, M = — £ Aft(Xi - X,). 40 Místo rovnice (3.16) řešíme rovnici nh Tuto rovnice lze také přepsat ve tvaru: h= nV}K)--. (3.17) Uvedenou rovnici lze řešit Newtonovou metodou, neboť derivaci funkce lze snadno spočítat užitím konvolucí. Řešení této rovnice označíme hn ■ Řešení rovnice (3.17) lze považovat za vhodnou aproximaci hopt0j„. Tato skutečnost je dokázána v následující větě [3]. Věta 5.3. Nectí Pak platí EP{h) = P(h) + o(h2k+1), vavP(h) = ^-V{K)V{f) + oin^h-1). nzh Poznámka 5.4. Rychlost konvergence pro metodu křížového ověřování: hcv - hopt,o,2 ^ o(n~1/10) hOpt,0,2 Rychlost konvergence pro iterační metodu hiT - hopt,o,2 ^ o(n~1/10) hOpt,0,2 Řády rychlosti konvergence pro CV metodu a iterační metodu jsou stejné, ale výhodou iterační metody je menší výpočetní náročnost. Příklad 5.5. Jádrový odhad hustoty dat z příkladu 2.1 s Epanečnikovým jádrem a vyhlazovacím parametrem určeným iterační metodou je uveden na obr. 3.10. Shrneme-li na závěr doposud vypočítané hodnoty vyhlazovacích parametrů pro simulovaná data z příkladu 2.1, můžeme vizuálně porovnat jednotlivé odhady - viz obrázek 3.11. Vt,o,2 = 0,5122 hREF = 0,9639 hMS = 1,0409 hcv = 0,5628 hIT = 0,5314 41 -2 -1 4 Obrázek 3.10: Odhad hustoty s h — 0,5314, odhad (—), původní funkce (—) 6 Automatická procedura Obdobně jako v případě regresní funkce můžeme nalézt podobnou formuli pro AMISE/(-, h), ve které budou jednotlivé parametry K, h, k separovány, což nám umožní navrhnout proceduru pro simultánní volbu těchto parametrů. Vyjdeme ze vztahu AMISE(W).W(^ + ^p) Ze vztahu pro /iopt,o,fe vypočteme V(f^) J0k_ 2knhopt0k a tuto hodnotu dosadíme do předchozího vztahu: (2k+l)ô0k AMISE(/iopt,o,fc) = T{K) 2nkh, opt ,0,/c Tento vztah je základem automatické procedury, označme jej L(k) ve shodě s označením u regresní funkce. Podobně množinu vhodných řádů k označme kn /(fco) = {2j,j = 0,...,[-^]} Krok 1 Pro k e I(ko) najděte optimální jádro Koptto,k <= Sok, které je dáno tabulkou 1.2, k němu příslušný kanonický faktor Sok- Krok 2 Pro k e I(ko) a Koptto,k <= Sofc najděte optimální vyhlazovací parametr hoptto,k- Krok 3 Pro k e I(ko) vypočtěte hodnotu výběrového kriteria L(k) s využitím hodnot získaných v krocích 1 a 2. Krok 4 Vypočtěte optimální hodnotu řádu k, které minimalizuje funkcionál L(k). Krok 5 Použijte parametry z předchozích kroků ke konstrukci optimálního jádrového odhadu hustoty, tj. opt,0,k i—l opt,0,k 42 (a) Odhad s ftopt,o,2 0.6i-1-1-1-1-1-1 0.6 -2 -1 0 1 2 3 4 -2 -1 0 1 2 3 4 (b) Odhad s Aref (c) Odhad s /ims (d) Odhad s hcv (e) Odhad s hn Obrázek 3.11: Srovnání odhadů pro data z příkladu 2.1 Příklad 6.1. Aplikace procedury na data z příkladu 2.1. Maximální řád jádra zvolme ko — 8, tedy množina možných řádů jader je 1(8) — {0, 2, 4, 6, 8}. Pro tyto řády spočítejme hodnoty z kroků 1-3. V toolboxu Matlabu, který je doprovodným materiálem těchto skript, je jako implicitní metoda pro odhad vyhlazovacího parametru automatickou procedurou použita iterační metoda (podrobněji např. [3]). Proto při výpočtu optimálních parametrů použijeme mezivýpočty z tohoto toolboxu. 43 k K opt,o,k Sok h L(K) 2 -|(x2-l) 1,7188 0,5314 0,0141 4 i§(x2 - l)(7x2 - 3) 2,0165 1,0734 0,0131 6 -i2|(x2-l)(33a;4-30a;2 + 5) 2,0834 1,6460 0,0125 8 J^(x2-l)(715x6-100Lr4 + 385:r2 - 35) 2,1021 2,1367 0,0126 Z tabulky vidíme, že optimální řád jádra je k — 6. Výsledný odhad je uveden na obrázku 3.12. 0.6 0.1 0.5 0.4 0.3 0.2 4 Obrázek 3.12: Simulovaná data (x) s jádrovým odhadem hustoty při použití procedury (—) a původní funkcí (—) Při bližším pohledu na odhadnutou hustotu je patrný vliv použití optimálního jádra vyššího řádu. Jádra vyšších řádů mohou nabývat záporných hodnot a tím ovlivnit výslednou odhadnutou funkci - viz obrázek 3.13. V takovém případě je vhodné použít jinou metodu pro nalezení vyhlazovacího parametru, případně použít jiné jádro. Lze doporučit jádra třídy S02, např. kvar-tické jádro: K{x) — j|(l — x2)2I[_ltl](x), nebo jádro triweight: K{x) — ||(1 — x2)3I[_ltl](x). Datový soubor obsahuje morfologická měření padesáti exemplářů od obojího pohlaví a obou barevných forem (oranžová a modrá) krabů rodu Leptograpsus.1 Pro odhad hustoty nám postačuje jeden druh měření, vybrali jsme délku podél středové osy krunýře, která byla měřena v milimetrech. Data jsou uvedena v tabulce 6.7. Užitím výše uvedených metod pro odhad vyhlazovacího parametru jsme (při použití Epaneč-nikova jádra) dostali následující hodnoty: U automatické procedury je v toolboxu implicitně nastavena iterační metoda pro odhad vyhlazovacího parametru, proto jsme tuto volbu ponechali i zde, af má čtenář možnost porovnání při vlastních výpočtech. Při použití procedury vyjde vyhlazovací parametr roven hpľoc — 31,7329 s optimálním jádrem Kopt,o,8- Výsledné odhady hustoty na jsou zachyceny na obrázku 3.14. Celý datový soubor je dostupný v programu R. 7 Aplikace na reálná data /iREF = 5,7856, hMS = 6,2480, hCv = 8,1317, hIT = 6,8263. 44 (a) Jádro A"optio,6 (b) Odhad hustoty Obrázek 3.13: Jádro třídy S06 a k němu příslušný odhad hustoty při použití procedury (—) a původní funkcí (—) 45 46 Shrnutí Odhad hustoty pravděpodobnosti / v bodě x je tvaru ' x — Xi nh i=l Asymptotická střední kvadratická chyba jádrového odhadu hustoty pravděpodobnosti s jádrem řádu k je součtem asymptotického tvaru rozptylu (AIV) a druhé mocniny vychýlení (AISB) AMISE(/0 = + J^h2kfá(K)V(fW). aiv " aisb ' Optimální vyhlazovací parametr vzhledem k AMISE pro odhad hustoty pravděpodobnosti s jádrem řádu k je tvaru . 2fe+i °ofc y^-i °Pt'°'fe 2knV(f(k))1 tj. Ziopt.o.fc = 0(n-1/(2fe+D), AMISE(/lopt,o,fc) = 0(n-2fe/(2fc+D). Metody pro odhad optimální hodnoty vyhlazovacího parametru h • metoda referenční hustoty f22kkrj^\^\ __i_ %EF = {-l2kjU-) Ô0k(7n 2k+1' • metoda maximálního vyhlazení hMS = *n-W°+»bk, • metoda křížového ověřování r 2 ™ - hCv = arg min CV(/i) = / f2(x, h) dx--V" f-l(Xl, h) i=i iterační metoda hrr = pevný bod funkce: h = 2k^J^Xi-Xj) Automatická procedura pro simultánní volbu optimálního jádra, vyhlazovacího parametru a řádu jádra je dostupná v toolboxu Matlabu. Dodatky a cvičení 1. Dokažte vztah (3.4) pro tvar chyby AMISE. 2. Dokažte vztah (3.13). 3. Spočítejte (3.10). 47 4. Spočítejte /ims pro • Epanečnikovo jádro K(x) — |(1 — x2), je-li = (I/gp- = 15, . kvartické jádro Eľ(x) = ±§(1 - x2)2, je-li = ^ = 35. 5. Aplikujte metody pro odhad vyhlazovacího parametru a automatickou proceduru na simulovaná i reálná data. 48