Úvod Vyhlazování je statistická technika pro rekonstrukci reálné funkce na základě pozorovaných nebo naměřených dat. Cílem vyhlazování je nalezení takového odhadu neznámé funkce, aby byly odfiltrovány náhodné výkyvy a bylo možné lépe poznat strukturu dat. K tomuto úkolu lze přistoupit dvěma způsoby - parametricky a neparametrický: • Parametrické odhady jsou založeny na předpokladu, že neznámá funkce patří do třídy funkcí závislých na parametrech, a cílem je odhadnout tyto parametry. • Neparametrické odhady nepředepisují datům „Prokrustovo lože" parametrizace, ale nechávají „hovořit samotná data". V tomto učebním textu se zaměříme na neparametrické odhady, a to zejména na jádrové odhady, které patří mezi efektivní neparametrické odhady. Budeme se zabývat jádrovými odhady regresní funkce, hustoty, distribuční funkce a také odhadem dvourozměrné hustoty. Všechny jádrové odhady závisí na jádře, které má roli váhové funkce, a na vyhlazovacím parametru, který řídí hladkost odhadu. Budeme zabývat následujícími otázkami: • Jaké jsou statistické vlastnosti jádrových odhadů. • Jaký vliv má tvar jádra na odhad. • Jaký vliv má šířka vyhlazovacího okna na odhad. • Jak lze tuto šířku stanovit v praxi. Volba vhodného vyhlazovacího parametru je zásadním problémem ve všech typech jádrových odhadů a tomuto problému budeme věnovat značnou pozornost. Všechny uvedené metody jsou implementovány v Matlabu, příslušný toolbox je dostupný na adrese: https://www.math.muni.cz/veda-a-vyzkum/vyvijeny-software/ 2 7 4-matlab-toolbox.html Na konci textu (kapitola 7) jsou uvedeny soubory dat pro samostatnou práci studentů. Tyto soubory již byly zpracovány v příslušných kapitolách a studenti si tak mohou ověřit správnost svých výsledků. Definice základních statistických pojmů a jejich vlastností lze najít např. v elektronických skriptech Pravděpodobnost a statistika I autorů M. Forbelské a J. Koláčka (jsou dostupná na Elportálu Informačního systému). Na tomto místě bych ráda poděkovala Mgr. Kamile Hasilové, Ph.D., za pomoc při sazbě tohoto textu a za příspěvek ke kapitole 5 a kapitolám o reálných datech. 1 Obsah 1 Jádrové funkce a jejich vlastnosti 6 1 Základní pojmy a definice................................. 6 1.1 Jádra s minimálním rozptylem.......................... 7 1.2 Optimální jádra................................... 8 2 Jádrové odhady regresní funkce 11 1 Motivace.................. ......................... 11 2 Základní typy neparametrických odhadů........................ 12 3 Statistické vlastnosti jádrových odhadů......................... 16 4 Volba jádra.......................................... 21 5 Volba vyhlazovacího parametru.............................. 22 5.1 Metoda křížového ověřování........................... 22 6 Automatická procedura .................................. 24 7 Aplikace na reálná data................................... 25 3 Jádrové odhady hustoty 30 1 Motivace........................................... 30 2 Základní typy neparametrických odhadů........................ 30 3 Statistické vlastnosti jádrových odhadů hustoty .................... 32 3.1 Odhad derivace hustoty.............................. 36 4 Volba jádra.......................................... 37 5 Volba vyhlazovacího parametru.............................. 38 5.1 Metoda referenční hustoty............................. 38 5.2 Metoda maximálního vyhlazení ......................... 39 5.3 Metoda křížového ověřování........................... 41 5.4 Iterační metoda................................... 43 6 Automatická procedura .................................. 45 7 Aplikace na reálná data................................... 49 4 Jádrové odhady distribuční funkce 53 1 Motivace.................. ......................... 53 2 Základní typy neparametrických odhadů........................ 54 3 Statistické vlastnosti odhadu ............................... 55 4 Volba jádra.......................................... 59 5 Volba vyhlazovacího parametru.............................. 59 5.1 Metody křížového ověřování........................... 59 5.2 Princip maximálního vyhlazení.......................... 59 5.3 Plug-in metoda................................... 60 6 Aplikace na reálná data................................... 61 2 5 Jádrové odhady dvourozměrných hustot 66 1 Motivace........................................... 66 2 Základní typy odhadů................................... 67 3 Statistické vlastnosti jádrových odhadů hustoty .................... 68 4 Volba jádra.......................................... 71 5 Volba vyhlazovacího parametru.............................. 71 5.1 Metoda referenční hustoty............................. 71 5.2 Metoda křížového ověřování........................... 72 6 Aplikace na reálná data................................... 73 6 Návody ke cvičením a odpovědi na otázky 78 7 Datové soubory 85 8 Přílohy 99 3 Seznam použitého značení h H to(-) /(•) F(-) h a m f F I K(-) V (K) V(g) h{K) f*9 W(-) Ck[0,1] NW PC LL GM vyhlazovací parametr matice vyhlazovacích parametrů regresní funkce hustota pravděpodobnosti distribuční funkce odhad vyhlazovacího parametru h odhad směrodatné odchylky a odhad regresní funkce m odhad hustoty / odhad distribuční funkce F značí integrál , pokud není uvedeno jinak jádrová funkce (jádro) V(K) = / K2{x) dx V(g)=fg2(x)dx fik(K) = / xkK(x) dx konvoluce funkcí / a g, (/ * g)(x) — J f (ť) g (x — t) dt integrál z jádra, W(x) — f* K {ť) dt prostor funkcí, které mají spojité derivace až do řádu k včetně na intervalu [0,1] Nadarayovy-Watsonovy odhady Priestleyovy-Chaovy odhady lokálně lineární odhady Gasserovy-Múllerovy odhady 4 MSE střední kvadratická chyba (mean square error) MISE střední integrální kvadratická chyba (mean integrated square error) AMISE asymptotická střední integrální kvadratická chyba (asymptotic mean integrated square error) AIV asymptotický tvar integrálu rozptylu (asymptotic integrated variance) AISB asymptotický tvar integrálu druhé mocniny vychýlení (asymptotic integrated square bias) AMSE střední průměrná kvadratická chyba (average mean square error) CV metoda křížového ověřování REF metoda referenční hustoty MS metoda maximálního vyhlazení PI plug-in metoda IT iterační metoda Symbolika O, o Symboly O a o se často používají pro vyjádření chyb matematických výrazů. (Nazývají se také řádové vztahy.) Nechť / je funkce definovaná v okolí bodu a. Symbol znaci, ze Podobně symbol znaci, ze f(x) — 0(g(x)) pro x —>• a lim sup ———- < oo. x^a g[x) f(x) — o(g(x)) pro x —>• a lim M = 0. >a g(x) Dodatek „pro x —>• a" se často vynechává, pokud je jasné, o které a se jedná. Je to zejména v případech a — 0 či a — oo. Často používaným výraz je 0(hk), resp. o(hk), kde g(h) — hk, přičemž zpravidla h —>• 0. Pro počítání s výrazy obsahujícími symboly O a o platí následující pravidla: 0(g(x)) + 0(g(x o(g(x)) + o(g(x 0(g(x)) + o(g(x 0{g{x)) ■ 0{h{x o(g(x)) ■ o(h(x 0(g(x)) ■ o(h(x o(g(x 0(g(x)) o(g(x)) 0(g(x)) 0(g(x) ■ h{x)) o(g(x) ■ h(x)) o(g(x) ■ h(x)) 0{g{x)) Pozor! Tyto rovnice nejsou symetrické, platí jen zleva doprava. Např. poslední rovnice značí, že je-li funkce f(x) — o(g(x)), pak je f(x) — 0(g(x)). Opačně to ovšem neplatí. 5 Kapitola 1 Jádrové funkce a jejich vlastnosti 1 Základní pojmy a definice V úvodu bylo uvedeno, že všechny jádrové odhady závisí na jádrové funkci (jádře), a proto se v této kapitole budeme zabývat jádrovými funkcemi. Nyní uvedeme definici jádra a jeho vlastnosti. Definice 1.1. Nechf v, k jsou nezáporná celá čísla, 0 < v < k, nechf K je reálná funkce s těmito vlastnostmi 1. K splňuje Lipschitzovu podmínku na intervalu [—1,1], tj. \K(x) — K(y)\ < L\x — y\ pro Vx,ye [-1,1],L>0, 2. nosič(.?f) = [-1,1], tj. K = 0 vně intervalu [-1,1], 3. K splňuje momentové podmínky i [(-1) 1/1 3 = v a xkK{x) dx ^ 0, tuto hodnotu označíme f3k{K). Taková funkce K se nazývá jádro řádu k a třída všech těchto funkcí se označuje Svk-Definice 1.2. Označme V(K) — j\ K2{x) dx, K e Suk pro jádro e ívfc a definujme veličinu ô,. k Tuto veličinu budeme nazývat kanonický faktor. Příklad 1.1. V tabulce 1.1 jsou uvedeny příklady několika jader společně s jejich grafy. Funkce I[-i,í\ (x) je indikátorová funkce intervalu [—1,1], tj. 1 pro x e [— 1,1], 0 jinak. Poznámka 1.1. V praxi se také používá Gaussovo jádro K(x) = ^=e"a;2/2, xeR. Jde sice o jádro řádu 2, avšak nemá vlastnost 2 z definice 1.1. 6 Tabulka 1.1: Jádra S 02 obdélníkové jádro K(x) = l(l-x2)I{_1A](x) Epanečnikovo jádro K(x) = &il-X2)2I[_ll]ix) kvartické jádro Nyní uvedeme dva typy jader - jádra s minimálním rozptylem a optimální jádra. 1.1 Jádra s minimálním rozptylem Předpokládejme, že K e S^fc, 0 < ^ < £; — 2, ^ a A; jsou obě sudá nebo lichá1. Uvažujme funkcionál V{K) — K2(x) dx a zabývejme se problémem najít takové jádro K e Svk, pro které tento funkcionál nabývá minimální hodnoty, tj. řešíme variační úlohu min V(K) za předpokladu K e S^fc. Řešení této úlohy se nazývají jádra s minimálním rozptylem, což jsou polynomy stupně k — 2 na intervalu [—1,1]. Tyto polynomy jsou sudé funkce pro k sudé a liché funkce pro k liché a mají A; — 2 různých kořenů v intervalu (—1,1). Obecný vztah pro jádra s minimálním rozptylem lze nalézt v [3]. Příklad 1.2. Jádra s minimálním rozptylem: Sí)2'- K(x) = 2J[-1,1 (x) Sl3- K(x) = 1,1] 0) = -f(l - 32ľ2)/[_M] (x) Poznámka 1.2. Jádra s minimálním rozptylem mají skoky v koncových bodech intervalu [—1,1], což negativně ovlivňuje hladkost výsledného odhadu. obvykle se používá pojem parita, tedy v a k mají stejnou paritu. 7 Tabulka 1.2: Jádra s minimálním rozptylem (MR) a optimální jádra (OPT) v k typ vzorec (na intervalu [—1,1]) P{K) T{K) 0 2 OPT 1(1-x2) 0,600 0,2000 0,3491 MR 1 2 0,500 0,3333 0,3701 0 4 OPT ±§(x2-l)(7x2-3) 1,250 -0,0476 0,6199 MR |(3 -hx2) 1,125 -0,0857 0,6432 1 3 OPT ^x(x2 - 1) 2,143 -0,4286 0,7477 MR ~~ 2X 1,500 -0,6000 0,8137 1 5 OPT -l§x(x2 - í)(9x2 -5) 11,93 0,1515 2,1675 MR ~g~ x (7 x 5) 9,375 0,2381 2,3278 2 4 OPT -^f (x2 - l)(5x2 - 1) 35,00 1,3333 6,6846 MR f (l-3x2) 22,50 1,7143 7,2622 2 6 OPT ^(x2-l)(77x4 - 58a;2 + 5) 381,6 -0,6293 27,162 MR ^(-45x4 + 42x2 - 5) 275,6 -0,9091 29,503 1.2 Optimální jádra Při vyšetřování statistických vlastností se setkáváme s následujícím funkcionálem T{K) xkK(x) dx 2v+l K2(x) dx k—u\ 2fe+l V(K) který lze zkráceně psát T(K) = (\pk(K)\2v+íV{K)(k-v>>)2/{2k+í). Jádra, pro která tento funk-cionál nabývá minimální hodnoty, se nazývají optimální jádra. Jde o polynomy stupně k, které mají k — 2 různých kořenů v intervalu (—1,1) a body —1,1 jsou rovněž kořeny těchto polynomů. Obecný vzorec pro tvar optimálních jader lze nalézt např. v [3]. Příklad 1.3. Optimální jádra: Sl3 S24 Kopt,0,2(X) = K1 -a;2)/[-l,l](:r) Kopt,l,3(x) = f1^2 ~~ 1)I[-1,1](X) KoptaA{x) = -1§ (x2 - l)(5x2 - l)/[_i,i](a;) Jádra -ř4T0pt,i,fc a Koptt2,k se používají pro odhad první a druhé derivace hustoty (viz kapitola 3.1) a z toho důvodu pro ně zavedeme dodatečné označení respektive K^2\ Přehled vybraných jader s minimálním rozptylem a optimálních jader je uveden v tabulce 1.2 Poznámka 1.3. Pro jádra s minimálním rozptylem a optimální jádra platí následující tvrzení, jehož důkaz je v [3]. Nechť K e ív+i,fc+i je jádro s minimálním rozptylem a Kopt^tk <= Svk je optimální jádro. Pak platí K'opt,vAx) = K(x)> ze [-1,1] Jako příklad lze uvést jádro Koptfit2(x) — 1(1 - x2) e S02 a K'opt02(x) — Klt3(x) — —|x e 5i3. 8 Poznámka 1.4. Optimální jádra jsou spojité funkce na celé reálné ose, což znamená, že jsou „hladší" než jádra s minimálním rozptylem. Odhadovaná funkce „zdědí" hladkost jádra a to znamená, že hladší jádra produkují hladší křivku. Nejčastěji používaným jádrem je Epanečnikovo jádro. Shrnutí Funkcionály závisející na jádře K 1 r-l pk{K) = J xkK{x)áx V{K) = J K2(x)dx T(K) = (\(3k(K)r+1V(K)^)^ 5vk = (jTfô) 2fe+1 Jádra s minimálním rozptylem minimalizují funkcionál V(K). Optimální jádra minimalizují funkcionál T(K). Podrobnější popis jader a vzorců s nimi souvisejících lze nalézt v toolboxu Matlabu. Výstupy z výukové jednotky Student • zná základní třídy jádrových funkcí, jejich vlastnosti a metody jejich konstrukce Dodatek Nechf K e Suk, ^ <= N, pro 7 > 0 položíme Jádra K a K7 nazýváme ekvivalentní. Nechf jsou dány funkce / a g, pro které platí /2(x)dx 0, přesněji rh{x, h) E Y~il[x-h,x+h](xi) i=l_ n E I[x-h,x+h}{%i) (2.3) Obr. 2.4 ilustruje aplikaci této metody na simulovaných datech příkladu 2.1. Uvedené metody Obrázek 2.4: Klouzavý průměr (červená, plná) pro simulovaná data z příkladu 2.1 s původní funkcí (modrá, čárkovaná) patří mezi nejjednodušší neparametrické vyhlazovací metody. Jádrové odhady lze považovat za zobecnění těchto metod. 13 Připomeňme zde základní myšlenku vyhlazování tak, jak ji formuloval R. Eubank v r. 1988: Jestliže předpokládáme, že m je hladká funkce, pak pozorování v bodech xi blízko bodu x obsahují informace o hodnotě m v bodě x. Bylo by tedy vhodné užít lokálních průměrů dat blízko bodu x, abychom získali odhad m(x). Obecně lze jádrové odhady regresní funkce m v bodě x definovat takto n m(x,h)=J2W>(x>h)Y» (2-4) i=l kde funkce Wi(x,h), i — 1,..., n, se nazývají váhy, nezávisí na hodnotách Yi, ale závisí na kladném čísle h, které se nazývá vyhlazovací parametr (nebo také šířka vyhlazovacího okna). Speciální, velmi užitečný typ vah, závisí na jádrové funkci K. Nechť K e Sofc, k je sudé číslo, položme Kh(-) — x^(tí)' Mezi nejznámější typy jádrových odhadů regresní funkce patří ([8]): 1. Nadarayovy-Watsonovy odhady (1964) n Y, Kh(x - xí)Yí mNW(x, h) = l-=^-, Y, Kh(x - Xi) i=l 2. Priestleyovy-Chaovy odhady (1972) 1 " mPC(x, h) = - V" Kh(x - xí)Yí, n '-^ i=l 3. lokálně lineární odhady (Stone 1977, Cleveland 1979) 1 ^2 {^(z, h) - š\(x, h)(x - Xi))Kh(x - xí)Yí i=l kde n ^—^ Š2(x, h)šo(x, h) — š±(x, h)2 1 " šr(x) — - y^(x - Xi)rKh(x - Xi), r — 0, 1,2, n '-^ i=l 4. Gasserovy-Múllerovy odhady (1979) Si n X r(x,h) = J2Y* J Kh(x-t)dt, mGM( kde So — 0, Si — (xí + Xj+i)/2, sn — 1. Tento odhad je konvolučním typem odhadu. Úmluva. Uvedené jádrové odhady budeme zapisovat ve tvaru n i=l kde index A značí příslušný typ odhadu NW, PC, LL, GM s danou váhovou funkcí. 14 V mnoha aplikacích je užitečný zejména Nadarayův-Watsonův odhad trmw- Popíšeme nyní jeho konstrukci a budeme ilustrovat vliv vyhlazovacího parametru na kvalitu odhadu. Pro daný bod xq, h < xq < 1 — h, jsou váhy Nadarayova-Watsonova odhadu dány vztahem Wi(x0,h)= Kh(£°-Xi) J2WÁ*o,h) = l. Obrázek 2.5 ilustruje konstrukci odhadu v bodě x0, který je založen na pěti pozorováních (xi, Yi), (x5,15) (černé křížky). Parabola reprezentuje Epanečnikovo jádro Kh a kroužky znázorňují hodnoty vah Wi(x0, h) — Kh(x0 — x i) / Y^=i Kh(xo — x i) pro i — 1,..., 5. Výsledný odhad regresní funkce fh v bodě x0 je označen hvězdičkou. Obrázek 2.5: Ilustrace Nadarayova-Watsonova odhadu v bodě xq OTÁZKA. Popište konstrukci Nadarayova-Watsonova odhadu, použijeme-li obdélníkové jádro místo Epanečnikova jádra. Vypočtěte váhy Wi(x0, h) pro odhad s obdélníkovým jádrem. (Odpověď viz kapitola 6.) Jádrový odhad není definován pro X)ľ=i ^h(x — x i) — 0. Jestliže nastane případ „0/0", pak klademe friNw(x, h) — 0. Omezíme se nyní na odhady funkce m v bodech plánu xí, í — 1,..., n. Pro malé hodnoty h je výraz 3-1 > 1 pro Xi ^ Xj, a tedy hodnota jádra v těchto bodech je rovna nule, a pro bod Xi dostáváme odhad mNW{x1, h) ->• K = Yi. To znamená, že při malé šířce vyhlazovacího okna (h —>• 0) odhad reprodukuje data (viz obr. 2.6(a)). Podobně pro velké hodnoty h je výraz 3-1 fe3'3 « 0, tedy pro všechny body plánu je hodnota jádrové funkce K \^Xl h^3 j « K(0) a dostaneme tak průměr dat n n Y.K{Q)Yj ^(O)E^ x „ E^(o) [)) 3=1 i=i Tedy velká šířka okna (h —>• 00) vede k přehlazení, a to k průměru dat (viz obr. 2.6(b)). Na obrázku 2.7(a) je znázorněn odhad s Epanečnikovým jádrem. Tento odhad se nejvíce blíží skutečné regresní funkci. Pokud jde o volbu vyhlazovacího parametru, je třeba si uvědomit, že konečné rozhodnutí o odhadované křivce je částečně subjektivní, neboí i asymptoticky optimální odhady obsahují poměrně značné ,/nnožství šumu" a to nechává prostor pro subjektivní posouzení. 15 (a) Podhlazený odhad, h = 0,02 (b) Přehlazený odhad, h = 0,4 Obrázek 2.6: Podhlazený a přehlazený odhad (červená, plná) regresní funkce (modrá, čárkovaná) z příkladu 2.1 (a) Odhad s Epanečnikovým jádrem a h = 0,08 (b) Odhad s obdélníkovým jádrem a h = 0,09 Obrázek 2.7: Odhady (červená, plná) regresní funkce (modrá, čárkovaná) z příkladu 2.1 s Epanečnikovým a obdélníkovým jádrem Poznámka 2.2. Intervaly spolehlivosti pro hodnotu regresní funkce to v bodě x jsou užitečné v mnoha aplikacích. Bodový interval spolehlivosti udává interval, v němž s pravděpodobností 1 — a leží hodnota funkce to v bodě x. Jsou definovány takto - IV(K)a2(x) _ IV(K)a2(x) m(x,h) -iíi-f y-^-,m{x,h) + u1-f Y-^- . kde Mi-f je (1 — a/2)-kvantil standardního normálního rozdělení a odhad rozptylu v bodě x je dán vztahem n d2(x) ^^2wi(x,h)(Yl-fh(x,h))2. i=l Ukázka intervalu spolehlivosti pro a — 0,05 je na obrázku 2.8. 3 Statistické vlastnosti jádrových odhadů Kvalitu jádrového odhadu lze lokálně popsat pomocí střední kvadratické chyby (MSE) odhadu to v bodě x, která je obecně dána vztahem MSEto(x, h) — E{rh{x, h) — to(x))2. 16 Obrázek 2.8: Interval spolehlivosti pro data z příkladu 2.1 při a — 0,05 (růžová, tečkovaná) se zobrazeným odhadem regresní funkce (červená, plná) a původní funkcí (modrá, čárkovaná) Upravíme tento vztah MSEm(x, h) — Efh2(x, h) — 2m(x)Em(x, h) + m2(x) — (£m(i, h) — m{x)Ý + Erh2{x, h) — (£m(i, h))2, (2-5) bias2 var což znamená, že střední kvadratická chyba může být vyjádřena jako součet rozptylu odhadu varm(i,/i) a čtverce vychýlení bias2 rh(x, h). Tento rozklad rozptyl-vychýlení usnadňuje analýzu vlastností odhadu. Všechny uvedené jádrové odhady regresní funkce fh^w, ínpc, t^ll, íncM jsou asymptoticky ekvivalentní (viz např. [8,14]). Z tohoto důvodu budeme dále podrobněji zabývat Priestleyovými-Chaovými odhady, které budeme psát bez uvedení označení PC, tedy: íh(x, h) a Wi(x, h). Připomeňme, že pro Priestleyovy-Chaovy odhady je váhová funkce tvaru Wi{x, h) = -kh{x - x,) = ±-k(^—^) . n nh \ h / Pro další výpočty budeme předpokládat: (i) Jádrová funkce K je sudou funkcí na intervalu [— 1,1], K e Sq^, (ii) vyhlazovací parametr h — h(n) je nenáhodnou posloupností kladných čísel splňující h —>• 0 a nh —> oo pro n —> oo, (iii) bod x, v němž počítáme odhad, splňuje nerovnost h < x < 1 — h pro všechna n > no, kde n0 je pevné, (iv) m e C2[0,1], (v) x, = i = 1,... ,n. Je zřejmé, že pro n —>• oo platí (jedná se o přibližný výpočet integrálu - viz Dodatek na str. 29)2 1™ 1™ f1 1 x t £m(i, /i) = - ^2Kh(x - xi)EYl = - ^(x - xí)m(xí) = / TK(-^— )m(*) dí + 0(n_1). n i=l n i=l ^° (2.6) Nechť u — 3Lff-, odtud dt — —h du, a tedy s využitím Taylorova rozvoje /x/h K(u)m(x — hu) du + 0(n_1) -{\-x)/h x/h r u2h2 i K(u) m(x) — uhm'(x) H--— m"(x) + o(h2) du + O^1). (2.7) -{\-x)/h L 2 -1 2 Symbolika O a o viz úvodní část na straně 5. 17 Podle výše uvedených předpokladů platí h < x < 1 — h, tedy i//i^ooa- (1 — x) / /i —>• — oo pro h —>• 0. Odtud, s využitím faktu, že nosičem funkce je interval [—1,1], plyne Em(x,h) = m(x) I K(u) du -hm'{x) I uK(u)du-\--m" (x) / u2K(u) du +o{h2) +0(n~1). =o =fo{K) Celkem dostaneme h2 bias íh(x, h) = —p2(K)m"(x) + o{h2) + O^1) Podobně pro rozptyl platí varm(x, h) — E(íh{x, h) — Efh{x, /i))2 /1 " 1 " = £ - V"^(aľ - Xi)Yi--^ Kh(x - xl)m(xl) \ n *—' n *—' , n *—' n v i=l i=l 1 / " = — E Kh(x - xt)(Yt -m{xi)) n vi=l z vlastností (2.2) plyne EKh(xi)Kh(xj)ei6j — 0 pro z 7^ j, tedy 1 " z—' i=l =^±*ra=£(/M2Tí)*+o<-'-,> Zde jsme opět použili přibližného výpočtu integrálu. Opět s využitím substituce u — a vztahu 0(n_1) — o((n/i)_1) můžeme pro n —>• 00 psát (7 2 /•! „2 /•! varm(i,/i) = — dw + o^n/i)"1) = — / K2(u) du + o(inh)-1) ^ —V(K) nh 7_i w/i 7-i n^ Tímto jsme dokázali následující větu o tvaru střední kvadratické chyby. Věta 2.1. Nectí jsou splněny předpoklady (i) — (iii), pak střední kvadratická chyba nabývá tvaru MSE m(x,h) = ^V(K) + h4p2(K)-(m''(x))2+o(h4 + (nh)-1). (2.8) nh v_4^__^ var bias2 Chyba MSE dává pouze lokální pohled na chybu odhadu, proto se častěji používá globální tvar chyby - AMISE - asymptotická střední integrální kvadratická chyba. AMISE je součástí střední integrální kvadratické chyby (MISE) a vztah mezi chybami MSE, MISE a AMISE je následující MISEm(-,/i) = j MSEm(x, h) dx = AMISE íh(-, h) + o(h4 + (nh)-1). Jo AMISE je tvaru 2 l4 AMISEm(-, h) = ^V(K) + —fá(K)V(m"), (2.9) AIV(h) AISB(Ti) 18 Obrázek 2.9: AMISE (růžová, plná) jako součet rozptylu AI V (červená, plná) a vychýlení AISB (modrá, čárkovaná) kde V(m") = JQ (m"(x)) dx a AIV značí asymptotický tvar rozptylu (asymptotic integrated variance) a AISB asymptotický tvar druhé mocniny vychýlení (asymptotic integrated square bias). Na obrázku 2.9 je znázorněn průběh AIV(h) a AISB(/i) a také výsledné chyby AMISE fh(-, h) Je vidět, že rozptyl AIV(h) nabývá velkých hodnot pro h malé, ale AISB(/i) klesá. Pro velké h je tomu naopak. Volba vyhlazovacího parametru je zřejmě klíčovým problémem jádrového vyhlazování. Naším cílem je minimalizovat AMISE fh(-, h), tzn. najít takovou hodnotu vyhlazovacího parametru h, pro kterou asymptotická střední integrální kvadratická chyba nabývá minimální hodnoty, a tedy odhad bude nejlepší ve smyslu AMISE. Užijeme metody matematické analýzy a spočítáme derivaci dAMISEm(-,/i) _ cr2V(K) dh nh2 položíme ji rovnu nule a vyjádříme h h3(322(K)V(m"), o2V(K) opt'°'2~ n(32(K)V(m")- [ Uj Poznámka 2.3. Tento výpočet vede k nalezení minima AMISE, protože platí d2 AMISE dh2 > 0. Vztah (2.10) má pouze teoretický charakter, protože hodnota /iOpt,0,2 závisí na neznámých veličinách a2 a m"(x), a tedy není užitečná pro praktické účely. Abychom odhadli optimální hodnotu vyhlazovacího parametru, musíme použít metody, které jsou založeny na datech (data-driven methods). Nejznámější z těchto metod bude uvedena v dalším odstavci. Vztah (2.10) pro optimální šířku vyhlazovacího okna ukazuje, že řád konvergence optimální šířky vyhlazovacího okna /iOpt,0,2 závisí na řádu jádra K, tedy pro jádra řádu 2 je 0(n_1/5). Do-sadíme-li (2.10) do vztahu (2.9) pro AMISE, dostaneme AMISE( V,,0,2) = Č.VIK) (,,^L„,) "* + i Um'^r *W = (a2)4/5 (V(K))4/5n-4/5 (Pl (K))1/b (V(m")Y'b + I((J2)4/5(-l/(^))4/5n-4/5(/32(^))l/5(-l/(m^))l/5 ^(ffV(ií))4/5'fl2'^"~'Ml/5 tj. AMISEm(-, /iopt.0,2) = 0(n-V5). l(v2V(K))4/5(p2(K)V(m"))1/5n-^, (2.11) 19 Poznámka 2.4. Jestliže jádro K náleží do třídy Sok, pak AMISE je tvaru AMISE m(-,/i) nh V(K) (k\)2 h2kl3l(K)V(m(k)) (2.12) a pro optimální vyhlazovací parametr h platí 2fe+1 a2V{K){k\f opt'°'k 2knPl{K)V{m(k)y (2.13) kde V(m^) — {rrSk\x)Ý dx, podrobněji např. [7]. Nyní uvedeme důležité lemma, které ukazuje zajímavou vlastnost vyhlazovacího parametru. Lemma 2.1. Pro /iopt,o,fe platí AIV(/iopt,o,fc) = 2kAISB(hopt^k). Důkaz. Viz cvičení. □ Lze ukázat, že pro jádra K e Sofc je AMISE m(-, h) — 0(n~ 2k+1). To znamená, že s rostoucím k se zvyšuje asymptotická rychlost konvergence. Ale není zcela jasné, zda tato zvýšená rychlost konvergence vede již k zlepšení pro konečné rozsahy výběrů, neboť ostatní veličiny se rovněž mění s k. Nevýhodou jader vyšších řádů je fakt, že pro tato jádra je optimální šířka okna větší, což může mít negativní dopad na hraniční efekty [9]. Na druhé straně, chování jádrových odhadů s jádry vyšších řádů je méně citlivé na volbu šířky okna, není-li určena zcela optimálně, neboť křivka AMISE m(-, h) je plošší. Poznámka 2.5. Vyšetřování kvality odhadu obvykle probíhá za předpokladu, že pracujeme s vnitřními body intervalu [0,1]. V hraničních oblastech, tj. v intervalech [0, h) U (1 — h, í], je kvalita odhadu ovlivněna negativně skutečností, že jádro K zde nesplňuje momentové podmínky (1.1). To je způsobeno tím, že blízko krajních bodů nosič jádra K zasahuje do oblasti, kde nejsou žádná data, což zhoršuje odhad - viz obr. 2.10. Hraniční efekty jsou také patrné na obrázcích 2.7(a) Obrázek 2.10: Hraniční efekt a 2.14, zejména u pravého okraje intervalu. Problém okrajových efektů lze překonat např. použitím hraničních jader (viz [9]) nebo reflexní metodou (viz [3]). Ukázkový příklad 2.2. Uvažujme simulovaná data generovaná regresní funkcí m(x) — sin2 ttx na intervalu x e [0,1] s chybami Ei ~ N(0; 0,252). Vypočítejme hodnotu optimálního vyhlazovacího parametru pro odhad s jádrem řádu 2. 20 Podle vztahu (2.10) potřebujeme spočítat výraz V(m"). m{x) — sin2 7ix m'{x) — 2 sin ttx cos ttxtt — 7rsin27ra; m"(x) — tt cos 2irx2ir — 27T2cos27ra; V(m") = / [m"(x)]2dx = í [27r2cos27rx]2dx Jo Jo 4tt / cos 2irx dx — 4-tt Jo Jo 2TT4. 1 + cos 4:irx ■ dx Výpočet /iopt,0,2 pro • Epanečnikovo jádro: V(K) — 3/5, ^(K) — 1/5, ^23/5 "opf.0,2 - n(!/5)227r4 - 27r4n' obdélníkové jádro: = l/2,ft(iř) = 1/3, ^5 _ a21/2 _ 9a2 n(l/3)227T4 ~ 47r4n' 'opt,0,2 Odhady s optimálním vyhlazovacím parametrem pro soubor o velikosti 50 hodnot jsou na obrázku 2.11. (Data jsou v tabulce 7.2.) Vidíme, že odhad s „hladším" Epanečnikovým jádrem generuje „hladší" křivku. 0.2 0.3 0.2 0.3 (a) Epanečnikovo jádro, hopt,o,2 = 0,1573 (b) Obdélníkové jádro, ftopt,o,2 = 0,1236 Obrázek 2.11: Odhad regresní funkce z ukázkového příkladu 2.2, odhad (červená, plná) a původní funkce (modrá, čárkovaná) 4 Volba jádra Volba jádra není z asymptotického hlediska podstatná, jak je zřejmé z faktu (2.11). Je vhodné zvolit optimální jádro, které minimalizuje funkcionál T(K), nebof tato jádra jsou spojitá na R a odhadovaná regresní funkce tak „zdědí" hladkost jádra. Vhodná jsou jádra třídy S02 a S04, lze je vybrat z tabulek 8.1 pro Sofc- 21 5 Volba vyhlazovacího parametru 5.1 Metoda křížového ověřování Jednou z nejrozšířenějších a nejpoužívanějších metod pro určení optimální hodnoty parametru h je metoda křížového ověřování (cross-validation method). Tato metoda je založena na odhadu regresní funkce (2.4), v němž vynecháme í-té pozorování: n fh-i{xi,h) — wÁxi, h)Yt, i — í,...,n. i=i Funkce křížového ověřování je definována takto 1 ™ CV(/i) = - Y,{™-í{?í, h) - Yi)2 (2-14) i=l a odhadem optimální hodnoty vyhlazovacího parametru je bod, v němž nastává minimum této funkce, tj. hopt,o,k = /icv = arg min CV(h). Hledáme tedy minimum na intervalu Hn — [akn~1^2k+1\bi!Ti~1^2k+1^], jehož tvar plyne ze vztahu (2.13), přičemž cifc, bk jsou konstanty (0 < ak < bk < oo), které ovšem neznáme. A proto pro ekvidistantní body plánu byl na základě zkušeností doporučen interval [^,2]. Poznámka 2.6. Někdy se místo chyby MISE používá průměrná střední kvadratická chyba AMSE (average mean square error) 1 ™ AMSEm(-,/i) = -E^2(m(xi,h) - m(xt)ý n i=l Využívá se zejména v případech, kdy není vhodné použít numerické integrování související s chybou, která se vyskytuje v MISE. Věta 2.2. Pro střední hodnotu funkce CV(h) platí 1 ™ ECV(h) = -E V(m(xi,/i) - m(xt))2+a2. n ^-^ i=l AMSE Důkaz. Funkci křížového ověřování lze rozepsat 1 ™ 2 CV(/i) = - ^(rn_j(xj, h) - m(xl) - (Yl - m(xj))) ^ n \ n 1 n = - y^(m-i(xi, h) - m(xl)) - 2- V"(m_j(x4, h) - m(xl))el + - V" e2 i—l i—l i—l ^ n ^ n n \ n — - ^2(m-i(xi, h) - m{xi))--£í Wt{x1,h)Yl - m(xl)j + - ^ e; n *—' n ů—' v-'—' / n i=l i=l i=l i=l Střední hodnota E CV(h) je rovna součtu tří veličin. Předpokládejme, že ín_j(xj, h) « fh(xi, h), pak první ze sčítanců je roven přímo AMSE fh(-, h): 1 ™ -E^2(m-i(xi,h) -m(xt))2 = AMSEm(-,/i). n i=l 22 Dále víme, že Yt — m(xi) + en, a tedy pro druhý sčítanec platí: 2 n n i=l £=1 £^i 2 n n =--E^2sl(^2We(xl,h)(m(xe) + et) - m(xt)J 71 i=l £=\ £^i Stejně jako pro druhý sčítanec, i pro třetí sčítanec využijeme vlastnosti (2.2): n n —^ i=l □ Tento výsledek znamená, že minimalizace EGV(h) odpovídá minimalizaci AMSE. Jestliže tedy předpokládáme, že minimum CV(/i) je blízko minima E CV(h), pak tato minimalizace dává dobrou volbu vyhlazovacího parametru - viz ilustrace na obr. 2.12. ......T" Obrázek 2.12: Porovnání minima AMSE (modrá, čárkovaná) a minima funkce křížového ověřování CV (červená, plná) pro simulovaná data z ukázkového příkladu Příklad 2.3. Použijeme metodu křížového ověřování pro nalezení vyhlazovacího parametru pro data z příkladu 2.1. Při použití Epanečnikova jádra získáme vyhlazovací parametr hcv — 0,1158. Na obrázku 2.13 je zobrazen odhad regresní funkce s tímto parametrem. Kromě metody křížového ověřování se také pro odhad optimálního vyhlazovacího parametru používají metody založené na ASE (average square error), metody plug-in, metody odvozené z Fourierovy transformace a bootstrapové metody (podrobněji např. [3, 6]). 23 Obrázek 2.13: Simulovaná data (x) s jádrovým odhadem regresní funkce (/icv — 0,1158) (červená, plná) a původní funkcí (modrá, čárkovaná) 6 Automatická procedura Z dříve uvedených odhadů chyb je zřejmé, že kvalita jádrového odhadu závisí na šířce okna h, na jádře K a na řádu jádra k, což je číslo, které odpovídá předpokládanému počtu derivací v odhadovaném modelu. Je zřejmé, že všechny tyto tři veličiny se vzájemně ovlivňují, a proto je třeba zabývat se jejich volbou současně. Pro simultánní volbu jádra, optimálního vyhlazovacího parametru a řádu jádra byla navržena automatická procedura (viz [3]), která odhadne všechny parametry tak, aby byla minimalizována AMISE. Procedura byla původně odvozena pro odhad hustoty pravděpodobnosti ([4]), ale lze ji aplikovat i pro odhad regresní funkce. Uvedeme zde její zjednodušenou verzi. Podle vztahů (2.12) a (2.13) je AMISE tvaru AMISE m(-, h) = ^rV{K) + -l-h2k^K)V(mw) nh {k}.)1 a hopt^k tvaru °Pt'°'fe 2knV{m(k)) ' Ze vztahu pro /iOpt,0,fe vypočteme V(m^) a dosadíme do vztahu pro AMISE. Dále použijeme vztahy Sok = (|^) ^ T{K) = (í3k(K)Vk(K)) ™& T(K) = 5fkft{K) a dostaneme vyjádření AMISEm(.,/loptj0,fc) = gyfc+1)^T(g); ve kterém jsou parametry K, h ak separovány, což umožňuje vybrat tyto parametry simultánně. Právě tento vztah je základem automatické procedury. Položme a množinu vhodných řádů k označme J(fco) = {2j,j = 0,...,[^]}, kde [z] značí celou část čísla z. Procedura pak probíhá v pěti krocích: 24 Procedura 1. Pro k G I(ko) najděte optimální jádro Koptto,k <= Sok, které je dáno tabulkou 8.1 a vypočtěte kanonický faktor Sok- 2. Pro k e I(k0) a Kopt0,k <= S'ofe najděte optimální vyhlazovací parametr hopt0,k- 3. Pro & e /(/ío) 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. 4. Vypočtěte optimální hodnotu řádu k, které minimalizuje funkcionál L(k). 5. Použijte parametry z předchozích kroků k získání optimálního jádrového odhadu regresní funkce, tj. n i=l Příklad 2.4. 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 kroku 2 jsme použili metodu křížového ověřování pro nalezení optimálního vyhlazovacího parametru hoptflj:- k Kopt,0,k ^Ofe hOpt,0,k L(k) 2 -i(-2- -1) 1,7188 0,1158 0,0648 4 15/x2 32 vx l)(7x2-3) 2,0165 0,2446 0,0575 6 105 12 256 v - l)(33x4 - 30x2 + 5) 2,0834 0,3575 0,0574 8 315 / 2 4096 v - l)(715x6 - -lOOlx4 + 385x2 - - 35) 2,1021 0,4543 0,0592 Z tabulky vidíme, že optimální řád jádra je k — 6. Výsledný odhad je uveden na obrázku 2.14. "*........................> , " * , i-7---v. „ < 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obrázek 2.14: Simulovaná data (x) s jádrovým odhadem regresní funkce při použití procedury (červená, plná) a skutečnou funkcí (modrá, čárkovaná) společně s 95% intervalem spolehlivosti (růžová, tečkovaná) 7 Aplikace na reálná data Nyní se vrátíme k motivačnímu příkladu. Datový soubor obsahuje měření úrovně hladiny vody v Hurónském jezeře. Měření byla prováděna ročně, v letech 1875 až 1972, tedy naměřené hodnoty jsou ekvidistantní. Data jsou shrnuta v tabulce 7.8 a na obrázku 2.15(a). 25 Vyhlazovací parametry jsme odhadli pomocí metody křížového ověřování (za použití Epa-nečnikova jádra) a také pomocí automatické procedury. Hodnoty vyhlazovacích parametrů jsou následující: /icv = 0,0204 (^,0,2), h = 0,1525 (^,0,12). Výsledná křivka fh(x) je odhadem funkce m(x) popisující průběh úrovně hladiny v letech 1875 až 1972. Odhady na obrázku 2.15 ukazují, že metoda křížového ověřování spíše podhlazuje. 1870 1880 1890 1900 1910 1920 1930 1940 1950 1960 1970 (a) Data 1900 1910 1920 1930 1940 1950 1960 1970 (b) hcv = 0,0204, Kopt^2 1900 1910 1920 1930 1940 1950 1960 1970 (c) h = 0,1525, Kopt,o,12 Obrázek 2.15: Úroveň hladiny Hurónského jezera - data a odhadnuté regresní funkce, na ose x jsou roky a na ose y je hladina vody ve stopách (snížená o 570 stop - viz poznámka u dat na str. 90) Druhým datovým souborem jsou měření axiální délky krystalů ledu. Měření byla prováděna v místnosti s konstantní teplotou — 5°C v časových intervalech 50 až 180 vteřin po přinesení krystalu do místnosti. V tomto případě nejde o ekvidistantní data, protože hodnoty se liší o pět či deset vteřin. Data jsou v tabulce 7.9 a na obrázku 2.16(a). Chceme odhadnout funkci, která vyjadřuje závislost axiální délky krystalů na čase. Odhady vyhlazovacích parametrů podle metody křížového ověřování a pomocí automatické procedury: /icv = 0,1865 (tfopt.0,2), h = 0,8826 (^opt,o,io)- Výsledné odhady regresní funkce jsou na obrázku 2.16. Je vidět, že metoda křížového ověřování dává spíš podhlazený odhad. Na druhou stranu, odhad pomocí procedury se může zdát již přehlazený. 26 (a) Data (b) hcv = 0,1865, Kopt,0,2 (c) h = 0,8826, Kopt,o,io Obrázek 2.16: Axiální délka krystalů ledu - data a odhadnuté regresní funkce, na ose x je vynesen čas ve vteřinách a na ose y délka krystalu v mikrometrech 27 Shrnutí Odhad regresní funkce m(x) v bodě x je tvaru i(x,h)='52Wi(x,h)Yi. i=l Wi(x, h) závisí na K a h. Vychýlení (bias) a rozptyl (var) odhadu s jádrem řádu 2 jsou h2 _ a2 biasm(x, h) ~ —j32{K)m (x), v&rm(x,h) « —^V(K). Asymptotická střední kvadratická chyba jádrového odhadu regresní funkce s jádrem řádu 2 AMISEto(-, h) = ^V(K) + ^h4p2(K)V(m"). Optimální vyhlazovací parametr vzhledem k AMISE pro k — 2 je tvaru 'AMISE — —jT, nß2{K)V{m") s řádem konvergence n-1/5. Metoda křížového ověřování pro odhad optimálního vyhlazovacího parametru 1 ™ CV(h) — - ^(rn_j(xj,/i) - Yi) => hCy — argminCV(/i). n i=l Automatická procedura pro simultánní volbu optimálního jádra, jeho řádu a vyhlazovacího parametru je dostupná v toolboxu Matlabu. Výstupy z výukové jednotky Student • zná základní typy jádrových odhadů regresní funkce a jejích derivací • je schopen analyzovat statistické vlastnosti odhadů • má přehled o metodách pro volbu vyhlazovacího parametru • se seznámil s automatickou procedurou pro simultánní volbu vyhlazovacího parametru, jádra a jeho řádu • je schopen analyzovat daný soubor dat a aplikovat uvedenou proceduru na tento soubor • je schopen použít příslušný toolbox v Matlabu a zkonstruovat odhad regresní funkce 28 Dodatek Vypočet integrálu, ti — ^, i — 1,..., n, n-l „tl+1 G(t) dt = J2 dí s využitím Taylorova rozvoje funkce G(ť) n-l „tl+1 = / (G(U+i) + (t- ti+1)G'(ti+1) + o(l)) dí n-l n-l »ti+1 = 53 G(íi+i)(íi+i -U)+ ^2 (t - ti+1)G'(ti+1) dí + o(l) - \ E G(^) + E 2 G'^) + °(!) n -. n — 1 —E^ + ^E^+iH^1)- i=l i=0 Za předpokladu |G"(íí+i)| < M, pak platí -i , n / G(t)dt = - ľGfcjlOr1). Jo Cvičeni 1. Odvoďte tvar Priestleyova-Chaova odhadu regresní funkce v bodě x pro obdélníkové jádro. 2. Pro odhad fh^w dokažte, že „množství" vyhlazení pomocí jádra K s vyhlazovacím parametrem h je stejné, jako „množství" vyhlazení jádrem Kg s parametrem h* — h/S,t). ŕh{x, h, K) — m{x, h*, Kg). 3. Uvažujte funkci m(x) — sin2 irx a regresní model Yi — m(x) + Si, í — í,dots, 100, kde Si ~ N(0; 0,25). Vypočtěte hodnotu /iopt,o,4 pro odhad s jádrem K(x) — |§(x2 — l)(7x2 — 3). Pomůcka: V(m^) = 32tt8. 4. Vypočtěte hcv pro data z ukázkového příkladu 2.2 a porovnejte s hopt. 5. Dokažte vztah (2.11) pro obecné k, tj. že platí AMISEÍ/^^n-2*/^1)^ V / V / (2fc)5fe+T 6. Dokažte vztah z lemmatu 2.1. 7. Aplikujte metodu křížového ověřování a automatickou proceduru na simulovaná data z ukázkového příkladu 2.2. 29 Kapitola 3 Jádrové odhady hustoty 1 Motivace Velká část pozorování a měření prováděných v biologii i jiných vědách má výsledek vyjádřený reálným číslem. Tato čísla jsou hodnotami nějaké reálné náhodné veličiny. Jako příklad můžeme vzít měření obsahu cholesterolu v krevní plazmě pacientů. Jde o soubor 371 měření - hodnoty jsou na obrázku 3.1 - která byla provedena u pacientů, kteří si stěžovali na bolest na hrudi. Data pocházejí z rozsáhlé studie, v níž autoři zkoumali vliv lipidů a jiných látek na nemoci srdce [10,11]. -+ +i- +-m +■ + + + + - 100 150 200 250 300 350 400 Obrázek 3.1: Měření cholesterolu v krevní plazmě 371 pacientů Chceme zjistit, jak často se v populaci vyskytuje dané množství cholesterolu v krvi? Případně, jaká je pravděpodobnost, že pacient bude mít zvýšený obsah cholesterolu? Rozdělení pravděpodobnosti je popsáno reálnou funkcí, která se nazývá hustota pravděpodobnosti a značí se /. Hustota pravděpodobnosti je základním pojmem ve statistice [1,2]. Funkce f(x) hustotou spojité náhodné veličiny, jestliže platí • f(x) je nezápornou funkcí pro všechna x, • J f (x) dx — 1. 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 30 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. • Histogram je schodovitá funkce, ale přitom předpokládáme, že neznámá hustota je spojitá. Příklad 3.1. Mějme dán datový soubor generovaných ze směsi dvou normálních rozdělení N(0; 0,25) a N(2; 0,25) s hustotou 1 2 1 0-2)2 f(x) = 0,3 ;„ , e~5Ť5 +0,7 ,„ _ e — který má rozsah n — 100. (Data jsou v tabulce 7.3.) Z obrázku 3.2 je patrné, že histogram nevystihuje korektně hustotu pravděpodobnosti dat. -0.5 0 0.5 1 1.5 2 2.5 (a) 7 intervalů o šířce h = 0,55 -0.5 0 0.5 1 1.5 2 2.5 (b) 13 intervalů o šířce h = 0,29 Obrázek 3.2: 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 e M je definovaný vztahem [14]. 1 n Y 1 n i=l (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ě. Na obrázku 3.3(a) jsou čárkovaně zobrazena jednotlivá Epančnikova jádra v bodech Xi (kroužky) a plnou čarou je pak zobrazen odhad hustoty. OTÁZKA. Popište konstrukci odhadu s obdélníkovým jádrem. Jak bude tento odhad vypadat? 31 Nyní uvedeme ještě vztah pro jádrový odhad ^-té derivace hustoty. Budeme předpokládat, íei) • 0 a n/i —>• oo pro n —>• oo. Pak platí mise/(,/*) = TO + ^^(jfiyf/Wj + ofr + H-1), kdeV(f^) = j(/(fe)(x))2da;. Důkaz. Nejprve vypočteme střední hodnotu 1 7./2ľ - y Ef(x,h) = (Kh*f)(x) = j Kh(x-y)f(y)dy = j -K(^-)j{y) dy = / K(z)f(x - hz) dz dále použijeme Taylorův rozvoj: f(x — hz) — f(x) — f'(x)hz + • • • + ^ f(k\x)hkzk + o{hk) = y K{z) [f(x) - f(x)hz +■■■ + L$LfW(x)hkzk + o(hk)] dz = /(x) / K(z)dz-f'(x)h í zK(z)dz+--- + ^^f(-k\x)hk í zk K (z) dz +o(hk) =1 =0 =0k(k) f{x) + (-l)krk\x)hkMK) + 0{hky Tedy vychýlení odhadu je tvaru bias/(•,/*) = tf/fo/i) - /(*) = 1 ' ^ 1 V/?fc(g) + o(hk), =o(l) 33 a tedy Ef(x,h) = f(x) + o(l). Nyní dokážeme vztah pro rozptyl. Víme, že 1 var/(x, h) = -((K2 * /)(*) - * ff{x)) {f{x)+o{l)Y a dále počítáme ^ I *T2(z) /(* - /iz) dz - + o(l))2 =/(z)+o(l) | K2(z)(f(x) + o(l))dz-±-(f(x)+o(l))2 Tedy o((«fe)-i) ^ í K2(z)dz + o((nh)-1). nh . MSE/(x, /i) = / A"2(*)dz + o((n/0_1) + i- ,, V ' hk/3k(K) + o(hk) nh J \ kl I{X) ľ K2(z) dz + o{{nh)^) + /(fc>(x)£^(ií) nh J w ' ' ' ' ■> y '(/d)2' jfc! =o(>2fe) a pak využijeme faktu, že j f(x) dx — 1 MISE/(-,/i) = / MS~Ef(x,h)dx ™ + ^^AW(/(fe)) + o (h2k + (nh)-1) □ Důsledek. Nectí h —>• 0, n/i —>• oo pro n —>• oo, paA: / je konzistentním odhadem /, ŕ/. E1/ —>• / fl var / —>• 0. 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) = + J^h2k(32(K)V(f^). (3.3) 34 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) = AISB(/0 = J^h2kfá(K)V(f^), tedy AMISE/(■, h) = AIV(/i) + AISB(/i). Užitím vztahů T(K) = (V(K)k pk(K))V(2k+1) a S^k+1 — jtj^j pro K e S0k lze AMISE zapsat ve tvaru AMISE/(,A)=T(/0(^ + ^p). 0.4) Důkaz viz cvičení. Odtud je zřejmé, že vyhlazovací parametr, pro nějž AMISE nabývá minimální hodnoty, je dán vztahem h2k+l _ °0k \K-) /o c\ noPt,o,k- 2knV(fWy 1 ' tj. Äopt.o.fc = Oín-V^+D). Vypočtěme hodnotu AMISE při dosazení optimálního parametru hopt0j„: AMISE/(, W) = T(A-)y(/(fc))V(2fc+i)n-2fc/(2fc+i) (2fc(^2+;(2fc+1), (3.6) tj. AMISE/(■,/iopt,o,fc) = 0(n-2fe/(2fe+D). I v tomto případě, podobně jako u odhadhu regresní funkce, platí vztah mezi asymptotickým rozptylem ATV(h) a vychýlením AISB(/i): AIV(/iopt,0,fe) = 2fcAISB(/iopt,0,fc)- (3.7) Nyní uvedeme zajímavou vlastnost vyhlazovacího parametru. Poznámka 3.1. Nechť K e So2- Pak optimální hodnota vyhlazovacího parametru je i)5 h5 - 02 opt'0'2 nV(f")' Počítejme derivace AMISE /(■, h) dané rovnicí (3.3) pro k — 2 d»AMISE/(,,)^j|ffi + 6)ife2(A,)v,(r) Řešením rovnice d3AMISE /(■, /i)/ d/i3 = 0 je V(K) 6? 02 n/32(^)l/(/") nV(f") °^0'2' tj. /lopt,0,2 také realizuje minimum d2AMISE /(■, /i)/ d/i2. Lze ukázat, že d2AMISE(/(-,/i)) d/i2 / _ 2fc-2 , 35 0.1 0 h_2 h_4 h_6 Obrázek 3.4: AMISE/(•, h) pro jádra vyšších řádů s vyznačenými minimálními hodnotami pro jádra řádu 2,4, 6 a to znamená, že pro jádra vyšších řádů je minimum AMISE f(-,h) plošší a tedy volba h blízká optimální hodnotě /iOpt,0,fc nevede k velkému růstu AMISE /(•, h). Na obrázku 3.4 jsou zobrazeny body minima funkce AMISE f(-,h) pro hustotu normálního rozdělení A^(0; 1) se sto prvky. Vztah pro optimální hodnotu vyhlazovacího parametru poskytuje informaci, že asymptoticky )eh — 0(n_1/(2fe+1)). Ale vztah má pouze teoretický charakter, protože optimální parametr závisí na neznámé hustotě /. Je zde tedy opět problém s volbou tohoto parametru. Metodám pro odhad vyhlazovacího parametru je věnován odstavec 5. 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 kde ak, bk jsou konstanty, 0 < ak < bk < oo. Ukázkový příklad 3.2. Máme k dispozici data, která pocházejí z rozdělení s hustotou f(x) — ||(1 — x2)3 pro x G [—1,1]. Vypočítejme hodnotu optimálního vyhlazovacího parametru pro odhad s jádrem řádu 2. Podle vztahu (3.5) potřebujeme spočítat výraz V(f"). /(x) = §(l-x2)3 nn = Jjf"(x)]2dx = 35 Výpočet /i0pt,o,2 pro Epanečnikovo jádro: V(K) — 3/5, ^(K) — 1/5, tedy 5q2 — ^(k) ~ ^ 5 _ 15(2!)2 _ _3. opt'0'2" 2-2-n-35 " In Tedy pro soubor 50 hodnot bude /iOpt,0,2 = 0,8441 • 50-1/5 = 0,3860. Odhad s optimálním vyhlazovacím parametrem pro tento datový soubor (viz tabulku 7.4) je na obrázku 3.5. 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. i=l 36 Obrázek 3.5: Odhad hustoty z ukázkového příkladu 3.2, odhad (červená, plná) a původní funkce (modrá, čárkovaná) při použití Epanečnikova jádra a /iOpt,0,2 — 0,3860 Předpokládejme nyní, že platí Q < v < k — 2, vak mají stejnou paritu, linin^oo h — 0, lim^oo nh2v+1 = oo, / e Cfeo (A: < k0) a F(/(fe)) = J(f^\x)fdx < oo. Pak lze ukázat, že asymptotická střední kvadratická chyba AMISE f^\-, h) je tvaru AMISE/»(,/*) = ^-g^ + J^h^-^ftiK^VVM). Důkaz je založen na použití vhodného Taylorova rozvoje hustoty /, podobně jako u důkazu tvaru AMISE u odhadu hustoty [3]. Optimální hodnota vyhlazovacího parametru je dána vztahem 2fe+1 ^+1(2^+l)(fcQ2 kdp ,2fc+i V(K^) noPt,u,k- 2n(k-v)V(f^) ' vk ~ p&KMy 1 ' Tento vzorec umožňuje výpočet optimálního vyhlazovacího parametru pro pomocí /iOpt,0,fc a hopt,i,k- Předpokládejme nejdříve, že ľ a k jsou sudá čísla. Pak opt,u,k _ ({2v+l)k\1/[2k+1) 6; vk 'ioPt,o,k \ k — v ) Sok Pro vak lichá platí h_1 je standardní normální kvantilová funkce a číslo A[3„/4], respektive X[„/4], je horní, respektive dolní výběrový kvartil. Je vhodné volit mm{a s d iqr} ■ Poznámka 3.4. Pokud za jádro K zvolíme Gaussovo jádro (k — 2) pak dostaneme jednoduchý vztah [11,12] / 4 \i/5 /iref = y—J cr. (3.15) Příklad 3.3. Použijme odhad vyhlazovacího parametru pro data z příkladu 3.1 metodou referenční hustoty. Pro Epanečnikovo jádro, které je řádu k — 2, se vztah (3.12) zjednoduší na tvar /8V5F\ . ^ _1/5 1/5 V 3 " Dále odhadneme směrodatnou odchylku: agjj — 1,0325, ctiqr — 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.7 je vykreslen odhad hustoty s tímto parametrem. Obrázek 3.7: Odhad hustoty s /iref — 0,9639, odhad (červená, plná), původní funkce (modrá, čárkovaná) OTÁZKA. Jak hodnotu bude mít odhad vyhlazovacího parametru pomocí metody referenční hustoty pro hustotu z ukázkového příkladu 3.2? 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. 39 Věta 3.3 (Terrell 1990). Mezi všemi hustotami f s nosičem [—1,1] má hustota rozdělení Beta(k+2, k + 2) (_(2fc+3)!_(l-x2)k+1 \x\< 1 10 jinak, nejmenší hodnotu integrálu (/(fe) (x))2 dx. Lze ukázat, že platí 2. Pro r > 0, J(rg(k\rx)^2 dx — r2k+1 j(g(k^ (x))2 dx pro každou hustotu, pro kterou integrál existuje. 3. Jestliže hustota g má rozptyl a2, pak hustota ^fg{^-x) má rozptyl a2. Jestliže / je neznámá hustota s rozptylem a2 a gk je hustota rozdělení Beta(/j + 2, k + 2), pro kterou je V(g^) minimální, pak využitím faktu 1-3 lze ukázat, že ({k\)2Y1+i ° (v((k)^ Hodnotu cr lze odhadnout pomocí dříve uvedených vztahů aa-k — 2fc1|_5. Hodnotu V{gkk^) lze vypočítat pomocí speciálních ortogonálních polynomů [5]: V(o^) ľ (o^(x)Ýdx _(2fc + 3)!(2fc + 2)!_ (3 17) V(9k '-J^k dx~ 22><+2(2k+l)(2k + 5)(k + l)\2- (ÓA/) Použijeme-li poslední vyjádření (3.17), dostaneme horní hranici pro vyhlazovací parametry hopt,o,k < hMS = dn-1/{2k+1)bkl pncemz , _ /KJ-^(22k+2V{K){2k + \){2k + 5)(fc + l)2(fc!)2\ fc_V + 1 P2{K){2k + mik + 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 3.4. Určeme hodnotu Iims Pro odhad hustoty s jádrem řádu k — 2, tj. K e So 2-Podle věty 3.3 je hustota g% tvaru 35 92{x) = —(1- x2)3, a; e [-1,1], a dále z vlastností funkce g? a ze vztahu (3.17) plyne 1 ŕ 2 x2g2(x) dx — — a y (g"(x)) dx = 35. Pak pro Koptfl2 - -1/5243 \ VŠ 40 Příklad 3.5. Pro data z příkladu 3.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 protože a — 1,0325 a n — 100. Výsledný odhad je vidět na obr. 3.8. Obrázek 3.8: Odhad hustoty s /ims — 1,0409, odhad (červená, plná), původní funkce (modrá, čárkovaná) Poznámka 3.5. Hodnota /ims 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 — [h^, /ims]/ kde ha je nejmenší vzdálenost mezi po sobě jdoucími body Xi, í — 1,..., n. OTÁZKA. Jak hodnotu bude mít odhad vyhlazovacího parametru pomocí metody maximálního vyhlazení pro hustotu z ukázkového příkladu 3.2? 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=£ J (f(x,h)-f(x)Yáx = E í f2(x,h)dx -2E í f(x,h)f(x)dx+ I f2(x)dx. MISE/(-,/i) - J f2(x)dx = EÍ^J f2(x,h)dx-2 J f(x,h)f(x)dx Definujme funkci křížového ověřování ľ 2 n CV(h) = j f2(x, h)dx--J2 f-ÁXi, h), (3.18) i=l kde f-i (Xi, h) je odhad v bodě Xi bez použití tohoto bodu. Věta 3.4. Platí #CV(/0 = MISE/(-,/0-y f(x)dx, tj. CV(h) je nevychýleným odhadem El / f (x,h)áx — 2 \ f(x,h)f(x)dx 41 Důkaz. Spočítejme střední hodnotu funkce CV(/i), tj. ľ 2 " ECV(h)=E / f(x,h)dx--Ej2f^(X,h). J n 2=1 Dále upravme druhý člen tohoto vyjádření: 1 n 1 n 1 n E- Y, f-i(Xi, h) = E-Y-7 E - xj) 2 = 1 2 — 1 J—1 1 ' 2=1j = l Kh(x -y) f (y)J (x) dxdy = E J f (x, h) f (x) dx E f (x,h) Odtud ECV(h)=E J f2(x,h)dx-2E J f(x,h)f(x)dx. Odhad hoptfl,k je dán vztahem □ /icv — arg min CV(/i). h£Hn Odtud plyne, že CV(/i) + J f2 (x) dx je pro každé /i nevychýleným odhadem MISE f(-,h). Protože J f2 (x) dx nezávisí na h, minimalizace E CV(h) odpovídá minimalizaci MISE. Jestliže předpokládáme, že min CV(/i) ~ min E CV(h), dostaneme dobrou aproximaci optimální hodnoty h. Obrázek 3.9: Porovnání minima MISE (modrá, čárkovaná) a minima funkce křížového ověřování CV (červená, plná) pro simulovaná data z příkladu 3.1 OTÁZKA. Jak hodnotu bude mít odhad vyhlazovacího parametru pomocí metody křížového ověřování pro hustotu z ukázkového příkladu 3.2? Poznámka 3.6. Předpokládejme, že k — 2. Pak vychýlení odhadu může být velké, jestliže (/")2 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.10. Příklad 3.6. Jádrový odhad hustoty dat z příkladu 3.1 je zobrazen na obr. 3.11. Pro rekonstrukci bylo použito Epanečnikovo jádro a vyhlazovací parametr určený metodou křížového ověřování hCY = 0,5628. 42 0.5 - 0.4 - 0.3 - 0.2 - 0.1 - 0 ^1-3-2-101234 Obrázek 3.10: Zahlazení vrcholů a údolí při odhadu hustoty směsi normálních rozdělení, odhad (červená, plná), původní funkce (modrá, čárkovaná) 0. Obrázek 3.11: Odhad hustoty s hcv — 0,5628, odhad (červená, plná), původní funkce (modrá, čárkovaná) 5.4 Iterační metoda V tomto odstavci pojednáme o další metodě pro odhad vyhlazovacího parametru. Metodu pouze popíšeme a nebudeme se zabývat statistickými vlastnostmi získané aproximace. Podrobnou analýzu lze najít v monografii [3]. Tato metoda je založena na vztahu (3.7): AIV(/iopt,o,fc) = 2k AISB(/iopt,o,fc). Přepíšeme tuto rovnici 1 " nh """" (A;!)2 YVp. -2kh2k?^V(fM) = 0. (3.19) Problém nalézt hoptto,k, 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 vSrf(x, h) = \ í K2(y)f(x - hy) dy nh , bias/(x, h) = (Kh * f)(x, h) - f(x, h) = / f(x - hy, h)K(y) dy - f(x, h). Odtud plyne AW(h) = (3.20) nh 43 2 AISB(/i) = J f(x-hy,h)K(y)dy-f(x,h)j dx Výraz lze upravit pomocí konvolucí a dostaneme = -3- V (K*K*K*K-2K*K*K + K*K)^X,~X n2h ^ 1 ^-^ {Xi — X n2/i -^-^ \ /i Místo rovnice (3.19) řešíme rovnici V(K) 2k 2 dx. Funkce A(z) = (K * K * K * K — 2K * K * K + K * K)(z) má tyto vlastnosti zJA(z) dz = 0, j = 0,l,...,2Ä;-l, z2feA(z)dz = (2%l (3.21) AISB/(-, h) je vychýleným odhadem AISB, a proto budeme uvažovat _ i n AISB(/i) = Ah(-X* - (3'22) nh n Rovnici lze také zapsat ve tvaru: *-2frVnW7% x,=0. (3.23) 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.23) lze považovat za vhodnou aproximaci /iOpt,0,fe- Tato skutečnost je dokázána v následující větě [3]. Věta 3.5. Nectí P(h)^^§-2kh2^V(^) a V(K) 2k 44 Pak platí EP{h) = P(h) + o(h2k+1), vavP(h) = ^-V{K)V{f) + oin^h-1). nzh Poznámka 3.7. 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 hyř — hoptfl, hOpt,0,2 = O(„-i/i0) Rády rychlosti konvergence pro CV metodu a iterační metodu jsou stejné, ale výhodou iterační metody je podstatně menší výpočetní náročnost. Příklad 3.7. Jádrový odhad hustoty dat z příkladu 3.1 s Epanečnikovým jádrem a vyhlazovacím parametrem určeným iterační metodou je uveden na obr. 3.12. Obrázek 3.12: Odhad hustoty s hn — 0,5314, odhad (červená, plná), původní funkce (modrá, čárkovaná) 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/(-,/iopt,o,fc) =T(K) Ze vztahu pro /iopt,o,fe vypočteme V(f^) 2knh2ok+qk a tuto hodnotu dosadíme do předchozího vztahu: (2k + í)ôok AMISE/(-,/iopt,o,fc) =T(K) 2nkh, opt,0,k 45 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 I(ko) = [2 j, 3 = 0, ~ko~ _ 2 . Procedura 1. Pro k e I(k0) najděte optimální jádro Kopt0j„ e Sok, které je dáno tabulkou 8.1, k němu příslušný kanonický faktor Sok- 2. Pro k e I(ko) a Koptto,k <= Sok najděte optimální vyhlazovací parametr hoptto,k- 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. 4. Vypočtěte optimální hodnotu řádu k, které minimalizuje funkcionál L(k). 5. Použijte parametry z předchozích kroků ke konstrukci optimálního jádrového odhadu hustoty, tj. nflopt,0,k i=l ^ opt,0,fc' Příklad 3.8. Aplikace procedury na data z příkladu 3.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. k Kopt,0,k (>0fc hOpt,0,k L(K) 2 -i(-2- -1) 1,7188 0,5314 0,0141 4 ^■(x2 -32 vx l)(7x2-3) 2,0165 1,0734 0,0131 6 105 (2 256 [- - l)(33x4 - 30x2 + 5) 2,0834 1,6460 0,0125 8 315 (2 4096 v - l)(715x6 - -lOOlx4 + 385x2 - - 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.13. 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.14. 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[_iti](x), nebo jádro triweight: K(x) — ||(1 — x2)3I[_iti](x). Shrneme-li na závěr vypočítané hodnoty vyhlazovacích parametrů pro simulovaná data z příkladu 3.1, můžeme vizuálně porovnat jednotlivé odhady - viz obrázek 3.15. hopt,o,2 = 0,5122, ^ref = 0,9639, hMS = 1,0409, hCv = 0,5628, hIT = 0,5314, h = 1,6460 (^opt,o,6)- 46 Obrázek 3.13: Simulovaná data s jádrovým odhadem hustoty při použití procedury, odhad (červená, plná), původní funkce (modrá, čárkovaná) (a) Odhad s hopt,0,2 = 0,5122 (b) Odhad s hREF = 0,9639 (e) Odhad s hIT = 0,5314 (f) Odhad s h = 1,6460 a K e S0e Obrázek 3.15: Srovnání odhadů pro data z příkladu 3.1 48 (c) Odhad s hcv = 30,55 (d) Odhad s hiT = 28,18 (e) Odhad podle procedury s h = 140,46 a A" e Sos Obrázek 3.16: Grafy odhadnutých hustot pro obsah cholesterolu 7 Aplikace na reálná data Vraťme se k úvodnímu příkladu, který obsahuje měření množství cholesterolu v krevní plazmě u 371 pacientů. Data jsou shrnuta v tabulce 7.11. Skupina pacientů obsahovala dvě podskupiny, a to pacienty, u nichž bylo potvrzeno zúžení tepen (320), a pacienty zdravé (51). Máme-li k dispozici tuto informaci, mohli bychom očekávat, že odhadovaná hustoty může být bimodální. Avšak na druhou stranu, pokud při rekonstrukci hustoty neodhalíme tuto strukturu, neznamená to, že tam není, může být skrytá. S použitím všech uvedených metod pro odhad optimálního vyhlazovacího parametru jsme vypočítali tyto hodnoty: /iref = 29,28, hMS = 31,62, hCv = 30,55, hIT = 28,18. Výsledné odhady s Epanečnikovým jádrem jsou zobrazeny na obrázku 3.16, kde je také odhad hustoty při použití automatické procedury, u níž vypočítáme odhad optimálního vyhlazovacího parametru h — 140,46 a jádra Kos- 49 (a) Odhad s hREF = 5,7856 (b) Odhad s hMS = 6,2480 (e) Odhad podle procedury sft = 31,73aA"e Sos Obrázek 3.17: Grafy odhadnutých hustot pro délku krunýře Druhý 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.2 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 7.10. 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: /iref = 5,7856, hMS = 6,2480, hCv = 8,1317, hIT = 6,8263. 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 h — 31,7329 s optimálním jádrem Koptto,8- Výsledné odhady hustoty na jsou zachyceny na obrázku 3.17. 2Celý datový soubor je dostupný v programu R. 50 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/(,/*) ^YÍ^ + J^h2k(3l(K)V(f^). AIVW " äísb(I) ' Optimální vyhlazovací parametr vzhledem k AMISE pro odhad hustoty pravděpodobnosti s jádrem řádu k je tvaru > 2fc+l _ °0fc \k-i °Pt'°'fe 2knV(f(k))1 tj. Ziopt.o.fc = 0(n_1/(2fe+1)), AMISE /(■, /ioptjo,fe) = 0(n-2fe/(2fc+D). Metody pro odhad optimální hodnoty vyhlazovacího parametru h • metoda referenční hustoty f22kkrj^\^\ __i_ hREF = {-l2kjU-) Sok(7n 2k+1' • metoda maximálního vyhlazení • metoda křížového ověřování /2(x, h)áx--y f-i(Xi, h) 2 — 1 • iterační metoda /iit = resem rovnice: h - -—i-rp-ttt = 0. Automatická procedura pro simultánní volbu optimálního jádra, vyhlazovacího parametru a řádu jádra je dostupná v toolboxu Matlabu. Výstupy z výukové jednotky Student • zná tvar jádrových odhadů hustoty pravděpodobnosti • je schopen analyzovat statistické vlastnosti těchto odhadů 51 • se seznámil s metodami pro volbu vyhlazovacího parametru • porozuměl 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í Cvičeni 1. Dokažte vztah (3.4) pro tvar chyby AMISE. 2. Uvažujte náhodný výběr, který pochází z rozdělení s hustotou f(x) — 20x(í — x)3 pro x e [0,1] a který obsahuje 50 prvků. Vypočítejte hodnotu hoptQ^. Je výsledná hodnota správná? 3. Vypočítejte hodnotu optimálního vyhlazovacího parametru pro odhad s jádrem řádu 2 pro soubor dat s hustotou K(x) — y|(l — x2)2 pro x e [—1,1]. Pomůcka: V(f") = 22,5. 4. Dokažte: a) a\ — §\x x2gk(x) dx — 2k+5r Pro hustotu gk(x) Beta(/j + 2, k + 2) rozdělení (rovnice (3.16)). b) Jestliže hustota g má rozptyl a2, pak hustota ^-g(^-x) má rozptyl a2. 5. Dokažte vztah (3.15). 6. Spočítejte /ims Pr° • Epanečnikovo jádro K(x) — |(1 — x2), • kvartické jádro K(x) — j|(l — x2)2. 7. Aplikujte metody pro odhad vyhlazovacího parametru a automatickou proceduru na simulovaná data z ukázkového příkladu 3.2. 52 Kapitola 4 Jádrové odhady distribuční funkce 1 Motivace Distribuční funkce F (x) popisuje rozložení pravděpodobnosti náhodné veličiny X (budeme předpokládat spojitost náhodné veličiny). Stejně jako při rekonstrukci hustoty z množiny pozorovaných dat lze distribuční funkci odhadnout parametrickými nebo neparametrickými metodami. Zaměříme se výhradně na neparametrické metody, kdy předpokládáme pouze jistou hladkost odhadované distribuční funkce. Distribuční funkce má použití v analýze přežití, v teorii spolehlivosti, klasifikaci diagnostických testů (ROC křivky) aj. Funkce F(x) — 1 — F(x) se nazývá funkce přežití a vyjadřuje pravděpodobnost, že daná událost (úmrtí, selhání, porucha) nenastane dříve než po nějakém čase x. S touto funkcí se můžeme setkat i v hydrologii, kde se nazývá křivkou zabezpečenosti. Křivka zabezpečenosti je jedním z nejdůležitějších indikátorů toku řeky během daného časového období. Tato křivka vyjadřuje míru, s jakou je zajištěn vybraný průtok. Máme k dispozici průměrné měsíční průtoky řeky Dyje v profilu Vranov1 (pod přehradou Vranov, plocha povodí asi 2 220 km2) v období 1931-1990. Zajímá nás stabilita měsíčních průměrných průtoků. T O-m^^^^^B^W^^MHH-MlMIIIIIM lili +-HHH-++++ ++ + -H- 0 10 20 30 40 50 60 70 Obrázek 4.1: Průměrné měsíční průtoky řeky Dyje Máme-li tedy náhodnou veličinu X, která je popsána svou distribuční funkcí F(x) a k ní 1 Data byla poskytnuta CHMU, pobočka Brno. 53 příslušnou hustotou f (x), pak platí • 0 < F (x) < 1, • lim^-oo F (x) = 0, lim^^oo F (x) = 1, • F(x) je zprava spojitá, • F(x) — f*f(t)dt,Fŕ(x)=f(x). Nejužívanějším neparametrickým odhadem distribuční funkce je empirická distribuční funkce Fn. Ovšem Fn je schodovitá funkce i v případě, že F je spojitá. Nadaraya (1964) navrhl „hladkou" alternativu k f1, a to jádrový odhad F, který se získá integrací známého jádrového odhadu hustoty (3.1). 2 Základní typy neparametrických odhadů Nechf Xi,..., Xn jsou nezávislé náhodné proměnné, které mají tutéž spojitou hustotu / a distribuční funkci F. Nejjednodušší neparametrický dohad distribuční funkce F je empirická distribuční funkce Fn definovaná v bodě x vztahem 1 " F„(x) = - ^/(_00í2.](Xí). n i=l Tento odhad má sice dobré statistické vlastnosti, ale je to schodovitá funkce (viz obr. 4.2, a proto se budeme zabývat postupy, které umožní zkonstruovat „hladký" odhad distribuční funkce F. Příklad 4.1. Mějme dán náhodný výběr o velikosti n — 100 ze směsi dvou normálních hustot N(0; 4/9) a N (2; 1) s hustotou (viz kapitola 3 o hustotě) (Data jsou v tabulce 7.5.) 3 9a:2 1 0-2)2 f(x) = 0,5—7= e~V +0,5^ e—^- . Z obrázku 4.2 je patrné, že schodovitá funkce nevystihuje plně charakter distribuční funkce. Obrázek 4.2: Empirická distribuční funkce Fn (červená, plná) a skutečná distribuční funkce F (modrá, čárkovaná) pro data z příkladu 4.1 Nejznámější postup, jak odvodit neparametrický odhad distribuční funkce, spočívá v integraci jádrového odhadu hustoty, t.j. /X -t n í-X fit, h)dt=—Y,l K -co n" i=1 J — co t-Xi dt. 54 Užijeme-li substituce y — (t — Xi) jh, dostaneme -r—X- To znamená, že odhad F v bodě x e R je definován takto ' x — X, x - Xi F{x,h) = -Y,W W(x) = / K (t) dt. (4.1) Zde předpokládáme, že K e So2,K(x) > Oprox e [—1,1]. Níže jsou uvedeny základní vlastnosti funkce W: 1. W(x) — 0 pro x e (—oo, —1] a W(x) — 1 pro x e [1, oo), 2- /-i W2{x) dx < J\ W{x) dx = 1, 3. /_\ W{x)K{x) dx = i, 4. J\ xW{x)K{x) dx=\(l- J\ W2{x) dx). OTÁZKA. Lze použít obdélníkové jádro? Jaký bude tvar a vlastnosti funkce W(x)? Příklad 4.2. Použijeme-li Epanečnikovo jádro K(x) — | (1 — x2)I[_11] (x), pak funkce W je tvaru W(x) = < \(-x3 + 3x + 2) \x\<í x < -1, |x| < ] x > 1. Obrázek 4.3: Epanečnikovo jádro K (vlevo) a k němu příslušná funkce W (vpravo) Pro data z příkladu 4.1 je jádrový odhad distribuční funkce prezentován na obrázku 4.4. 3 Statistické vlastnosti odhadu Kvalitu jádrového odhadu lze lokálně popsat pomocí střední kvadratické chyby MSE: MSE F(x, h) = E(F(x, h) - F{x))2 = (EF(x, h) - F{x))2 + E(F(x, h))2 - (EF(x, h))2 . Spočítejme nejdříve hodnotu EF(x, h) v bodě ieR: EF(x,h)= íw(^-)f(y)dy cl f- OO h I W(t)f(x-ht)dt + h / W(t)f(x - ht)dt. > Ji 55 Obrázek 4.4: Jádrový odhad distribuční funkce s parametrem h — 1,5, odhad (červená, plná), původní funkce (modrá, čárkovaná) Předpokládejme dále, že F e C2. Označme první integrál li a druhý I2- Integrál ii počítáme metodou per partes a využijeme vlastnosti funkce W(t) h = h í W(t)f(x - ht) dt J — OO u = W(t) v! = W'(t) = K(ť) v1 = /(X - /rf)^ í; = - /ií) = [-W(t)F(x - ht)]\ + J F(x- ht)W'(t) dt = -F(x - h) + J K (t) F (x - ht) dt. (4.2) Dále použijeme Taylorův rozvoj F(x — ht) — F (x) — htF'(x) + ^-F"(x) + o(h2), tedy = -F(x -h) + F (x) + ^F"(x)h2p2(K) + o(h2). Počítejme nyní integrál I2: /OO W(t)f(x - ht) dt, uvažujeme-li substituce x — ht — z, dostaneme x: r-x—h f(z)dz= / f(z)dz = F(x-h). (4.3) x — h J — 00 Vychýlení odhadu je tedy tvaru Has F (x, h) = ^F"(x)h2p2(K) + o(h2). Poznámka 4.1. Vztahy (4.2) a (4.3) dávají zajímavý vztah pro vychýlení EF(x, h) - F(x) = -F(x - h) + J K (t) F (x - ht) dt + F(x -h)- F (x). Odtud plyne EF(x,h) = J K (t) F (x - ht)dt. 56 A dále (z Taylorova vzorce) EF(x, h)= K(ť) F(x) - htF'(x) + -h2t2F"{x) + o{h2) dt F{x) + -h2p2{K)F"{x)+o{h2) =o{h) = F{x) + o{h). Nyní dokážeme tvar rozptylu. 1 var F(x,h) = - ( EW . n \ \ h x-X -EZW x-X Zde E2W (3^2£) = {EW (3^2£))2 = (F(x) + o(h))2. Počítáme tedy jen integrál (označme jej J3): h = 1 |substituce: x — y — th\ 1 n VK2(í)/(x-/ií)dí + /i / f{x-hť)út =F(a;-7i) První integrál počítáme metodou per partes a máme 1 2 /* 1 J3 = - [-F(x - ht)W2(t)} * + - / F(x - /ií)VK(í) W(í) dí + -F(x - h) 71 L J ^ Ti I 71 = --F(x - /i) + - / F(x - /iíW(í)lf (í) dí + -Fix - h), n n n použijeme nyní Taylorův rozvoj funkce F = - / W(t)K(t) (F(x) - /iíF'(x) + o{h)) dí n = -F(x)\ W{t)K{t)dt--hF'(x) tW(t)K(t)dt +o^k = i(vlasnost3) =i (l —/11 (a;) da;) (vlastnost 4) užitím vlastností funkce W a F'(x) — f(x) dostaneme F(x) - hf(x) 1 - / W2(t)dt h o | - n Rozptyl je tedy tvaru var Fix, h) — — n F (x) - h f (x) ( 1 - / W2(t)dt i + o[-)--(F(x)+o(h))2 n J n = -F(x)(l - F(x)) --f(x) (1 - f W2(t) dt) + o (-Výše uvedené výsledky můžeme nyní zformulovat v následující větě: 57 Věta 4.1. Nectí F e C2, h —>• O, n/i —>• oo pro n —>• oo. PflA: MSEF(i,/í) = -F(i)(l - - -/í/(i) I 1 - / W2(t)dt ^ (4-4) + 1-(F"(x))2h4(322(K) + o(-+h4" 4 \n Globální pohled na kvalitu odhadu lze získat prostřednictvím střední integrální kvadratické chyby (MISE). Věta 4.2. Nectí F e C2, V {F") = / (F" (x))2 dx < oo, K e S02, lim„^oo /i = 0 fl lim^oo n/i = oo. MISEF(-,/i) = - í F(x)(í-F(x))dx-Cl- + c2h4+ o( - + h4) , (4.5) n J n \n J kde ci = 1 - / T^2(í) dí, c2 = ±fá(K)V(F"). Naším cílem je nalézt takovou hodnotu vyhlazovacího parametru, pro kterou bude MISE nabývat minimální hodnoty. Ale uvedený tvar MISE není pro takovou analýzu vhodný, a proto (stejně jako při odhadu hustoty a regresní funkce) budeme uvažovat asymptotickou střední integrální kvadratickou chybu AMISE, která v tomto případě je tvaru: AMISE F(-,/i) = - / F(x)(l-F(x))dx-Cl-+ c2h4 . (4.6) n J n AIV(h) AISB(Ti) Nyní už lze standardními metodami matematické analýzy nalézt takovou hodnotu h, pro kterou AMISE nabývá minimální hodnoty. Je snadné ukázat, že /iopt.o.2 = n"1/3 (^) = Oln-W) (4 7) a pak AMISEF(; hoptfi,2) = -/ F(x)(í-F(x))dx- — (^) n-4'\ (4.8) n J c2 Poznámka 4.2. Optimální hodnota vyhlazovacího parametru pro odhad distribuční funkce je řádu 0(n_1/3), zatímco pro odhad hustoty s jádrem K e So 2 je vyhlazovací parametr řádu 0(n_1/5). Ukázkový příklad 4.3. Předpokládejme, že známe tvar distribuční funkce F(x) — 5x7 + 2íx5 — 35x3 + 35x + 16) pro x e [—1,1]. Vypočítejme hodnotu optimálního vyhlazovacího parametru pro odhad s jádrem řádu 2. Podle vztahu (4.7) potřebujeme spočítat hodnoty c\ a c2. S Epanečnikovým jádrem je ci = 1 / W2(x)dx = 1 - / — (-x3 + 3x + 2)2da; = 0,2571, 7-i 7-1 16 C2 = i/32(^)y(0 = TO .3,1818. Pak /ioptj0j2 = 1,2642-n-Vs. Na obrázku 4.5 je odhad s optimálním vyhlazovacím parametrem pro náhodný výběr o 50 pozorování, která pochází z rozdělení s uvedenou distribuční funkcí (data jsou v tabulce 7.4). 58 0.2 Obrázek 4.5: Odhad distribuční funkce z ukázkového příkladu 4.3, odhad (červená, plná) a původní funkce (modrá, čárkovaná) za použití Epanečnikova jádra a /iopt,o,2 — 0,3432 I v tomto případě je volba jádra méně důležitá než volba vyhlazovacího parametru. Lze doporučit jádra třídy So 2/ např. • Epanečnikovo K(x) — |(1 — x2)I[_11](x), • kvartické K(x) — y|(l - x2)2I[_1 aj(x), • triweight K(x) = ||(1 - x2)3/[_ia] (x). 5 Volba vyhlazovacího parametru 5.1 Metody křížového ověřování Metody křížového ověřování patří k nejužívanějším metodám pro volbu vyhlazovacího parametru. Zde uvedeme pouze metodu navrženou A. Bowmanem (1998). Funkce křížového ověřování je v tomto případě tvaru kde F-i(x, h) je jádrový odhad distribuční funkce s vynecháním bodu A4. Pak /icv — arg min CV(/i), h£Hn přičemž Hn — [an-1/3, brT1^] pro vhodná 0 < a < b < 00. 5.2 Princip maximálního vyhlazení Myšlenka této metody je stejná jako pro odhad hustoty. Užijeme-li faktu, že 4 Volba jádra můžeme aplikovat Terrelovu větu 3.3 pro k — 1. V tomto případě je 59 a tedy kde (Ti — J x2g1(x) dx — j, V(g[) — y. Odtud plyne, že ti je odhadem a (viz rovnice (3.13) a (3.14)). Poznámka 4.3. Hodnota h-^s může sloužit jako horní hranice pro množinu vyhlazovacích parametrů volených podle metody křížového ověřování. Tedy Hn — [h^, hyis], kde ha je nejmenší vzdálenost mezi po sobě jdoucími body Xt■, í — 1,..., n. Příklad 4.4. Pro data z příkladu 4.1 zvolme Epanečnikovo jádro. Pak hodnoty potřebné pro odhad vyhlazovacího parametru metodou maximálního vyhlazení jsou následující: i ŕ n =100, £ = 1,3426, l32(K) = -, Cl = l - J W2(x) dx — 0,2571. Pak platí /ims — 1,1037 a na obrázku 4.6 je zobrazen odhad distribuční funkce. Obrázek 4.6: Odhad distribuční funkce s /ims — 1,1037, odhad (červená, plná), původní funkce (modrá, čárkovaná) 5.3 Plug-in metoda Společným cílem metod typu plug-in (PI) je odhadnout V(F"). Za předpokladu dostatečné hladkosti funkce / užitím metody per partes dostaneme vztah V{F") = J {F"{x)fdx = -j f'(x)f(x)dx. Tudíž se budeme dále zabývat odhadem funkcionálu v>i=y f(x)f(x)dx. Je zřejmé, že ipi — Ef"(X), což vede k metodě založené na odhadu druhé derivace hustoty /. Vztah (3.9) použijeme k odhadu druhé derivace s jádrem — Kopt,2,4 € S24. A = n-1 J2 f"(Xi, h) = n-2h-z J2 E ^(2) Xi — Xj 60 kde podle vztahu (3.11) je #04 hQpt,2,4 — lO1''9"?—hopt, 0,4- Pak c2 = -i/3|(^)íi. Shrnutím předchozích úvah dostaneme proceduru pro odhad distribuční funkce F: Procedura 1. Najděte optimální vyhlazovací parametr /iOpt,0,4 pro odhad hustoty s optimálním jádrem KoptflA e S04. 2. Najděte optimální vyhlazovací parametr /i0pt,2,4 pro odhad druhé derivace hustoty podle vztahu (3.11) s k — 4 a optimálním jádrem e S24. 3. Vypočtěte odhad funkcionálu i/>i s využitím hodnoty /i0pt,2,4 získané v kroku 2. 4. Vyčíslete optimální hodnotu vyhlazovacího parametru hPI = n-^ {-^-) ' 5. Použijte parametry z předchozích kroků ke konstrukci optimálního jádrového odhadu distribuční funkce F(x, h) s daným jádrem K e So2- Příklad 4.5. S použitím funkce toolboxu zjistíme, že pro data z příkladu 4.1 je vyhlazovací parametr určený plug-in metodou roven hpi — 0,5717. Na obrázku 4.7 je odhad distribuční funkce společně se skutečnou distribuční funkcí. Obrázek 4.7: Odhad distribuční funkce s /ipi — 0,5717, odhad (červená, plná), původní funkce (modrá, čárkovaná) 6 Aplikace na reálná data Znovu se podívejme na soubor dat z povodí řeky Dyje (viz tabulky 7.13 a 7.14). Podle metod pro odhad vyhlazovacích parametrů jsme pro Epanečnikovo jádro vypočítali tyto odhady hMS = 1,9808, hpi = 0,3636. 61 (c) Empirická distribuční funkce Obrázek 4.8: Odhadnuté křivky zabezpečenosti F (x) řeky Dyje v období 1931-1990, na ose x je velikost průtoku v m3/s Na obrázku 4.8 jsou zobrazeny odhady křivky zabezpečenosti (F (x) — 1 — F (x)) společně s empirickou křivkou zabezpečenosti. Průběh křivky zabezpečenosti ukazuje například, že průtok 8m3/s a vyšší se vyskytuje s pravděpodobností asi 0,3, průtok 20m3/s a vyšší se vyskytuje s pravděpodobností 0,1. Provoz přehrady tedy stabilizuje měsíční průtok v rozmezí 5 — 8m3/s s pravděpodobností 0,9 — 0,3. Následující datový soubor pochází z rozsáhlé studie, v níž autoři studovali vliv substituentů v 2,4-diamino-5-(substituovaný benzyl)pyrimidinech. Biologická aktivita při inhibici dihydro-folát reduktázy byla měřena pomocí asociační konstanty. Data jsou v tabulce 7.12 a jsou dostupná na osobních stránkách Dennise D. Boose2, kde je také odkaz na původní článek Jonathana D. Hirsta z roku 1994. Užitím výše uvedených metod jsme (při použití Epanečnikova jádra) dostali následující hodnoty vyhlazovacích parametrů: hMs = 0,1139, hpi = 0,1931. Na obrázku 4.9 jsou uvedeny odhady distribuční funkce s těmito parametry a také je zde pro srovnání uvedena empirická distribuční funkce. 2 http://www4.stat.ncsu.edu/~boos/var.select/pyrimidine.html 62 (a) Odhad s /ims = 0,1139 (b) Odhad s hpi = 0,1931 (c) Empirická distribuční funkce Obrázek 4.9: Odhadnuté distribuční funkce pro vliv substituentů v pyrimidinech 63 Shrnutí Odhad distribuční funkce F(x) v bodě x je tvaru n*) = \ E w (^ŕ)' w{x) = /lm áf- Asymptotická střední kvadratická chyba jádrového odhadu distribuční funkce je součtem asymptotického tvaru rozptylu (AIV) a druhé mocniny vychýlení (AISB) 1 ľ h AMISEF(-,/i) = - / F(x)(í- F(x))dx-Cl-+ c2K n J n AlV(h) kde ci = 1 - / W2(t) dt, c2 = -B2{K)V{F"). 4 ' AISB(Ti) 4' Optimální vyhlazovací parametr vzhledem k AMISE pro odhad distribuční funkce je tvaru hopt,o,2 = n-v3 (Jj-J , t.j. /iopt.o.2 = C^n-1/3). Metody pro odhad optimální hodnoty vyhlazovacího parametru h • metoda křížového ověřování hcv — argmin^g^ CV(/i) CVW = ^E / (l(-oo,x](xi) - F-i(x>h))2 i—l • metoda maximálního vyhlazení ^=""/3(iáwí)"s^ • plug-in metoda Výstupy z výukové jednotky Student • zná základní typy jádrových odhadů distribuční funkce a jejich statistické vlastnosti • získal přehled o metodách pro volbu vyhlazovacího parametru • je schopen navrhnout a implementovat proceduru pro zpracování reálných dat • se naučil používat příslušný toolbox v Matlabu a dokáže zkonstruovat jádrový odhad distribuční funkce pro daná reálná data 64 Cvičeni 1. Odvoďte tvar funkce W(x) pro kvartické jádro K(x) — y|(l — x2)2I[_11](x). 2. Dokažte vlastnosti 2,3 a 4 funkce W. 3. Dokažte vztahy (4.7) a (4.8). 4. Odvoďte tvar vyhlazovacího parametru podle metody maximálního vyhlazení pro kvartické jádro. 5. Aplikujte metodu maximálního vyhlazení a plug-in metodu na simulovaná data z ukázkového příkladu 4.3. 65 Kapitola 5 Jádrové odhady dvourozměrných hustot 1 Motivace Uvažujme soubor dat ze studie The State of the Worlďs Children (UNICEF), který obsahuje míru úmrtnosti dětí od narození do pěti let věku počítanou na 1000 živě narozených dětí a očekávanou délku života narozeného dítěte (s ohledem na úmrtnost v dané populaci v době jeho narození) v 72 zemích, které měly v roce 2001 hrubý národní produkt menší než 1000 amerických dolarů na osobu a rok. Chceme vědět, jaká je pravděpodobnost, že dítě narozené v zemi s vysokou Obrázek 5.1: Studie UNICEF - míra úmrtnosti (osa x) a očekávaná délka života u dětí (osa y) mírou úmrtnosti, se dožije vyššího věku. Případně, zda se v datech nevyskytují shluky, tj. zda jde o vícemodální hustotu. V tomto případě se jedná o odhad dvourozměrné hustoty. Odhady vícerozměrných hustot se budeme zabývat v této kapitole. 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) z kapitoly 3, a zaměříme se zejména na odhad dvourozměrné hustoty. Tedy naším cílem je rekonstruovat hustotu pravděpodobnosti z náhodného výběru. Poznámka 5.1. Jádrové odhady dvourozměrných hustot se obvykle znázorňují pomocí vrstevnic, které umožňují snazší náhled na odhadnutou funkci. 66 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 5.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) (Data jsou v tabulce 7.6.) f(x v) = lJLe-l(*2Hy+tr)+±±e-h(*2+(y-2f)+L±e-l(^Hy-4)2) Jy 3 2tt 3 2tt 3 2tt Z obrázků 5.2(a) a 5.2(b) je patrné, že histogram nepostihuje charakteristické rysy hustoty pravděpodobnosti dat, která je zobrazena na obrázku 5.2(c). 1 1 X2\ 1 3 '-"X 1~ (a) 7 x 7 obdélníků (b) 5 x 10 obdélníků (c) Hustota simulovaných dat Obrázek 5.2: Histogramy s různými počty třídicích obdélníků a hustota pro data z příkladu 5.1 Přejdeme nyní k jádrovým odhadům dvourozměrné hustoty. Předpokládejme, že máme k dispozici náhodný výběr ([Xi, Y\],..., [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 -t n i n f(x, y, H) = - Y, Kn(x - Xity - Y) = YK{H^(X ~ X^ V ~ Y^ (5.1) K je dvourozměrné jádro a H je pozitivně definitní matice typu 2 x 2 a |H| značí její determinant. Prvky matice H se nazývají vyhlazovací parametry a matice se pak zkráceně označuje jako vyhlazovací matice. 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), Používáme zde zkrácený zápis pro dvourozměrnou hustotu normálního rozdělení, a to N(fii, [12', a\, a\, Ě>) 67 • sféricky symetrické jádro Ks{x, y) — ck XK\ (\Ac2 + y2), ck — j j Ki (\Ac2 + y2) dx dy. Příklad 5.2. Epanečnikovo jádro, které je v jednorozměrném případě tvaru K(x) — |(1 — x2), 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 sféricky symetrického Epanečnikova jádra v závislosti na třídě matic. 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, Yj] 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/(x, y, H) = ^-((K^ * f)(x, y) - (Ku * f)2(x, y)) +((KU * f)(x y) - f(x,y))2, 68 Tabulka 5.1: Třídy vyhlazovacích matic S V h\ 0 0 h2 h\ h12 h12 h\ 2 -2 2 -2 Obrázek 5.4: Konstrukce jádrového odhadu hustoty nebo globálně pomocí střední integrální kvadratické chyby MISE/(-,H) = // MSE/(x,y,H)dxdy. 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ýpočty pro matici třídy T jsou velmi náročné, a proto se v dalších úvahách omezíme na odhady s diagonální maticí. Věta 5.1. Předpokládejme, že funkce f, jádro K a diagonální matice vyhlazovacích parametrů H — 'h\ 0 diagOí,^) 0 h2 splňují následující předpoklady: (i) Nectí H — H„ je posloupnost matic takových, že {nh\h\) 1 a prvky h\ a h\ konvergují k nule pro n —> oo. (íí) Dále nectí všechny druhé parciální derivace funkce / jsou ohraničené, spojité a integrovatelné se čtvercem. (ííí) Jádro K splňuje xK{x,y) áx áy — J J y K(x, y) dx dy — 0, x2 K (x, y) dxdy — // y2 K (x, y) dxdy — /32(K). Pak platí MISE/(•, H) = AMISE /(■, H) + o(hj + h\) + o((/i1/i2n)-1), 69 kde AMISE/(, H) = + -^2{K)(h\V{fxx) + 2h2h2V(fxy) + h42V(fyy)), (5.2) AIV(H) AISB(H) přičemž označení je ve shodě s předchozími kapitolami, tj. V(g) — JJ g2(x, y) dx dy. 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\, h2, pro které AMISE nabývá minimální hodnoty jsou dány vztahy: ' W\fyy)V{K) X 1/6 Kopt \nft{K)V*/\fxx) [V(fxy) + V^(fxx)V^(fyy)} 1 ' (53) V3'\fxx)V{K) X 1/6 h2'°pt \n^{K)Vy\.fyy) [V{.fxy) + V^(fxx)V^(fyy)}) ■ {5A) Z těchto vztahů plyne, že množina přípustných vyhlazovacích parametrů je tvaru ain-1/6 < hi < birT1^, a2rnT1^ < h2 < b2nrxl^ pro vhodné konstanty 0 < a\ < b\ < oo, 0 < a2 < b2 < oo. Ukázkový příklad 5.3. Předpokládejme, že známe tvar hustoty f(x, y) — ^ e-3^' +v ^2 pro [x, y] e M2. Vypočítejme h Epanečnikovým jádrem. 2vr [x, y] eM2. Vypočítejme hodnoty optimálních vyhlazovacích parametrů pro odhad se součinovým Podle vztahů (5.3) a (5.4) potřebujeme spočítat výrazy 2 V(fxx) = jj fxx(x, y)dxdy = JJ fxxxx(x, y)f(x, y) dxdy = 34 • (2tt) \ V(fxy) = U J2xy{x,y)dxdy = II fxxyy(x,y)f(x,y)^33 ■(2tt)-\ v(fyy) = II fyy(x,y) dxdy= j j fyyyy(x,y)f(x,y) dx dy = 34 • (27i-)-1. Pro součinové Epanečnikovo jádro K(x, y) — y|(l — x2){\ — y2) platí V{K) = II K2(x,y) dxdy = 0,36, P2{K)= íí x2K(x,y)dxdy= U y2 K (x, y) dx dy = 0,2. Pak pro optimální vyhlazovací parametry máme (34 • (27t)-1)3/4 • 0,36 ILl,opt ^__(34 ■ (27t)-1)3/4 ■ 0,36 0,22(34 • (2tt)-1)3/4[33 • (27t)-1 + (34 • (2tt)-1)1/2(34 . (27r)-i)i/2] L2,opt 0,22(34 • (2tt)-1)3/4[33 • (27t)-1 + (34 • (2tt)-1)1/2(34 . (2tt)-1)1/2 n Tedy hx = 0,8978 • n"1/6, h2 = 0,8978 • n"1/6. Oba parametry jsou totožné, což se dalo očekávat, protože hustota f(x,y) je symetrická jak podle osy x, tak podle osy y. Tedy pro soubor o 100 prvcích bude vyhlazovací matice rovna Hopto,2 — diag(/i2,/i2,) — diag(0,776, 0,776). Odhad s optimální vyhlazovací maticí pro tento datový soubor (viz tabulku 7.7) jsou na obrázku 5.5. 2Použili jsme metodu „per partes" pro vícerozměrné integrály. 70 Obrázek 5.5: Odhad hustoty z ukázkového příkladu 5.3, odhad (červená, plná) a původní funkce (modrá, čárkovaná) při použití součinového Epanečnikova jádra 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 5.3. V literatuře se také využívá Gaussovo jádro K(x, y) — (2tt)^1 e~(x +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 1 _ *2 _ y2 f(?,V) = ô-e ^ ^ 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 tento odhad vyhlazovací matice H RE F Tfn-1/3 0 0 a2n-V3 kde a2 jsou odhady rozptylů pro jednotlivé proměnné. Obecnější tvar tohoto vztahu, který je znám také jako Scottovo pravidlo, lze nalézt v [13]. 71 Příklad 5.4. Pro simulovaná data z příkladu 5.1 je matice vyhlazovacích parametrů podle metody referenční hustoty rovna s využitím Gaussova jádra takto: H REF 0,1292 0 0 0,9945, Na obrázku 5.6 je vykreslen odhad hustoty s pomocí Gaussova jádra a porovnání odhadu se skutečnou hustotou. V X : * > (a) Data -2-1 0 1 2 3 (b) Odhad s href -1 o 1 (c) Skutečná hustota Obrázek 5.6: Jádrový odhad dvourozměrné hustoty z příkladu 5.1 s Gaussovým jádrem - metoda referenční hustoty 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, Yj] s vynecháním tohoto pozorování. Funkci metody křížového ověřování CV můžeme zapsat ve tvaru CV(H)= /7(/(x,y,H))2dxdy-^/_í(Xí,yí,H). i—l kde 1 " f-i(Xi, li, H) = —j K^ - - yj) n 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.2. Funkce CV je nevychýleným odhadem MISE, tj. platí E CV(H) = MISE /(-,H) - J J f (x, y) dxdy. Důkaz. Vypočtěme střední hodnotu CV: ECV(n) = E í í f2{x1y1K)dxdy--YJEf_l{XllYllll) •* •* n i=l = E U f2(x, y, H) dxdy - 2EKU(X1 - X2,Yľ - Y2) = E ľ (x, V, H) dxdy - 2 //// Ku(x-u,y-v)f(x,y)f(u,v)dxdydudv 72 a úpravou MISE MISE/(.,H) = £ JJ(f(x,y,H)-f(x,y))2dxdy = eJJ f2(x,y,H)dxdy-2E J f(x,y,H)f(x,y)dxdy + J J f2(x,y) dxdy = E JI f2(x,y,n)dxdy - 2 jjjj KH(x - u, y - v) f (x, y) f (u, v) dxdy du dv f (x, y) dxdy. 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). Příklad 5.5. Použijeme-li součinové Epanečnikovo jádro, pak pro simulovaná data z příkladu 5.1 dostanem matici vyhlazovacích parametrů určenou podle metody křížového ověřování v následujícím tvaru: 11,4532 0 H-cv — \ 0 0,1431, Na obrázku 5.7 je vykreslen odhad hustoty s touto maticí. A* -3-2-10123 (a) Data (b) Odhad s hc (c) Skutečná hustota Obrázek 5.7: Jádrový odhad dvourozměrné hustoty z příkladu 5.1 se součinovým Epanečnikovým jádrem - metoda křížového ověřování 6 Aplikace na reálná data Vraťme se k datovému souboru z motivačního odstavce. Jak už jsme zmínili, data pochází ze studie Tlie State of the Worlďs Children (UNICEF) a obsahuje míru úmrtnosti dětí od narození do pěti let věku počítanou na 1000 živě narozených dětí a očekávanou délku života narozeného dítěte (s ohledem na úmrtnost v dané populaci v době jeho narození) v 72 zemích, které měly v roce 2001 hrubý národní produkt menší než 1000 amerických dolarů na osobu a rok.3. Data jsou shrnuta v tabulce 7.15. 3http://www.unicef.org/sowc0 3/tables/tablel.html 73 Obrázek 5.8: Vrstevnicové grafy odhadnutých hustot pro data ze studie UNICEF - na ose x je míra úmrtnosti dětí do pěti let (na 1000 živě narozených dětí) a na ose y očekávaná délka života narozeného dítěte Odhady vyhlazovacích matic podle uvedených metod jsou tyto /ll45,2 0 \ A.74,14 0 \ íIref — ticv — \ 0 24,9y y 0 59,35y Na obrázku 5.8 jsou znázorněna data, vrstevnice jádrového odhadu s Gaussovým jádrem a maticí určenou metodou referenční hustoty a vrstevnice odhadu se součinovým Epanečnikovým jádrem a maticí určenou podle metody křížového oěřování. jádrem. Následující 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 a u nichž bylo potvrzeno zúžení tepen. Data jsou shrnuta v tabulkách 7.16 a 7.17. Matice vyhlazovacích parametrů určené podle metody referenční hustoty a metody křížového ověřování (s Epanečnikovým jádrem) jsou následující: 1270,5 0 \ /l797,2 0 \ íIref — tlcv — \ 0 1516,5 y y 0 892,7y Na obrázku 5.9(a) jsou znázorněna data, na obrázku 5.9(b) jsou vrstevnice jádrového odhadu s Gaussovým jádrem a maticí určenou metodou referenční hustoty a na obrázku 5.9(c) jsou zobrazeny vrstevnice odhadu se součinovým Epanečnikovým jádrem a maticí určenou podle metody křížového oěřování. jádrem. 74 Obrázek 5.9: 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) 75 Shrnutí Odhad dvourozměrné hustoty pravděpodobnosti f(x, y) v bodě [x, y] je tvaru 1 n /(*, y>n) = ^\ J2K(n^^x -x»y- Yi)q n' ' i=i Dva typy jader: • součinové jádro: Kp(x,y) — Ki(x) ■ K\(y), sféricky symetrické jádro: Ks{x, y) — ck xK\{^Jx2 + y2), Ck — jj K± (\px?~-~Vy2) dx dy. Asymptotická střední integrální kvadratická chyba dvourozměrného jádrového odhadu AMISE/(-, H) = + \p2{K)(h\V{fxx) + 2h2h2V(fxy) + h2V(fyy)). Optimální vyhlazovací parametry vzhledem k AMISE \ 1/6 _V3/\fyy)V(K) '°pt [n(32(K)V^(fxx) [V(fxy) + y/V(fxx)V(fyy)\ hl'°pt \np^K )V*'*VXX) [V [fxy) + V V [fxx)V {.fyy)\ J 1/6 _ /_V3/\fxx)V(K)_ l2'°pt " \nP2(K)V^(fyy)[V(fxy) + ^V(fxx)V(fyy)\ Metody pro odhad optimálních hodnot matice vyhlazovacích parametrů H — diag(/ii, h^) • metoda referenční hustoty /ij.REF = ^n-1/6, 2 = 1,2, • 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 Výstupy z výukové jednotky Student • zná součinová a sférická dvourozměrná jádra pro odhady dvourozměrných hustot • porozuměl principu vyhlazování dvourozměrných hustot • pochopil 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 76 Dodatek Symetrická matice A je pozitivně definitní právě tehdy, když všechny její hlavní minory (subde-terminanty) jsou kladné, tj. platí A f «11 «12 «21 «22 \«nl «n2 «ln «2n Mu «ii > 0, M22 «11 «12 «21 «22 > 0,..., Mnn = \A\ > 0. Cvičeni 1. Určete součinové a sféricky symetrické dvourozměrné jádro odvozené z kvartického jádra ^(x) = lf(l-x2)2. 2. Odvoďte vztahy (5.3) a (5.4) pro optimální vyhlazovací parametry. 3. Ukažte, že pro optimální vyhlazovací matici vzhledem k AMISE /(■, H) platí AIV(Hamise) = 2 • AISB(Hamise), kde Hamise — arg niinHex> AMISE. 4. Aplikujte metodu referenční hustoty a metodu křížového ověřování na simulovaná data z ukázkového příkladu 5.3. 5. Ověřte, že matice A je pozitivně definitní maticí A = 6-12 -1 4 -3 2-3 9 77 Kapitola 6 Návody ke cvičením a odpovědi na otázky Kapitola 1 Cv. 1 Je zřejmé, že jádro K1 má nosič [—7,7]. Do funkcionálu T dosadíme jádro K1(ť): T(KS) k —h K2s{ť) dt -8 tk K s (t) dt -8 k — i 2v+l tk-^K (f- I dt -8 ôv+1 \ ô 2v+l _i 52v+2 ô K2 (x) dx k—v 1 £fc 2^+1 x ——SK(x) dx íí2(aľ) dx k —u xkK(x) dx 2v+l 1 5(2^+l)(fe-1.) ^(fc-^^+l) = T(K) Cv. 4 Výraz (X~ * K)(x)x3 dx — J J K (z) K (x — z)x3 dz dx spočítáme postupně pro jednotlivé hodnoty j. j = 0 : J J K(z)K(x - z) dz dx = y (z) y K(x - z) dx dz — I substituce vnitřního integrálu: x — z — w\ — J K (z) J K(w) dw dz = J K{z) dz ■ J K(w) dw = 1 Podobně pro j — 1 a j — 2, přičemž využíváme vlastností jádra (definice 1.1). 78 Cv. 5 Ověříme podmínky jádra třídy S02, tj. • K(x) dx — 1 —>• odtud spočítáme hodnotu konstanty A — ^, • xK{x) dx — 0 —>• metodou per partes spočítáme nebo úvahou o integrování liché funkce přes symetrický interval odvodíme, že vztah platí, • j\ x2K{x) dx — 0 —>• dvojím použitím metody per partes vypočítáme hodnotu /?2 (Ä") = 1 — ^ = 0,1894, která je nenulová, tudíž jde o jádro třídy S'o2- Nakonec, s využitím vztahu pro druhou mocninu funkce cos x, tj. cos2 x — (1 + cos 2x)j2, spočítáme hodnotu V(K) = j\ K2(x) dx = = 0,6169. Kapitola 2 Otázka na str. 15 Obdélníkové jádro: Kh(z) — ^ pro z G [—/i, /i]. Váhy Nadarayova-Watsonova odhadu: EJ=i Kh(x0 - xj) YTj=i éi^-hM^o - xj) _ "éhh-hM (x0 ~ Xi) _ h-h,h] (x0 ~ xi) Jkh-hM (xo - Xj) nl[_h^ (x0 - Xj) Tedy NW odhad s obdélníkovým jádrem je totožný s klouzavým průměrem (viz obr. 2.7(b)). Cv. 1 Obdélníkové jádro: Kh(z) — ^ pro z e [—h, h]. Priestleyův-Chaův odhad: ^ n 1^1 1 n fhpc(x,h) = - ^Kh(x - xl)Yl = - ^2 2fol[-h,h\(x ~ xí)yí = ^Jí ^ YlI^h^(x ~ x^m i—l i—l i—l Princip je podobný jako u Nadaraya-Watsonových odhadů (str. 15). Cv.2 yn 1K (x-Xl\y V™ -J-K (X-X*)Y fhNW(x, h, K) = U=lh,[/h JT =\h = h*S\ = ™ V 2^1=1 h1^ V h ) 2^1=1 h*S V h*8 ) ET=1 (V) * Eti ^ (2^) Yi 2-,i=i h* áxM <5 ; Cv. 3 Do vztahu (2.13) dosadíme hodnoty V(K) = 1,250 a /34(FT) = -0,0476 (z tabulky 1.2) a z regresního modelu cr2y(iT)(fc!)2 opt'°'k ~ 2knP2(K)V(m(k)) h2kf}^ 0,25 • 1,25 • (4!)2 A y ' - 3,2705 • 10~4, optflA 2 • 4 • 100 • (-0,0476)2 • 32tt8 hopt,o,4 — 0,41. Cv. 5 Ve vztahu (2.7) použijeme Taylorův rozvoj i(x - hu) = m(x) - hum'ix) H-----h , , hkuk(x) + o(hk). 79 Využijeme vlastnosti jádra K e Sofc a dostaneme hk Em(x, h) = mix) + i-í)k—l3k(K)m^k\x) + o(hk) + Oŕn"1). Dále postupujeme jako při odvození vztahu (2.11). Cv. 6 Spočítejme nejdříve derivaci AMISE to(-, h) (viz vztah (2.12)) dAMIS!f(-'fe) = + >»(m«). d/i n/i2 (k!)2 Jelikož hledáme minimum funkce AMISE íh(-,h), položíme tuto derivaci rovnu nule Nechť hoptfl,k je řešením této rovnice. Dále vynásobíme tuto rovnici parametrem /iOpt,0,fe ct2V(K) 2k neboli nhoptfl,k (kl) ^K1^2k- _ ■ 2fc \\2 opt,0 ^2{K)V{m^) = 0, 1 j^h2kt^kPl{K)V{m^) nhopt,0,k AlV(hoptfi,k) = 2fcAISB(/iopt,0,fc) Kapitola 3 Otázka na str. 31 Obdélníkové jádro: Kh(z) — ^ pro z e [—/i, /i]. Odhad hustoty _^ n 1^1 1 n f{x,h) = -^Kh{x-Xi) = -^—I[_hM{x-Xi) ^^Hh-hMÍx-Xi)- i=l Tedy v každém bodě Xi sestrojíme obdélník se šířkou 2h a výškou (2nh)~1, tyto odhady se pak sečtou. Tomuto typu odhadu se také říká naivní odhad. Naivní odhad v jistém smyslu osvobozuje histogram od volby polohy třídicích intervalů. Tento odhad „klouže" přes data - jako histogram pro odhad hustoty odpovídá regreso-gramu pro odhad regresní funkce, tak výsledek odhadu hustoty s obdélníkovým jádem odpovídá klouzavému průměru u regrese. 1/(2h)- Cv.l AMISE/(, h) = + -±^h2kp2{K)V{.f^) = T(K) V(K) 1 V(K) nhT(K) (k\)2 Pomocné výpočty: T(K) = S2k(32(K), S2^+1 = V{K) _ V{K) Ť(K) ~ S2k(32(K) (32(K)^ (32(K) ^ 1 T(K) 52kMK) 52k k2krk(K) T(K) C+1 _ g Ok, 80 Celkem AMISE/(-./i) =T(K) ( Sok nhT(K) (k\)2 1 h2kV(f(kK x2k °0k C v. 2 Pro danou hustotu f(x) — 20x(l — x)3 spočítáme její čtvrtou derivaci: (x) — —480, pak V(f^) — 4802. Dosazením do vztahu (3.5) s využitím vlastností jádra K e So4 dostaneme €+1(^!)2 h2k+l _ 9 = ^094(4!)2 = réoTw^2 opt.0,4 2 • 4 • 50 • 4802 2 • 4 • 50 • 4802 hoptflA — 0,5326 = 0,0034 Tato hodnota je příliš velká na to, aby byla optimální hodnotou vyhlazovacího parametru pro hustotu, která je definována na intervalu [0,1]. Ve výpočtu samotném není chyba, avšak byl porušen předpoklad o spojitosti derivací funkce až do řádu 4 včetně - viz větu 3.2. Už první derivace této funkce hustoty není spojitá (nakreslete si její graf). Cv. 3 Postupujeme podobně jako v ukázkovém příkladě na str. 36, /iOpt,0,2 — 0,9221 • n-1/5. Cv. 4 a) a2 = J\ x2gk(x) dx, gk(x) = Ak(í - x2)k+1 a Ak = {k^%+3 pak < = Ak J x2(l-x2)k+1dx = Ak J x-x(l -x2)k+íd. v! = 1 v' — x(í — x2) 2^fe+l = Ak Ai -1 2(fe+2) 1 2(k + 2) (í-x2) 2\fc+2 1_C1- r2)k+2 1 + Ak 1 2(k + 2) (í-x2)k+2dx 2(k + 2) J_i (1 - x2)k+2 dx Víme, že Odtud plyne, že a tedy Ak+1(í - xz)k+z dx — 1 (gk+i je hustota). =9k+i(x) (1 - x2)k+2 dx Ak+i' = Ak 1 1 2(fe + 2) Ak+1 (2fc + 3)! 1 ((fc + 1) + l)!222(fc+l)+3 (/j + 1)!222fe+3 2(k + 2) (2(Jfc+l) + 3)! (2fc + 3)!(fc + 2)!222fc+5 _ 1 (k + l)!222fe+32(/j + 2)(2fc + 5)! ~ 2/j + 5' 81 b) —x] —g(—x)dx — a ) a„ \ a x2—g (—x) dx -a \ a ) í t2—g(t)—dt =% I t2g(t) dt = a J 1. >0 = £ Er=i = § + i Er=i Vlastnosti T^(x): x + 1 1 T^2(x)dx = T (X + 1)2 dx = - (x+1)3 dx x + 1 , 1 -dx = — 2 2 (x + 1)2 Cv. 1 Podle definice stačí integrovat polynom, který odpovídá kvartickému jádru, tj. i|(l-í2)2dt= 1| i (1_2í2+í4)dí -i 16 16 7_! 15 16 2 „ 1 . t — -t3 + -í5 15 (x - ^x3 + \xb + ^) = 7^(3x5 - 10x3 + 15x + í 16 15' 16v Cv. 2 • Vlastnost 2: platí W'(x) = K(x), W(x) dx v! — 1 u — x xW^x)]1^ - / xif(x) dx = 1 • VK(1) - (-1) • W(-í) = 1. =o 82 Protože 0 < W(x) < 1, pak platí j\ W2(x) dx < f\ W(x) dx = 1. • Vlastnost 3: W(x)K(x) dx = u = W(x) u' = K(x) v' = K(x) v = = [^(x)]1!! - / ET(x)VK(x) dx, = 1 a W{-1) = 0 VK(x)ET(x)dx = -. =T xW(x)K(x) dx u — x u' — 1 í/ = K{x)W{x) v = ±VK2(x) [^^2W]-i - f \w2(x) te=\-\ £ W2(x) dx 1- / W2{x)dx C v. 3 Vztah (4.6) pro AMISE zderivujeme vzhledem k h, položíme roven nule a vypočítáme h: dAMISEF(-,/i) 1 t ,3 n ah, n Toto vypočítané h dosadíme do rovnice (4.6) a upravíme. Cv. 4 Kvartické jádro: K(x) = ±§(1 - x2)2, fi2(K) = |a/'j Wr2(x)dx = 0,7835. Dosadíme do vztahu (4.9) a dostaname 1/3/71a^^y/3^.^4;5089.,.n-1/3. V 15 • 4g / Kapitola 5 Cv. 1 Hodnota c* — J J K\ (\/x2 + y2) dxdy pro kvartické jádro je ir/3. Cv. 2 Zderivujeme vztah (5.2) pro AMISE a dostaneme dAMISE/(-,H) V{K) 1 2 0 = 0 = d/ii nh\h2 ' 4 dAMISE/(-,H) 1 2 + -B22{K)(Ah\V{fxx) + Ah^lViUy)) + -B22{K)(Ah\h2V{fxy)+AhlV{fyy)) dh2 nhih2 4 První rovnici vynásobíme hlr druhou rovnici vynásobíme h2 a odečteme je 0 = \ŕ2(K)(4hjV(fxx) +4h2h2V(fxy)) - \r2(K)(4h2h2V(fxy) + 4h42V(fyy)), odtud plyne ±h\V{fxx) + 4h2h2V(fxy) = 4h2h2V(fxy) + 4h4V(fyy), hj ^ V(fyy) h\ v(fxxy Nyní dosadíme do rovnice (6.1) h\ — h2- X\fV\> vypočítáme h2 a pak h±. (6.1) (6.2) 83 Cv. 3 Rovnici (6.1) vynásobíme h±, rovnici (6.2) vynásobíme h2 a rovnice sečteme ^jffi = 4 • \pl{K) (h\V{.fxx) + 2h\h\V{]xy) + h42V(fyy)) 2-AI V Odtud už plyne tvrzení. Cv. 5 Ověříme, že hlavní minory jsou kladné: AISB Mu 6 > 0, M2 6 -1 -1 4 23 > 0, M- 33 -1 4 -3 2-3 9 84 Kapitola 7 Datové soubory Tabulka 7.1: Hodnoty simulovaných dat z příkladu 2.1 - regrese X y X y x y x y x y 0,01 0,3002 0,02 0,9792 0,03 -1,0372 0,04 0,5519 0,05 0,3070 0,06 -0,4816 0,07 -0,0225 0,08 0,3848 0,09 2,0187 0,10 1,6268 0,11 -0,4240 0,12 1,7734 0,13 0,6198 0,14 0,2228 0,15 0,6049 0,16 0,1343 0,17 0,1602 0,18 0,9489 0,19 0,8871 0,20 0,8664 0,21 0,4661 0,22 -0,5034 0,23 0,4270 0,24 0,8499 0,25 0,2444 0,26 0,4820 0,27 0,2926 0,28 -0,2577 0,29 0,0068 0,30 -0,5664 0,31 0,2407 0,32 -0,8052 0,33 -0,7914 0,34 -0,6835 0,35 -1,7689 0,36 0,4086 0,37 -0,1572 0,38 -0,7018 0,39 0,3614 0,40 -1,1739 0,41 -0,3584 0,42 -0,4120 0,43 -0,1105 0,44 -0,0875 0,45 -0,6454 0,46 -0,1926 0,47 -0,2206 0,48 0,2188 0,49 0,4979 0,50 0,5546 0,51 -0,3811 0,52 0,1413 0,53 -0,4521 0,54 -0,3496 0,55 0,2547 0,56 1,0735 0,57 -0,0313 0,58 0,5820 0,59 0,3219 0,60 1,0265 0,61 -0,0495 0,62 0,5318 0,63 0,8049 0,64 1,0842 0,65 1,3028 0,66 0,5615 0,67 -0,2485 0,68 0,0955 0,69 -0,1043 0,70 1,5522 0,71 0,0104 0,72 0,6246 0,73 0,0784 0,74 0,5351 0,75 -0,3824 0,76 -0,7979 0,77 -0,9098 0,78 -0,0599 0,79 -0,5005 0,80 -0,6184 0,81 0,0816 0,82 -0,5874 0,83 -0,7349 0,84 -0,1342 0,85 -1,4160 0,86 -0,7407 0,87 -0,7340 0,88 -1,3214 0,89 -1,1229 0,90 -1,8261 0,91 -1,8089 0,92 -1,1517 0,93 -0,7882 0,94 0,2239 0,95 -1,2950 0,96 -0,7328 0,97 -0,7041 0,98 -1,4370 0,99 -0,4688 1,00 -0,8973 85 Tabulka 7.2: Hodnoty simulovaných dat z ukázkového příkladu 2.2 - regrese X y X y x y x y x y 0 0,1848 0,0204 0,4321 0,0408 -0,0322 0,0612 -0,4980 0,0816 -0,1456 0,1020 0,4379 0,1224 -0,1272 0,1429 0,4285 0,1633 0,2718 0,1837 0,6568 0,2041 -0,1325 0,2245 0,3708 0,2449 0,1820 0,2653 1,2750 0,2857 0,8176 0,3061 1,0174 0,3265 0,4667 0,3469 0,6689 0,3673 0,7680 0,3878 1,1553 0,4082 0,8496 0,4286 1,1259 0,4490 0,4616 0,4694 0,9023 0,4898 0,7931 0,5102 0,6047 0,5306 1,1178 0,5510 1,0450 0,5714 0,9589 0,5918 0,5856 0,6122 1,1626 0,6327 0,9237 0,6531 0,7113 0,6735 0,7370 0,6939 0,6072 0,7143 0,1737 0,7347 0,4766 0,7551 0,2761 0,7755 0,1754 0,7959 0,0686 0,8163 0,1642 0,8367 -0,2599 0,8571 0,4293 0,8776 0,2708 0,8980 0,0943 0,9184 0,0556 0,9388 -0,1630 0,9592 0,2710 0,9796 -0,0292 1,0000 -0,1786 Tabulka 7.3: Hodnoty simulovaných dat z příkladu 3.1 - hustota 0,0916 -0,5149 0,4746 0,1535 0,0676 0,2576 0,1307 -0,4707 -0,0812 -0,0730 -0,2660 0,8411 -0,4379 -0,2419 -0,3560 -0,5871 -0,0961 -0,1370 0,7650 -0,1245 -0,5321 0,8017 0,6173 -0,1148 -0,7531 -0,2223 -0,0780 0,1380 -0,1306 0,2217 2,1959 1,3747 1,5260 1,6294 1,7461 1,8397 2,0062 0,4854 1,7715 2,6212 1,4666 2,4669 2,1752 1,9855 2,0912 1,2175 1,9577 2,8020 2,0492 2,0207 1,6329 1,9846 2,1162 2,2132 1,8136 1,8818 3,0118 0,8708 3,1147 2,1688 2,5000 1,1679 1,7050 1,8610 2,2114 1,1649 2,2358 1,3936 2,0331 2,3262 2,1635 2,5413 2,5030 1,6745 2,1285 1,5278 1,3391 2,4624 2,0000 1,9725 2,4556 2,2973 2,1751 2,6251 2,4649 2,1199 1,6548 1,6742 2,5961 1,1941 1,9878 1,0256 2,5102 2,4309 2,0006 1,9646 0,7569 2,2906 0,9038 0,8404 86 Tabulka 7.4: Hodnoty simulovaných dat z ukázkového příkladu 3.2 a 4.3 - hustota/distribuční funkce_ 0,3636 0,5101 0,1957 -0,6054 0,5500 -0,5962 -0,0499 0,1363 -0,4110 -0,0015 0,0054 0,1886 0,3520 -0,2376 0,1681 0,2509 -0,2134 -0,5599 -0,0562 -0,1096 0,1991 0,2475 0,5774 -0,1492 0,4272 0,5760 0,2348 -0,0995 -0,4479 0,3291 0,2596 0,2940 -0,2144 0,1691 0,0784 -0,2717 0,0432 -0,3811 0,1452 -0,3361 0,1843 -0,1722 -0,3160 -0,0094 0,1448 -0,3475 0,2437 -0,2367 -0,3658 -0,2341 Tabulka 7.5: Hodnoty simulovaných dat z příkladu 4.1 - distribuční funkce 0,5603 -0,5920 0,0667 -0,3630 0,2023 -0,4002 0,3266 0,4929 1,1413 -0,1294 -1,4256 -0,5597 0,9031 -0,7148 0,6406 0,0827 0,9578 -1,3073 -0,1318 -0,8052 1,9387 0,5501 0,9193 -0,7055 -0,3124 -0,1816 0,7323 -0,1852 0,4677 -1,3679 -0,2359 -0,5491 -1,0514 0,3386 0,1880 0,0223 -0,8891 0,7517 0,2335 -0,1994 0,0153 -0,1747 -1,1668 -0,1904 -0,5542 -0,6528 -0,7709 -0,3557 -1,3351 0,6428 2,5201 1,9800 1,9652 1,2018 3,0187 1,8668 1,2855 3,3514 1,7752 1,4110 1,7062 1,1521 0,8799 4,5260 3,6555 2,3075 0,7429 1,1345 1,8235 2,7914 0,6680 -0,3299 0,5509 2,3335 2,3914 2,4517 1,8697 2,1837 1,5238 2,8620 0,6383 2,4550 1,1513 1,6651 2,5528 3,0391 0,8824 3,2607 2,6601 1,9321 1,8048 1,7824 1,6969 2,0230 2,0513 2,8261 3,5270 2,4669 1,7903 2,6252 87 Tabulka 7.6: Hodnoty simulovaných dat z příkladu 5.1 - dvourozměrná hustota [x;y] [x;y] [x;y] [-0,9466; -0,9065] [0,3206; 3,7036] [-0,1005; -0,8364] [0,3740; 4,6725] [-1,2990; 0,2211] [1,5530; 2,2048] [-0,0255; 3,1416] [-0,6542; 1,8057] [-0,1060;-0,5221] [-0,4199; 2,2447] [1,0924; 2,4151] [0,0067; 3,2366] [-0,2661; 3,7164] [0,5607; 3,8076] [-2,5687; 1,1808] [0,3246; 3,6851] [-0,2423; -0,4698] [-0,3203;-0,5281] [0,3535; 4,0112] [0,5172; 1,0718] [0,4428; 3,0552] [1,2592;-0,9592] [-0,2283;-1,5876] [0,3677; -0,3296] [-0,4947; 3,9449] [0,3992; -0,5683] [0,4789; 4,7000] [0,8149; 3,9766] [0,7421; 1,8418] [0,5522; -0,5483] [-0,7812; 3,5773] [0,9459; -0,4850] [-0,2062; 4,6773] [-0,6864; 3,8091] [-0,4114;-1,9323] [1,3919; 4,3463] [-0,5084; -0,7163] [-0,0837;-1,6576] [0,3252; 4,4022] [0,5677;-1,0168] [1,0639; 2,7468] [-0,4263; -0,2424] [0,0356;-1,6677] [-0,3552; -0,4372] [0,8428; 1,3214] [-0,6772; 3,9948] [0,0065;-0,3730] [-1,1498;-0,6113] [0,9380; 4,2303] [-1,0416; 4,7185] [0,0741;-1,8323] [0,4332; 2,3411] [1,6562; 2,5772] [0,2473; 1,6168] [1,0169;-0,8521] [-0,1442;-1,4423] [-0,9036; 2,6207] [0,2175; 1,8817] [1,2094; 2,2929] [-0,1639;-0,8925] [-0,2548; 4,5990] [-0,0761;-0,3419] [1,6432; 4,2209] [-0,2305;-1,1597] [-0,5887; 3,8457] [-0,8657; 4,2465] [-0,4201; 4,1738] [0,7746; 1,5665] [0,0917; 2,5265] [-0,1711; 4,0729] [0,0750; 3,5324] [-0,4325; 2,2264] [-0,8325;-1,7409] [0,4678; 1,2271] [-0,4280; -0,6307] [0,5144;-1,3583] [-1,0391; 2,6486] [1,5309; 1,6144] [0,0852; 3,9299] [0,2001; 2,7036] [2,0835; 2,9561] [-0,6129; 3,7847] [-0,6465; 1,7759] [-0,2576; 1,2650] [-0,5027; -0,9906] [0,4180;-1,0515] [0,9001; 4,4942] [0,7250; 1,6023] [0,0828; 3,0127] [0,2756; 2,8724] [0,4105; 3,7973] [0,3790;-1,5109] [-0,5589;-1,1707] [-0,8010; 3,9085] [0,5950; 4,4687] [0,5173; 3,0249] [-2,0882; 3,6987] [-0,8476; 3,9186] [0,2027; 3,6723] [0,2096;-0,3292] 88 Tabulka 7.7: Hodnoty simulovaných dat z ukázkového příkladu 5.3 - dvourozměrná hustota [x;v] [0,7747; -0,2338] [-0,5707; 0,3048] [1,0496;-0,5814] [-0,2162; 0,6287] [-0,8382; 1,0305] [-0,3572;-0,1754] [0,5395;-0,0050] [0,6096; 0,2925] [0,0925; 0,6947] [0,1659; 0,3014] [0,3654; 0,2292] [-0,8424; -0,2788] [-0,3359; -0,1337] [-1,0566; 0,3541] [-0,2593; 0,9716] [0,5481; 0,3282] [0,4142;-0,6963] [1,3209; 0,2500] [0,0963;-0,0532] [-1,2451;-0,1409] [0,9754;-0,1265] [0,7403; -0,5079] [-0,3364; -0,1852] [0,1285;-0,4529] [0,4500;-0,2107] [0,2221; 0,0677] [0,4020; 0,1007] [-0,0651;-0,1245] [-0,0224; -0,0881] [0,0509; 0,0194] [-0,4559; 0,2646] [0,8215; 0,7400] [0,0037; 0,3580] [0,3963;-0,1655] [-0,4936; 0,3453] [-0,6208;-0,1418] [-0,0525;-1,0281] [-0,1459;-1,3552] [0,6898;-0,9893] [0,3499;-0,1369] [0,3121;-0,3577] [-0,8342;-0,4158] [-0,5587; 0,0235] [0,1167;-0,3805] [-0,2008; -0,3640] [0,7448; 0,3520] [0,7743; 0,4517] [-0,3353; 1,4068] [0,5053; 0,1746] [0,8057; 0,0337] [0,1853;-0,3315] [0,9373;-0,1127] [0,6134;-0,0292] [0,1236;-1,0137] [0,5062;-0,1486] [0,1122; 0,4327] [-0,2395; -0,3295] [0,2070; 0,2853] [0,0209; 0,5724] [-0,2105; 0,6219] [1,0225; 0,4485] [0,1278;-1,3047] [1,5764;-0,3258] [-0,1710; 0,5205] [0,3258; 0,2279] [0,9137; 0,0028] [1,5757; 0,2523] [0,1753; 0,6524] [-0,4563; 0,0888] [0,4638; -0,4380] [-0,7620; -0,1040] [-0,1581;-0,1200] [0,1570; 0,5177] [0,8600; 0,2380] [0,8297; 0,3161] [-0,0159; 0,0854] [0,5334; -0,2092] [-0,1855; 0,0353] [0,3817; 0,1251] [1,1058;-0,8072] [0,0905; 0,1033] [-0,1735; 0,5355] [-0,2887; -0,0636] [0,4137; 0,9078] [0,7721; 0,3236] [1,2273;-0,2427] [0,0312;-0,0889] [0,0941;-0,1589] [-0,3653; 0,1392] [0,9307; 0,4357] [-0,0436; -0,1685] [-0,2732; 0,2647] [1,2611; 1,0134] [0,4676; 0,5378] [0,4136; 0,4765] [-0,5806; -0,4704] [0,2506; -0,3084] [0,3003; 0,1400] [-0,6306;-0,0581] [-0,1304;-0,9382] 89 Tabulka 7.8: Hodnoty reálných dat z kapitoly 2, odstavce 7 - Hurónské jezero X V X V X V X V X V 1875 10,38 1876 11,86 1877 10,97 1878 10,80 1879 9,79 1880 10,39 1881 10,42 1882 10,82 1883 11,40 1884 11,32 1885 11,44 1886 11,68 1887 11,17 1888 10,53 1889 10,01 1890 9,91 1891 9,14 1892 9,16 1893 9,55 1894 9,67 1895 8,44 1896 8,24 1897 9,10 1898 9,09 1899 9,35 1900 8,82 1901 9,32 1902 9,01 1903 9,00 1904 9,80 1905 9,83 1906 9,72 1907 9,89 1908 10,01 1909 9,37 1910 8,69 1911 8,19 1912 8,67 1913 9,55 1914 8,92 1915 8,09 1916 9,37 1917 10,13 1918 10,14 1919 9,51 1920 9,24 1921 8,66 1922 8,86 1923 8,05 1924 7,79 1925 6,75 1926 6,75 1927 7,82 1928 8,64 1929 10,58 1930 9,48 1931 7,38 1932 6,90 1933 6,94 1934 6,24 1935 6,84 1936 6,85 1937 6,90 1938 7,79 1939 8,18 1940 7,51 1941 7,23 1942 8,42 1943 9,61 1944 9,05 1945 9,26 1946 9,22 1947 9,38 1948 9,10 1949 7,95 1950 8,12 1951 9,75 1952 10,85 1953 10,41 1954 9,96 1955 9,61 1956 8,76 1957 8,18 1958 7,21 1959 7,13 1960 9,10 1961 8,25 1962 7,91 1963 6,89 1964 5,96 1965 6,80 1966 7,68 1967 8,38 1968 8,52 1969 9,74 1970 9,31 1971 9,89 1972 9,96 Tato data jsou přístupná i v programu r, kde jsou původní hodnoty úrovně hladiny. V této tabulce je hodnota úrovně hladiny snížena o 570 stop. 90 Tabulka 7.9: Hodnoty reálných dat z kapitoly 2, odstavce 7 - krystaly ledu X V X V X V X V X V 50 19 60 20 60 21 70 17 70 22 80 25 80 28 90 21 90 25 90 31 95 25 100 29 100 30 100 33 105 32 105 35 110 28 110 30 110 30 115 30 115 31 115 36 120 25 120 28 120 36 125 28 130 31 130 32 135 25 135 34 140 26 140 33 145 31 150 33 150 36 155 33 155 41 160 30 160 37 160 40 165 32 170 35 180 38 Tabulka 7.10: Hodnoty reálných dat z kapitoly 3, odstavce 7 - krabi 16,1 18,1 19,0 20,1 20,3 23,0 23,8 24,5 24,2 25,2 27,3 26,8 27,7 27,2 27,4 26,8 28,2 28,3 27,8 29,2 31,3 31,9 31,4 32,4 32,5 32,3 33,0 35,8 34,0 33,8 34,9 36,0 35,6 35,7 38,1 36,2 37,3 36,4 36,7 37,6 38,7 39,7 39,2 42,1 41,6 40,9 41,9 43,2 42,4 47,1 14,7 19,3 18,5 19,2 19,6 20,4 20,9 21,3 21,7 22,5 22,5 22,8 24,7 24,6 23,7 24,9 26,0 24,6 25,4 26,1 27,1 26,7 27,9 27,3 27,6 27,9 28,4 28,6 30,0 30,1 30,1 31,7 32,8 31,8 31,9 31,7 33,9 32,6 32,4 33,4 32,8 33,9 33,6 34,5 34,5 34,2 36,6 38,2 38,6 40,9 16,7 20,2 20,7 22,7 23,2 24,2 26,0 27,1 26,6 27,5 29,2 28,9 29,1 28,7 28,7 27,8 29,2 29,9 29,0 30,2 30,9 30,2 31,7 32,3 31,6 35,0 36,1 34,4 34,6 36,0 36,9 36,7 38,8 37,9 37,8 36,9 37,2 39,2 39,1 39,8 40,6 42,8 42,9 45,5 45,7 43,4 45,4 44,6 47,2 47,6 21,4 21,7 24,1 25,0 25,8 27,0 28,8 28,1 29,6 30,0 30,1 31,2 31,6 31,0 31,0 31,6 31,4 31,6 32,3 33,1 34,5 34,5 33,3 34,0 34,7 37,9 35,1 35,6 36,5 37,0 34,7 35,8 36,3 37,8 37,9 39,9 39,4 40,1 40,4 39,8 39,4 40,0 41,5 39,9 43,8 41,2 41,7 42,6 43,0 46,2 91 Tabulka 7.11: Hodnoty reálných dat z kapitoly 3, odstavce 7 - cholesterol 184 215 221 210 208 197 250 180 212 297 168 208 180 268 219 319 250 285 221 227 224 172 181 215 179 245 193 242 172 262 243 211 219 173 308 249 294 266 169 260 267 270 213 131 218 225 263 233 131 251 284 216 243 208 193 232 197 220 254 248 159 171 196 184 204 197 209 174 191 228 218 191 332 175 190 211 261 249 233 260 227 258 167 217 204 199 228 188 178 233 194 280 185 212 211 175 231 230 175 386 230 150 417 191 191 245 200 194 298 228 276 196 223 192 185 245 279 207 194 138 144 178 185 209 220 258 168 194 208 249 184 207 187 160 172 269 252 185 271 221 232 185 171 265 200 236 169 239 172 119 176 171 233 244 306 171 165 193 278 221 206 186 234 248 195 244 194 331 171 177 348 131 178 140 208 218 206 206 304 218 198 170 184 163 173 239 313 184 258 197 240 230 181 178 240 171 283 239 232 236 175 229 211 211 251 283 210 242 264 139 243 206 105 235 222 165 194 168 164 187 185 245 198 210 140 257 222 149 203 216 230 168 240 198 164 230 185 188 189 242 191 179 253 196 189 260 251 195 264 185 140 178 226 201 237 246 271 191 201 267 231 299 230 208 151 171 159 242 242 229 209 259 238 194 239 222 231 176 198 230 233 213 200 180 323 217 208 220 247 237 254 256 214 245 157 197 185 219 239 162 247 197 223 193 227 258 274 250 287 165 221 222 262 189 232 198 139 273 142 232 195 170 234 158 219 155 243 168 237 150 222 167 266 207 209 207 205 200 116 217 190 238 221 201 228 157 234 156 168 178 190 169 194 190 187 210 289 180 178 130 178 149 208 201 193 251 206 265 147 160 168 92 Tabulka 7.12: Hodnoty reálných dat z kapitoly 4, odstavce 6 - pyrimidin 0,571 0,900 0,639 0,100 0,763 0,619 0,897 0,745 0,579 0,642 0,561 0,763 0,584 0,602 0,589 0,621 0,700 0,716 0,772 0,805 0,833 0,582 0,547 0,568 0,613 0,619 0,560 0,584 0,720 0,619 0,624 0,534 0,628 0,595 0,628 0,634 0,717 0,734 0,587 0,549 0,516 0,900 0,859 0,540 0,900 0,893 0,632 0,451 0,554 0,628 0,646 0,545 0,649 0,661 0,741 0,749 0,742 0,634 0,538 0,531 0,893 0,838 0,674 0,569 0,572 0,738 0,638 0,829 0,675 0,568 0,665 0,671 0,753 0,756 93 Tabulka 7.13: Hodnoty reálných dat z kapitoly 4, odstavce 6 - Dyje - část 1. 16,9 13,2 11,4 15,4 40,4 18,9 5,44 2,15 1,96 2,20 5,74 6,15 8,69 10,4 23,2 3,53 5,66 5,68 3,55 5,59 3,14 3,02 0,937 1,49 4,33 4,90 5,41 5,00 5,03 5,13 14,4 15,1 11,1 11,7 6,34 4,82 7,02 7,86 5,61 8,87 33,1 14,5 8,09 4,79 5,06 5,28 9,12 12,0 8,14 10,2 12,6 14,1 7,34 6,14 5,00 4,75 4,82 8,31 25,5 6,74 10,9 11,0 16,6 18,0 19,4 12,4 20,7 29,4 14,7 17,8 15,3 17,4 24,1 49,4 16,4 18,4 36,4 31,8 29,5 12,1 15,6 12,0 8,29 7,87 7,48 7,24 9,30 14,5 64,7 69,2 26,0 29,9 16,5 13,9 13,7 8,62 22,9 13,0 20,1 22,0 15,0 45,0 16,2 14,4 10,9 9,43 11,5 6,06 5,17 5,14 6,09 6,96 6,13 5,42 5,61 5,93 5,97 6,56 5,28 4,07 4,24 4,34 4,09 4,98 8,58 31,1 14,3 16,8 12,1 8,91 7,87 6,48 9,65 17,6 15,2 21,7 29,6 16,4 8,52 6,23 4,75 3,82 4,09 4,87 7,26 7,56 7,31 20,2 29,1 5,64 4,03 4,25 4,10 4,38 7,11 5,43 6,71 7,09 8,84 16,1 39,4 27,7 8,03 7,22 5,35 4,95 5,37 3,98 3,89 5,38 20,9 50,5 20,9 9,01 4,90 4,90 4,53 5,57 5,39 9,48 9,31 5,34 5,10 3,97 3,10 2,49 1,96 2,69 2,67 8,43 10,4 12,8 8,67 4,95 7,11 6,07 5,85 7,34 6,66 8,69 6,57 3,46 1,59 1,92 5,56 5,27 10,2 15,8 12,2 6,41 28,4 6,89 8,04 5,36 9,47 11,2 14,0 3,94 3,75 6,55 19,4 14,2 11,3 6,15 4,28 3,65 3,87 4,10 4,84 5,34 15,4 9,31 6,93 4,87 3,69 5,21 7,20 10,3 9,35 5,14 4,87 4,94 4,87 5,06 4,39 4,06 4,78 4,94 10,6 4,67 5,14 5,24 7,29 9,40 23,90 10,50 18,40 25,20 9,01 5,57 6,01 15,9 9,10 5,90 5,88 6,26 9,41 20,1 16,8 13,0 17,9 8,87 6,81 6,38 7,58 7,36 9,97 8,72 9,18 14,8 12,6 13,0 7,66 6,30 5,72 6,38 7,82 7,43 8,51 8,91 7,94 6,69 9,70 18,4 8,60 6,84 5,87 5,48 6,97 5,99 8,55 10,9 5,92 7,49 6,71 5,87 5,83 6,55 9,04 11,8 8,11 8,11 6,58 6,15 7,59 6,35 20,9 12,8 4,97 9,47 6,22 25,0 9,99 13,1 11,3 8,80 10,5 18,6 11,3 8,41 8,02 24,0 9,17 7,40 7,15 8,41 6,92 9,77 7,67 25,2 94 Tabulka 7.14: Hodnoty reálných dat z kapitoly 4, odstavce 6 - Dyje - část 2. 7,60 6,09 24,2 11,8 7,63 6,52 7,88 8,27 6,56 9,86 5,25 3,92 5,32 21,6 11,3 8,90 6,84 4,53 4,71 4,74 4,78 4,53 4,67 4,64 4,57 5,04 5,49 5,02 4,82 4,70 4,58 5,17 5,71 6,70 7,88 7,00 34,7 37,2 37,5 69,9 18,4 7,79 5,37 4,65 6,16 6,20 10,0 31,2 8,34 7,36 5,33 5,45 5,42 11,4 22,2 7,26 5,90 14,3 22,4 23,0 24,8 14,3 5,55 27,8 7,97 4,40 4,62 5,34 5,05 4,58 17,5 16,7 12,7 9,37 5,05 6,10 4,61 4,17 5,26 5,84 5,68 7,01 5,65 5,56 23,7 18,9 7,60 8,16 6,12 4,42 4,33 4,20 4,40 11,5 5,28 4,70 17,5 27,5 7,72 6,47 5,48 5,23 5,14 5,19 4,89 10,5 8,18 8,68 8,00 8,05 6,82 8,88 6,67 5,71 5,74 4,86 5,13 5,07 5,45 4,90 4,87 9,28 24,5 9,70 10,0 7,18 4,17 4,55 4,71 4,66 4,49 4,18 4,67 12,0 8,80 5,61 4,40 4,70 4,70 4,08 3,60 3,92 3,69 3,62 3,60 4,67 4,14 3,48 8,85 5,65 3,88 4,78 4,10 27,6 18,4 9,33 3,64 8,10 8,28 4,62 19,9 5,55 4,38 3,66 4,07 4,10 32,5 11,1 10,4 6,39 5,21 5,58 5,23 4,61 4,69 4,35 5,14 15,5 8,87 41,5 38,9 14,0 8,39 6,73 6,40 5,81 5,35 4,96 4,91 4,16 4,65 6,23 15,5 6,79 6,93 6,09 5,54 6,03 6,51 5,11 3,82 3,85 4,42 4,00 9,98 27,0 11,0 6,18 7,21 6,77 5,30 6,65 5,76 21,6 9,67 19,2 4,55 26,7 18,9 7,75 6,29 9,09 8,43 7,83 6,80 4,55 4,66 11,2 16,4 6,72 4,97 5,27 4,59 5,81 4,11 7,52 12,5 22,5 19,5 22,8 17,2 14,6 6,88 5,96 5,69 6,64 6,05 5,57 5,48 4,24 7,19 16,7 17,7 23,3 9,60 9,37 8,98 7,41 4,96 3,90 3,47 2,85 3,03 3,67 2,95 6,80 10,4 5,67 5,50 5,66 5,39 4,91 5,14 5,98 5,69 9,23 22,5 20,2 20,7 16,9 8,35 28,9 8,81 6,25 5,02 16,7 22,7 12,7 17,7 14,0 8,13 13,6 11,1 8,27 7,68 5,41 4,73 4,46 29,8 25,2 39,3 42,9 24,9 29,5 24,5 10,5 7,57 7,90 8,77 12,5 10,6 9,64 44,2 29,5 8,45 8,42 8,02 7,04 5,83 4,13 3,83 4,99 12,7 5,56 6,72 6,44 10,1 7,63 8,41 7,66 7,28 6,14 4,79 3,64 2,85 3,15 3,93 7,51 6,70 7,76 6,60 4,88 3,44 3,77 95 Tabulka 7.15: Hodnoty reálných dat z kapitoly 5, odstavce 6 - studie UNICEF [x;y] [x;y] [x;y] [x;y] [x;y] [x;y] [19; 73] [235; 53] [43; 69] [257; 43] [109; 56] [225; 48] [205; 52] [28; 71] [99; 67] [77; 63] [38; 66] [39; 71] [143; 42] [72; 63] [19; 72] [20; 68] [108; 51] [153; 51] [45; 67] [105; 72] [95; 62] [175; 48] [29; 73] [24; 69] [94; 57] [155; 50] [35; 73] [68; 69] [132; 44] [260; 45] [123; 53] [123; 43] [138; 54] [93; 64] [107; 61] [109; 60] [38; 69] [76; 63] [169; 48] [32; 67] [79; 60] [77; 60] [158; 54] [183; 52] [122; 50] [107; 56] [126; 47] [202; 42] [100; 54] [100; 57] [183; 52] [61; 68] [124; 45] [138; 56] [141; 52] [165; 51] [180; 44] [136; 53] [91; 59] [183; 40] [197; 47] [197; 39] [231; 52] [200; 46] [111; 52] [72; 68] [183; 40] [265; 46] [211; 45] [316; 40] [172; 44] [190; 41] 96 Tabulka 7.16: Hodnoty reálných dat z kapitoly 5, odstavce 6 - koncentrace lipidů - část 1. [184; 145] [215; 168] [221; 432] [210; 92] [208; 112] [197; 87] [250; 118] [180; 80] [212; 130] [297; 232] [168; 208] [208; 262] [180; 102] [268; 154] [219; 454] [319; 418] [250; 161] [285; 930] [221; 268] [227; 146] [224; 124] [172; 106] [181; 119] [215; 325] [179; 126] [245; 166] [193; 290] [242; 179] [172; 207] [262; 88] [243; 126] [211; 306] [219; 163] [173; 300] [308; 260] [249; 146] [294; 135] [266; 164] [169; 158] [260; 98] [267; 192] [270; 110] [213; 261] [131; 96] [218; 567] [225; 240] [263; 142] [233; 340] [131; 137] [251; 189] [284; 245] [216; 112] [243; 50] [208; 220] [193; 188] [232; 328] [197; 291] [220; 75] [254; 153] [248; 312] [159; 125] [171; 78] [196; 130] [184; 255] [204; 150] [197; 265] [209; 82] [174; 117] [191; 233] [228; 130] [218; 123] [191; 90] [332; 250] [175; 246] [190; 120] [211; 304] [261; 174] [249; 256] [233; 101] [260; 127] [227; 172] [258; 145] [167; 80] [217; 227] [204; 84] [199; 153] [228; 149] [188; 148] [178; 125] [233; 141] [194; 278] [280; 218] [185; 115] [212; 171] [211; 124] [175; 148] [231; 181] [230; 90] [175; 489] [386; 162] [230; 158] [150; 426] [417; 198] [191; 115] [191; 136] [245; 120] [200; 152] [194; 183] [298; 143] [228; 142] [276; 199] [196; 103] [223; 80] [192; 101] [185; 130] [245; 257] [279; 317] [207; 316] [194; 116] [138; 91] [144; 125] [178; 84] [185; 100] [209; 89] [220; 153] [258; 151] [168; 126] [194; 196] [208; 201] [249; 200] [184; 182] [207; 150] [187; 115] [160; 125] [172; 146] [269; 84] [252; 233] [185; 110] [271; 128] [221; 140] [232; 258] [185; 256] [171; 165] [265; 156] [200; 68] [236; 152] [169; 112] [239; 154] [172; 140] [119; 84] [176; 217] [171; 108] [233; 127] [244; 108] [306; 408] [171; 120] [165; 121] [193; 170] [278; 152] [221; 179] 97 Tabulka 7.17: Hodnoty reálných dat z kapitoly 5, odstavce 6 - koncentrace lipidů - část 2. [x;y] [x;y] [x;y] [206; 133] [186; 273] [234; 135] [248; 142] [195; 363] [244; 177] [194; 125] [331; 134] [171; 90] [177; 133] [348; 154] [131; 61] [178; 101] [140; 99] [208; 148] [218; 207] [206; 148] [206; 107] [304; 149] [218; 96] [198; 103] [170; 284] [184; 184] [163; 156] [173; 56] [239; 97] [313; 256] [184; 222] [258; 210] [197; 158] [240; 196] [230; 162] [181; 104] [178; 100] [240; 441] [171; 170] [283; 424] [239; 92] [232; 131] [236; 148] [175; 153] [229; 242] [211; 91] [211; 122] [251; 152] [283; 199] [210; 217] [242; 85] [264; 269] [139; 173] [243; 112] [206; 201] [105; 36] [235; 144] [222; 229] [165; 151] [194; 400] [168; 91] [164; 80] [187; 390] [185; 231] [245; 322] [198; 124] [210; 95] [140; 102] [257; 402] [222; 348] [149; 237] [203; 170] [216; 101] [230; 304] [168; 131] [240; 221] [198; 149] [164; 76] [230; 146] [185; 116] [188; 220] [189; 84] [242; 144] [191; 115] [179; 126] [253; 222] [196; 141] [189; 135] [260; 144] [251; 117] [195; 137] [264; 259] [185; 120] [140; 164] [178; 109] [226; 72] [201; 297] [237; 88] [246; 87] [271; 89] [191; 149] [201; 92] [267; 199] [231; 161] [299; 93] [230; 137] [208; 77] [151; 73] [171; 135] [159; 82] [242; 180] [242; 248] [229; 296] [209; 376] [259; 240] [238; 156] [194; 272] [239; 38] [222; 151] [231; 145] [176; 166] [198; 333] [230; 492] [233; 142] [213; 130] [200; 101] [180; 202] [323; 196] [217; 327] [208; 149] [220; 172] [247; 137] [237; 400] [254; 170] [256; 271] [214; 223] [245; 446] [157; 59] [197; 101] [185; 168] [219; 267] [239; 137] [162; 91] [247; 91] [197; 347] [223; 154] [193; 259] [227; 202] [258; 328] [274; 323] [250; 160] [287; 209] [165; 155] [221; 156] [222; 108] [262; 169] [189; 176] [232; 583] [198; 105] [139; 54] [273; 146] [142; 111] [232; 161] 98 Kapitola 8 Přílohy Tabulka 8.1: Optimální jádra pro v — 0 a v — 1 i/ = 0 k K