Ú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 fa f F I K(-) V (K) V(g) PÁK) f*g w(-) Ck[0,l] 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) = j'K2(x)dx V(g) = Jg2(x)dx fik{K) = J xkK(x)dx konvoluce funkcí / a g, (/ * g)(x) = J f (ť) g (x — t) dt integrál z jádra, W(x) = J* K (i) 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 znaci, ze Podobně symbol znaci, ze 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 f(x) = 0(g(x)) pro x ->• a lim sup—< oo. f(x) = o(g(x)) pro x ->• a x^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)) = 0(g(x)) °{g{x)) + °{g{x)) = °{g{x)) 0{g{x)) + o{g{x)) = 0{g{x)) 0{g{x)) ■ 0{h{x)) = 0(g(x) ■ h{x)) °(g(x)) ■ °(h(x)) = °(g(x) ■ h(x)) 0(g(x)) ■ o(h(x)) = o(g(x) ■ h{x)) o{g{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. Nechť v, k jsou nezáporná celá čísla, 0 < v < k, nechť 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,y e [-1,1], L > 0, 2. nosič(ŕf) = [-1,1], tj. K = 0 vně intervalu [-1,1], 3. K splňuje momentové podmínky xJK(x) dx = < 0<3>)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: S 02 S13 S 24 KoPt,iAx) = f1^2 - iKi-i.i]^) K, opt,2,4 (x) 16 (x2-í)(5x2-í)I{_hl](x) Jádra -KT0pt,i,fc a Kopt^,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 Sjy+i^+i je jádro s minimálním rozptylem a Kopt,v,k <= S„fc je optimální jádro. Pak platí KPt,„,k(x) = K(x), a; e [-1,1] Jako příklad lze uvést jádro Koptfl^{x) = |(1 - x2) e S02 a K'opt02(x) = Klt3(x) = —|x G Sis. 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 í3k{K) = J xkK{x)áx V(K)= J K2(x)dx T(K) = (\pk(K)\2"+1V(K)V°-"^ 5vk = (J^j 2fc" 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 Nechť K e 5„fc, veN, pro 7 > 0 položíme Jádra K a K1 nazýváme ekvivalentní. Nechf jsou dány funkce /ag, pro které platí J2 (x) dx < 00 a j g2 (x) dx < 00. Konvoluci f * g definujeme vztahem (f*9)(x)= í f(t)g(x-t)dt. Cvičeni 1. Dokažte, že funkcionál T(K) je invariantní vzhledem k transformaci K7(-), tj. T(K) = T(K1). 2. Vypočtěte kanonické faktory 5yk a hodnoty funkcionálu T(K) pro jádra z tabulky 1.2. 3. Ukažtem že pro konvoluci f * g platí • f*g = g*f, 9 • / *{g*h) = {f *g)*h, • f*(g + h)=f*g + f*h. 4. Nechť K e S02■ Ukažte, že platí následující vztah 1 j = O, (K*K)(x)xJdx = { O j = 1, 2/32(íO j = 2. 5. Jak je třeba zvolit konstantu A, aby funkce K{x) A cos x e [—1,1], 0 jinak, byla jádrem třídy S02? Dále dopočítejte hodnoty f3(K) a V(K) pro toto jádro. 10 Kapitola 2 Jádrové odhady regresní funkce 1 Motivace Uvažujme datový soubor, který obsahuje měření úrovně hladiny Hurónského jezera. Hurónské jezero je druhé největší jezero v systému pěti velkých jezer v Severní Americe. Jezerem prochází státní hranice mezi Kanadou a USA.1 Měření byla prováděna ročně, v letech 1875 až 1972, a výsledky měření jsou zobrazeny na obrázku 2.1. Naším cílem je najít funkci popisující úroveň hladiny v uvedených letech. Vidíme, Obrázek 2.1: Úroveň hladiny Hurónského jezera že pouhý pohled na tento dvourozměrný bodový diagram obvykle nestačí k tomu, abychom určili tento funkční vztah. Statistická úloha, kterou se budeme zabývat, spočívá v proložení vhodné křivky těmito body tak, aby byly odfiltrovány náhodné výkyvy a bylo možné lépe poznat strukturu dat. Tuto křivku nazýváme regresní křivkou. Formalizujme nyní tuto úlohu: Uvažujme standardní regresní model Yj = m(xi) + Bi, i = l,...,n, (2.1) kde m je neznámá regresní funkce, x i, í = 1,..., n, jsou body plánu, Yi, i = 1,..., n, jsou hodnoty závisle proměnné a. e i, i = 1,..., n, jsou chyby měření, o nichž se předpokládá, že jsou nezávislými, identicky rozdělenými náhodnými veličinami splňujícími podmínky Eeí = 0, var£i=er2, i = l,...,n. (2-2) Poznámka 2.1. Jsou-li body plánu uspořádaná nenáhodná čísla, mluvíme o regresním modelu s pevným plánem. V případě, že body plánu X\,..., Xn jsou náhodné veličiny se stejnou hustotou 1 „Great Lakes from space" od SeaWiFS Project, NASA/Goddard Space Flight Center, and ORBIMAGE. - http: // visibleearth. nasa. gov/view_rec. php?id=7 93. Licencováno pod Public domain via Wikimedia Commons. 11 /, jedná se o regresní model s náhodným plánem (podrobněji např. [14]). Budeme se dále zabývat modelem s pevným plánem. Bez újmy na obecnosti budeme v dalším předpokládat, že pro body Xi,i = 1,..., n, platí 0 < x\ < X2 < • • • < xn < 1. Cílem regresní analýzy je nalézt vhodnou aproximaci fh neznámé funkce m. Tento proces odhadu regresní funkce se obvykle nazývá vyhlazování. K tomuto úkolu lze přistoupit dvěma způsoby - parametricky a neparametrický Příkladem parametrického odhadu regresní funkce je regresní přímka vyjadřující lineární závislost. Naopak u neparametrického přístupu nepředpokládáme, že funkce má nějaký předepsaný tvar, pouze předpokládáme jistou hladkost odhadované funkce (tj. dostatečný počet spojitých derivací). V první polovině dvacátého století byla věnována pozornost zejména parametrickým metodám. V posledních letech však zaznamenaly značný rozvoj neparametrické metody. Tento vývoj souvisí s rostoucími požadavky na zpracování dat, af už jde o rozsah souborů, rozmanitost těchto dat apod. Cistě parametrický přístup nevyhovuje vždy potřebám flexibility a nebývalý rozmach výpočetní techniky vytvořil dobré předpoklady pro rozvoj neparametrických metod. I přes tento vývoj si oba způsoby zachovávají své výhody a nijak si nekonkurují. Někdy je vhodné užít neparametrické metody a pak na výsledný odhad použít parametrickou metodu. Příklad 2.1. Obr. 2.2 ilustruje na simulovaných datech nevhodnost aplikace parametrického přístupu. V tomto případě byla data generována podle vztahu sin Attxj Yi = _-__\- £i (1 + cos Qfi-nxi)2 ' kde body Xi = i/100, i = 1,..., 100, a chyby Eí, í = 1,..., 100, mají normální rozdělení n(0; 0,25). (Data jsou v tabulce 7.1.) Obrázek 2.2: Simulovaná data (x) s regresní přímkou (červená, plná) a původní funkcí (modrá, čárkovaná) Předpokládejme, že hledaná křivka je přímka a známou metodou nejmenších čtverců určíme rovnici této přímky. Obr. 2.2 znázorňuje přesnou funkci, generovaná data a výslednou přímku. Je zřejmé, že náš předpoklad, že hledaná funkce je přímka, není správný. 2 Základní typy neparametrických odhadů Pokud jde o historii neparametrických metod, připomeňme, že v r. 1857 saský ekonom Engel analyzoval data týkající se nákladů na domácnost a pro vyjádření závislosti použil schodovitou (tj. po částech konstantní funkci), kterou dnes nazýváme regresogram. Regresogram užívá stejné základní myšlenky jako histogram pro odhad hustoty. Tato myšlenka spočívá v rozdělení 12 množiny hodnot proměnné X na intervaly B j, j = 1,..., J, a za odhad v bodě x G B j se vezme průměr hodnot Y na tomto sub intervalu, tj. E YiI[Bj]{xi) ™(x) = ^-> i=l kde -f [Bj ] je indikátorová funkce sub intervalu B j. Výsledek aplikace regresogramu na simulovaná data z příkladu 2.1 je znázorněn na obr. 2.3. Vidíme, že tento odhad „vhodně" vystihuje tvar funkce, ale výsledný odhad je příliš hrubý. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obrázek 2.3: Regresogram (červená, plná) pro simulovaná data z příkladu 2.1 s původní funkcí (modrá, čárkovaná) Přirozeným zobecněním regresogramu je metoda klouzavých průměrů. Tato metoda používá lokálních průměrů hodnot Y, ale odhad v bodě x je založen na centrovaném okolí bodu x: [x — h,x + h], h > 0, přesněji m(x, h) E Yil[x-h,x+h] (Xi) i=l_ n E I[x-h,x+h] (xi) (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)Yl, (2.4) i=i 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 Sok, k je sudé číslo, položme Kh(-) = j^K^). Mezi nejznámější typy jádrových odhadů regresní funkce patří ([8]): 1. Nadarayovy-Watsonovy odhady (1964) Y, Kh(x - xt)Yt mNW(x,h) = ^—i- Y, Kh(x - Xi) i=l 2. Priestleyovy-Chaovy odhady (1972) I " ;(x,h) = - Kh(x - Xí)Yí, II < J mPC( n ■ i=l 3. lokálně lineární odhady (Stone 1977, Cleveland 1979) ^ ^ ^ _ í_ {h{x,h) - ši(x,h)(x - Xj)}Kh(x - xt)Yt i=l kde n ^—^ Š2(x,h)šo(x,h) — ši(x,h) 1 " šr(x) = - y^(x - xl)rKh(x - Xi), r = 0,1, 2 i=l 4. Gasserovy-Můllerovy odhady (1979) n s* r(x,h)=J2^ 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 1(x,h) = J2Wl(A)(x,h)Yl, 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 řriMw- Popíšeme nyní jeho konstrukci a budeme ilustrovat vliv vyhlazovacího parametru na kvalitu odhadu. Pro daný bod xq, h < xo < 1 — h, jsou váhy Nadarayova-Watsonova odhadu dány vztahem Wi(x0,h)= Kh(£°-Xi) Y.W1(xí)lh) = l. Ei=i Kh(x0 - xj) ^ Obrázek 2.5 ilustruje konstrukci odhadu v bodě xq, který je založen na pěti pozorováních (xi ,Yi), ..., (x5, y5) (černé křížky). Parabola reprezentuje Epanečnikovo jádro a kroužky znázorňují hodnoty vah Wi(x0,h) = Kh(x0 — x i) / Ei=i Kh(xo — x i) pro i = 1,..., 5. Výsledný odhad regresní funkce fh v bodě xq 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 (xo, h) pro odhad s obdélníkovým jádrem. (Odpověď viz kapitola 6.) Jádrový odhad není definován pro Eľ=i Kh(x — x i) = 0. Jestliže nastane případ „0/0", pak klademe rriNw{x, h) = 0- Omezíme se nyní na odhady funkce m v bodech plánu Xi,i = 1,..., n. Pro malé hodnoty h je výraz Xz f^3 > 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(xl,h) ->• = 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 Xz ^3 « 0, tedy pro všechny body plánu je hodnota jádrové funkce K ^a"' fea'J ^ « ^(0) a dostaneme tak průměr dat n n E*W- ^(O)E^ x n mNW(Xi,h) -i- ^4-= = ~J2Yr E K(0) [U) 3=1 i=i Tedy velká šířka okna (h —>• oo) 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í, nebol i asymptoticky optimální odhady obsahují poměrně značné ,pxnož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(a:,/i)-Mi-¥y ^ , m(x, h) + Ml_f ý ^ , kde 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) = J2 Wi(x, h) (Yi - m(x, h)f. i=i 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(i, h) = Efh2(x, h) — 2m(x)Eíh(x, h) + m2(x) = (Em(x, h) - m(x)f + Em2(x, h) - (Em(x, h)f, (2.5) ^ v ^ ^ v ^ bias2 var což znamená, že střední kvadratická chyba může být vyjádřena jako součet rozptylu odhadu varm(x, h) a čtverce vychýlení bias2 m(x, h). Tento rozklad rozptyl-vychýlení usnadňuje analýzu vlastností odhadu. Všechny uvedené jádrové odhady regresní funkce m^wr ffipci ™-LLi ™gm 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: fh(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 502/ (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=l,...,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 i " i " ŕ i /x-t\ Em(x,h) = -Y^Kh(x-xt)EYt =-Y^Kh(x-xt)m(xt)= / -K —— )m{t) dt + O^1). ti ti in n v n / i=i i=i j[> (2.6) Nechf u = ^j^, odtud dí = —hdu, a tedy s využitím Taylorova rozvoje f-x/h /x / a K(u)m(x — hu) du + 0(n_1) -(l-x)/h (l-x)/h x/h r u2h2 K(u) m{x) - uhm'(x) H--— m"(x) + o(h2) dM + 0(n_1). (2.7) -(l-x)/h L 2 i 2Symbolika 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— —>• — oo pro h —>• 0. Odtud, s využitím faktu, že nosičem funkce je interval [—1,1], plyne Efh(x, h) = m(x) J K(u) du — hm'(x) J u K (u) du+—m"(x) J u K(u)du+o(h )+0(n ). = 1 =0 =fo{K) Celkem dostaneme h2 biasm(x, h) = —h(K)m"(x) + o(h2) + C^n"1) « ^(32(K)m"(x). Podobně pro rozptyl platí var ŕh{x, h) = E(m(x, h) — Efh(x, /i))2 /l ™ 1 " = E[ Kh(x - xl)Yl--V]Kh(x - xl)m(xl) , n í—' n v i=l i=l 1 / " = —Ey ^Kh(x - xl)(Yl - m(xi)) n z vlastností (2.2) plyne EK^x^K^x^EíEj = 0 pro i ^ j, tedy 1 " nz z—' i=l i=l Zde jsme opět použili přibližného výpočtu integrálu. Opět s využitím substituce u = ^j^- a vztahu 0(n_1) = o((n/i)_1) můžeme pro n —>• oo psát cr2 f1 cr2 f1 cr2 varm(x,h) = — K2(u) du + oíính)-1) = — K2(u) du + oíính)-1) m —V(K). nh 7_i w/i 7_i Tímto jsme dokázali následující větu o tvaru střední kvadratické chyby. Věta 2.1. Nectí jsou splněny predpoklady (i) — (iii), pak střední kvadratická chyba nabývá tvaru a2 1 ? MSEm(i,/i) = — V (K) + h4 p2 (K)-(m" (x)) +o(h4 + (nhy1). (2.8) 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) = í MSEfh(x,h)dx = AMISEfh(-,h)+o(h4 + (nh)-1). Jo AMISE je tvaru u2 h4 AMISEm(-, h) = -^V(K) + —^2(K)V(m"), (2.9) AIV(h) AISB(h) 18 Obrázek 2.9: AMISE (růžová, plná) jako součet rozptylu AIV (č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 ATV(h) a AISB(/i) a také výsledné chyby AMISE fh(-, h) Je vidět, že rozptyl ATV(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) _ a2V(K) dh nh2 položíme ji rovnu nule a vyjádříme h h^i(K)V(m"), o2V(K) 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,o,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,o,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) = £„W (,,^U)"/5 + í (^ff^^WK) = (a2)ilb(V(K)ýlbn-ilb(p2(K))1lb(V(m"))1lb + ^(cr2)4/5(V(K))4/5n-4/5(p2(K))1/5(V(m"))1/5 = \(a2V(K))A/*(p2(K)V(rn"))1/5n-4í\ (2.11) tj. AMISEm(; ^,0,2) = Ofc-VS). 19 Poznámka 2.4. Jestliže jádro K náleží do třídy Sok, pak AMISE je tvaru AMISE m(-,/i) = —rV(K) nh (k\)2 h2kpl{K)V{m^) (2.12) a pro optimální vyhlazovací parametr h platí ' 2fc+l (j2V{K){k\)2 (2.13) opt'°'k 2knP2{K)V{m(k)y kde V(m(fe)) = (m^(x))2 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,fc pZařř AlVt/iopt.o,*) = 2A;AISB(/loptj0jfe). Důkaz. Viz cvičení. □ Lze ukázat, že pro jádra K e Sok je AMISE m(-, h) = O(n~ 2fc+!). To znamená, že s rostoucím 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, 1], 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 irx na intervalu x e [0,1] s chybami Ei ~ A(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 ttx m'(x) = 2 sin ttx cos ttxtt = 7rsin27ra; m"(x) = tt cos 2irx2ir = 2ir2 cos 2irx V(m") = í [m"(x)]2dx = í [27r2cos27rx]2dx Jo Jo = 4tt4 / cos2 2ttx dx = 4"7r Jo Jo = 2k\ Výpočet hoptfl,2 pro • Epanečnikovo jádro: V(K) = 3/5, ^(K) = 1/5, ít23/5 1 + cos AlTX ■ dx 11„ Iber2 opt'0'2 n(l/5)227T4 27T4n' obdélníkové jádro: V(K) = 1/2, ß2(K) = 1/3, (T2l/2 h opt,0,2 9<72 n(l/3)227T4 4ir4n 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.4 0.5 0.B 0.2 0.3 (a) Epanečnikovo jádro, hopt,o,2 = 0,1573 (b) Obdélníkové jádro, hopt,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 Sok- 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) = ^2 We(xi,h)Ye, i = í,...,n. i=i Funkce křížového ověřování je definována takto 1 ™ CV(h) = -J2(m-l(xl,h)-Yt)2 (2.14) i=i a odhadem optimální hodnoty vyhlazovacího parametru je bod, v němž nastává minimum této funkce, tj. hopt,o,k = hCv = arg min CV(h). h£Hn Hledáme tedy minimum na intervalu Hn = [afcn_1/(2fe+1), &fcn_1/(2fe+1)], jehož tvar plyne ze vztahu (2.13), přičemž a^, bk jsou konstanty (0 < < 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(/i) platí 1 ™ ECV(h) = -E^(m{xi,h) -m(xt))2+a2. 71 i=l AMSE Důkaz. Funkci křížového ověřování lze rozepsat 1 ™ 2 CV(h) = — ^(m_j(aľj, h) - m(xi) - {Y% - m(xi))) ^ n \ n I n = - y^(m_j(aľj, h) -m(xi)) - 2- y^(m^t(xu h) -m(xl))el + - e\ i=i i=i i=i = - ^(m_j(xj,/i) - m(xi))--^eJ^Weix^^Ye - m(Xi)) + - ^e: n *—' n *—' \*—' / n i=l i=l t=l i=l Střední hodnota E CV(h) je rovna součtu tří veličin. Předpokládejme, že m_i(xi, h) « m(xi, h), pak první ze sčítanců je roven přímo AMSEm(-, h): 1 ™ -E^2(m-i(xi,h) -m(xt))2 = AMSEm(-, h). n i=l 22 Dále víme, že Yi = m(xi) + e^, a tedy pro druhý sčítanec platí: i=i e=i 1^1 =--E^2ei[^2We(xi,h)(m(xe) + et) - m(xl) i=i e=i 1^1 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 E CV(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. 0I-1-1-1-1-1-1- 0 0.1 0.2 0.3 0.4 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) = ^V(K) + J^h2k (32(K)V(m^) a hopt,o,k tvaru h opt'°'k 2knV{m(k)) ^ Ze vztahu pro /iopt,o,fe vypočteme V(nrSk^) a dosadíme do vztahu pro AMISE. Dále použijeme vztahy Sok = (^fj)3^ T(K) = ((3k(K)Vk(K))^ T(K) = 52kp2{K) a dostaneme vyjádření AMISEm(.,/lopt,o,fe) = °y* + 1)S°kT(K), ve kterém jsou parametry K, h a k separovány, což umožňuje vybrat tyto parametry simultánně. Právě tento vztah je základem automatické procedury. Položme 2nkhoptfltk a množinu vhodných řádů k označme /(*:o) = {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 G Sok, které je dáno tabulkou 8.1 a vypočtěte kanonický faktor Sok- 2. Pro k G I(ko) a Koptto,k G Sofc najděte optimální vyhlazovací parametr /iopt,o,fe- 3. Pro G 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ů 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 hoptt0,k- k Kopt,0,k Sok hopt,o,k L(k) 2 -i(-2- -1) 1,7188 0,1158 0,0648 4 15 IJ2 32 vx l)(7x2-3) 2,0165 0,2446 0,0575 6 105/2 256 vx - l)(33x4- 30x2 + 5) 2,0834 0,3575 0,0574 8 315 (2 4096 ^ - 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. ........................% X Jť**^ Nil, ''*f • * /' jtlf....... ...... yrý//"" " 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 (Kopt^2), h = 0,1525 (^,0,12). Výsledná křivka rh(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. v". 1B70 1BB0 1B90 1900 1910 1920 1930 1940 1950 1960 1970 19B0 (a) Data (b) hcv = 0,0204, Xopt,0,2 (c) h = 0,1525, Koptfi,i2 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 (tfopf.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 27 Shrnutí Odhad regresní funkce m(x) v bodě x je tvaru n m(x, h) = ^ Wi(x, h)Yi. i=i 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), var fh{x, h) « — V(K). Asymptotická střední kvadratická chyba jádrového odhadu regresní funkce s jádrem řádu 2 AMISEm(-,/i) = ^V(K) + ^h4p2(K)V(m"). Optimální vyhlazovací parametr vzhledem k AMISE pro k = 2 je tvaru iL A 'AMISE - nß2{K)V{mll) s řádem konvergence rT1^. Metoda křížového ověřování pro odhad optimálního vyhlazovacího parametru 1 ™ CV(h) = -^(fh^ix^h) -Yi)2 hCv = 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 Výpočet integrálu, ŕ j = ^, i = 1,..., n, /•i /-ti+i / G(í) dt=J2 G(*) dí s využitím Taylorova rozvoje funkce G(t) n-l rti+1 = J2 / (G(íi+i) + (í - íi+i)G'(íi+i) + o(l)) dí i i /-ti+i = J2 G{U+i){U+1 -U) + J2 (* - *í+i)G'(*í+i) dí + o(l) = -H G {ti) + £ [U+\ U) G'(tl+1) + o(l) 1 n -. n — 1 i=l i=0 Za předpokladu |G'(íí+1)| < M, pak platí ŕ 1 " / G(t)dí=-VG(íí) + 0(n-1). Jo Cvičeni 1. Odvoďte tvar Priestleyova-Chaova odhadu regresní funkce v bodě x pro obdélníkové jádro. 2. Pro odhad ttimw dokažte, že „množství" vyhlazení pomocí jádra K s vyhlazovacím parametrem h je stejné, jako „množství" vyhlazení jádrem Ks s parametrem h* = h/S, tj. rh(x, h, K) = ih(x, h*, Ks). 3. Uvažujte funkci m(x) = sin27rx a regresní model Yi = m(x) + Ei, i = í,dots, 100, kde Ei ~ 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^o,,) = n-^/^Í^ViK))2'^2^ 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 i i i i i 0-1- -H- + HIM lili MIHNIlllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllll I III + + + +- 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, • Jf(x)dx=í. 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 f(x) = 0,3 ,„ , e~fe +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 3 3.5 (a) 7 intervalů o šířce h = 0,55 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 (b) 13 intervalů o šířce h = 0,29 Obrázek 3.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 R je definovaný vztahem [14]. 1 " nh ^-^ V h i=l x - Xi I " -y2Kh(x-Xi), II —' (3.1) i=l K e Sok 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 v-té derivace hustoty. Budeme předpokládat, iei) • 0 a nh —>• oo pro n —>• oo. PflA: p/flrř MISE/(■,/») = ^ + pj^W*') +o(/í2fe + (n/í)-1), kdeV(fW) = J(fik)(x))2 dx. Důkaz. Nejprve vypočteme střední hodnotu 1 /x - y Ef(x,h) = (Kh*f)(x) = j Kh(x-y)f(y)dy = j -K^-^)f(y)dy = / K(z)f(x - hz) dz dále použijeme Taylorův rozvoj: f(x — hz) = f(x) — f'{x)hz + • • • + fcV f^(x)hkzk + o{hk) = / m im- n*v«+- + *&ŕw+°«><->] * = /(x) / lf(z)dz-/'(x)/i / zif(z)dz + --- + i1^/(fe)(a;)/lfe /* zfeif(z) dz +o(/ife) =1 =0 =Pk(K) = f(x) + {~l)kf^){x)hkh{K) + o{hk). Tedy vychýlení odhadu je tvaru bmSf(;h)=Ef(x,h)-f(x) = ( ^^{x) hk (3k(K) + o(hk), =o(l) 33 a tedy Ef(x,h)=f(x) + o(í). Nyní dokážeme vztah pro rozptyl. Víme, že var/(a;, h) = ^((K2 * /)(*) - (Kh * /)2(x)) (/(:r)+o(l))2 a dále počítáme ^ I if2(z) /(* - hz) dz - i(/(a;) + o(l))2 =/(2;)+o(l) ^ | ^2(z)(/(x) + 0(l))dz-i(/(x) + 0(l))2 Tedy o((nh)-i) ^ /* Ä-2(z)dz + o((ra/i)_1). MSE/(a;,/i) = ^ / K2(z) dz + o^n/í)-1) + (L-^^hk(3k(K) + o(hk) nh J \ k\ = 4^/ K2(z)dz + o((nh)-í)+f(k\x)^t(K) , 0(-l)fe/(fe)W,fefl , „,,.2^ jfc! -hk/3k(K)(hk)+o(h2k) =o(h2k) a pak využijeme faktu, že J f (x) dx = 1 MISE/(-,/i) = J MSVf(x,h)dx □ Důsledek. Nectí h —>• 0, n/i —>• oo pro n —>• oo, paA: / ;'e konzistentním odhadem f, tj. E f —>• / « var / —>• 0. Stejně jako u odhadu regresní funkce má význam asymptotická integrální střední kvadratická chyba AMISE/(-, h): MISE/O, h) = AMISE/(-, h) + o(h2k + (nh)-1), kde AMISE je tvaru AMISE/(■,/») = ^ + 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(ft) = -±^h2kPl{K)VUW), tedy AMISE f(-,h) = AIV(h) + AISB(/i). Užitím vztahů T(K) = (V(K)kpk(K))2/(2k+1) a S^+1 = pro K e S0k lze AMISE zapsat ve tvaru AMISE f{.,h) = T(K)(^ + ^^l). (3.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(f^)' 1 ' tj. Kpt^k = 0(n-W°+V). Vypočtěme hodnotu AMISE při dosazení optimálního parametru hoptfl,k- AMISE f{,hMk) = r(g)v(/W) WD^^+i) (2fc(^2+;(2fc+1), (3.6) tj. AMISE/(-, Vt.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,o,fe) = 2A: AISB(/iopt,o,fe)- (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 Vt,0,2 - nV{f//) ■ Počítejme derivace AMISE /(■, h) dané rovnicí (3.3) pro k = 2 d2AMISE/'(-,/l) 2V(K) 3h2f32(K)V(f") dh2 nh? 6h(3l(K)V(f"). d3AMISE/(-,/i) _ -6V(K) , a 2 dh3 nh4 Řešením rovnice d3AMISE /(■, h)/ dh3 = 0 je 5 = V(K) = 5l2 = n(32(K)V(f") nV(f") op tj- hoptfl,2 také realizuje minimum d2AMISE /(■, h)/ dh2. Lze ukázat, že d2AMISE(/'(-,/l)) dh2 = 0{n 2fc+i) h—h0pt^Q^k 35 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,o,fe nevede k velkému růstu AMISE f(-,h). Na obrázku 3.4 jsou zobrazeny body minima funkce AMISE f(-,h) pro hustotu normálního rozdělení n(0; 1) se sto prvky. Vztah pro optimální hodnotu vyhlazovacího parametru poskytuje informaci, že asymptoticky je h = 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 Hn = [akn-W\bkn-W% kde aj;, 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 e [—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 V(f")= yV'(x)]2dx = 35 Výpočet hopt,o,2 pro Epanečnikovo jádro: v(k) = 3/5, ^(k) = 1/5, tedy 5q2 = jt^) = ^ 5 _ 15(2!)2 _ 3_ Vt,o,2 2-2-n-35 7n Tedy pro soubor 50 hodnot bude /iopt,o,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,o,2 = 0,3860 Předpokládejme nyní, že platí 0 < v < k — 2, v & k mají stejnou paritu, linin^oo h = 0, Ymin^conh2v+1 = oo, / e Ck° (k < k0) a V(f^) = J(f^(x)fdx < oo. Pak lze ukázat, že asymptotická střední kvadratická chyba AMISE f^u\-, h) je tvaru AMISE/»(■,/») = + j^h2^(32(K^)V(f^). 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^ + i)(fcQ2 opt'v'k 2n(k - v)V(fM) ' kde ^=V^ (3.8) Tento vzorec umožňuje výpočet optimálního vyhlazovacího parametru pro pomocí /iopt,o,fc a hopt,i,k- Předpokládejme nejdříve, že v a k jsou sudá čísla. Pak h opt,v,k h opt, 0, k Pro ľ a k lichá platí h opt, v, k h opt ,l,k (2, + í)k\1/^k+1)5^ Sok (2^ + l)(fc-l)\1/(2fc+1)^fc 3(k -v) ) Slk ■ Speciálně pro v = 2, k = 4 dostáváme velmi užitečný vztah (3.9) (3.10) h, opt,2,4 004 opt, 0,4; (3.11) pricemz K0ptflA{x) = g(x2-l)(7x2-3), í>04 í>24 2,0165, 1,3925. 4 Volba jádra Volba jádra není z asymptotického hlediska podstatná, jak je zřejmé z faktu (3.6). Je vhodné zvolit optimální jádro, které minimalizuje fukcionál T(K), nebof tato jádra jsou spojitá na M a hladkost jádra také „zdědí" odhadovaná hustota. Poznámka 3.3. Při výpočtech se používá nejčastěji Epanečnikovo jádro. 37 5 Volba vyhlazovacího parametru Volba vyhlazovacího parametru pro jádrový odhad hustoty je, stejně jako u regrese, zásadním problémem. I když tomuto problému byla a je věnována značná pozornost, doposud neexistuje univerzální přístup k řešení tohoto problému. Nejjednodušší metoda je „okometrická". Je účelné „nakreslit" několik křivek s různými vyhlazovacími parametry dříve než užijeme nějakou automatickou proceduru. Je třeba zdůraznit, že z hlediska analýzy všechny volby vyhlazovacího parametru vedou k užitečnému odhadu hustoty. Velká šířka okna charakterizuje globální strukturu hustoty a naopak malá šířka odhaluje lokální strukturu, která může nebo nemusí být přítomná v přesné hustotě. Tuto myšlenku ilustruje obrázek 3.6, na němž jsou zobrazeny odhady pro simulovaná data (f(x) ~ ^(0; 1), n = 100) s hodnotami vyhlazovacího parametru z intervalu [0,05,1]. Jednotlivé odhady hustoty příslušející těmto hodnotám jsou znázorněny tenkými čarami. Silná křivka znázorňuje odhad s optimální hodnotou h = 0,4217. Třída těchto odhadů ukazuje široký rozsah vyhlazení od podhlazení (modré křivky) až k přehlazení (červené křivky). 0.5 0.4 0.3 0.2 0.1 0 Obrázek 3.6: Volba vyhlazovacího parametru V dalších úvahách bude užitečná následující definice: Definice 3.1. Nechť h je odhad hoptto,k- Řekneme, že h konverguje k hoptto,k s relativní rychlostí ■nTa, jestliže h-hMk =0{n-ay hopt,o,k 5.1 Metoda referenční hustoty Nejčastěji se pro odhad neznámé veličiny V(f^) (viz rovnice (3.3)) používá parametrické třídy hustot. Jednou z možností je použít standardní normální hustotu / s rozptylem a1, tj. předpokládáme, že 1 *2 m — V tomto případě je odhad optimálního vyhlazovacího parametru tvaru pro K e Sok 38 Je třeba ještě odhadnout směrodatnou odchylku a. To lze dvěma způsoby: 1 n 1/2 1 n i=l i=l (3.13) (3.14) kde <3> 1 je standardní normální kvantilová funkce a číslo X[3„/4], respektive X[„/4], je horní, respektive dolní výběrový kvartil. Je vhodné volit mm{$sj}, piqr}. Poznámka 3.4. Pokud za jádro K zvolíme Gaussovo jádro (k = 2) 1 K{x) 2tt e 2 pak dostaneme jednoduchý vztah [11,12] _ / 4 \i/5 ~ \Šň) /íREF = ( — ) f. (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 h REF S02vn 1/5. Dále odhadneme směrodatnou odchylku: Ôsd = 1,0325, piqr = 1,3344, tedy a = 1,0325. Po dosazení počtu prvků n = 100 a parametru t>02 = 1,7 1 88 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,o,fc 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)^1 |x| 0, j (rg^ (rx))2 dx = r2fe+1 J(V'0-' (x))2 dx pro každou hustotu, pro kterou integrál existuje. 3. Jestliže hustota g má rozptyl a2, pak hustota ^-g(^-aľ) má rozptyl a2. Jestliže / je neznámá hustota s rozptylem a1 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 hopt,o,k < ^ofc Hodnotu a lze odhadnout pomocí dříve uvedených vztahů a cr^ = 2fc1|_5. Hodnotu V(g^) lze vypočítat pomocí speciálních ortogonálních polynomů [5]: V(o^) ľ Uk\x)Ýdx _(2fc + 3)!(2fc + 2)!_ (3 17) Použijeme-li poslední vyjádření (3.17), dostaneme horní hranici pro vyhlazovací parametry hopt,o,k < hMS = čn-VW+Vh, pricemz 22k+2V(K)(2k + í)(2k + 5)(k + í)2(k\)2\^ 6fe-V2fcT5, /32W(2A. + 3)!(2A. + 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 /ims Pro odhad hustoty s jádrem řádu k = 2, tj. K e So2-Podle věty 3.3 je hustota g2 tvaru 35 92(x) = -(l-x2f, a; e [-1,1], a dále z vlastností funkce g2 a ze vztahu (3.17) plyne x2g2(x)dx= — a / [g"(x)Y dx = 35. J-i ° pro Äoptj0,2 /32(^) 35 j-i v J-i Pak pro K,^„ u -i/5fV(K) 243\ 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 ^ _1/5/ 3/5 243 \ VŠ 0409, 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 h-^s může sloužit jako horní hranice pro množinu vyhlazovacích parametrů volených podle jiné metody, např. metody křížového ověřování. Tedy Hn = [hg, /ims]/ kde hn je nejmenší vzdálenost mezi po sobě jdoucími body Xi, i = 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/(-,/l) = S j (f(x,h)-f(x)Ydx = E f2(x,h)dx -2E / f(x,h)f(x)dx+ I f2(x)dx. mSEf(;h)-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)= f2(x,h)dx--J2f-*(X*,h), J 71 i=l kde f-i (Xi, h) je odhad v bodě Xi bez použití tohoto bodu. Věta 3.4. Platí EGV(h) =MISE/(-,/i) - J f(x)dx, tj. CV(/i) je nevychýleným odhadem (3.18) E[ / f2{x,h)dx-2 / f(x,h)f(x)dx 41 Důkaz. Spočítejme střední hodnotu funkce CV(h), tj. ľ 2 " ECV(h) = E / f2(x,h)dx--Ej2f-*(X*,h). J 71 2=1 Dále upravme druhý člen tohoto vyjádření: 1 n 1 n 1 n e- J2 f-i(Xi,h) = e-J2 —7 E K^x* - x^ n ^-^ n ^-^ n — 1 ^-^ 2=1 2 — 1 j —1 E , * É É Kh(Xl - X j) = EKhÍX! - X2) n(n — í) ^-^ ^-^ y > 2=1 j=l Kh(x-y)f(y)f(x)dxdy = E / f(x,h)f(x)dx E f {x,h) Odtud EGV(h) = E J f2(x,h)dx-2E J f(x,h)f(x)dx. Odhad hoptok je dán vztahem □ hcv = arg min CV(/i). Odtud plyne, že CV(/i) + J /2(a;) dx je pro každé h nevychýleným odhadem MISE f(-,h). Protože J f2(x)dx nezávisí na h, minimalizace ECV(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 (f")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 -4-3-2-1 0 1 2 3 4 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 hcy = 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): AlV(hopt^k) = 2kAISB(hopt^k). Přepíšeme tuto rovnici 2kh2k^-V{f^) = 0. (3.19) nh (k\)z 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 var/(x, h) = \ í K2(y)f(x - hy) dy nh bias/(x, h) = (Kh * f)(x, h) - f(x, h) = I f(x - hy, h)K{y) dy - f(x, h). Odtud plyne AW(h) = (3.20) nh 43 x - X, Xi — X j AISB/(-, h) je vychýleným odhadem AISB, a proto budeme uvažovat ....... 1 " AISB(ft) = — J2 ^ÁX.-X,). Místo rovnice (3.19) řešíme rovnici Rovnici lze také zapsat ve tvaru: 2 I dx. AISB(ft) = J (J f(x-hy,h)K(y)dy-f(x,h^j dx \ i—l ' i—l Výraz lze upravit pomocí konvolucí a dostaneme 1 " = — V (K*K*K*K-2K*K*K + K*K) n2h ^ n2h ^ \ h J Funkce A(z) = (K * K * K * K — 2K * K * K + K * K)(z) má tyto vlastnosti J zjA(z) dz = 0, j = 0,1,..., 2k - 1, z2kA(z)dz= (3.21) (3.22) _ nV(K) _ 2kY,tJ=iK{xl-x1)-[)- [Ó-ZÓ) 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 hopt0j„. Tato skutečnost je dokázána v následující větě [3]. Věta 3.5. Nechí 44 Pak EP(h) = P(h) +o(h2k+1), Rk2 varP(h) = -— V(A)V(f) +o(n-2h-1). nzh Poznámka 3.7. Rychlost konvergence pro metodu křížového ověřování: hcv — hoptflp h. opt,0,2 Rychlost konvergence pro iterační metodu hn - hoptfla = Q^n-i/w^ nopt,0,2 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 čárkovaná) 0,5314, odhad (červená, plná), původní funkce (modrá, 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/(-,/loptj0jfe) = T(^) Ze vztahu pro /iopt,o,fe vypočteme V(f^) nh, Sok , ^,0,fcn/(fc))\ S2k(kl)2 )■ opt,0,k a tuto hodnotu dosadíme do předchozího vztahu: AMISEf(;hoptt0tk)=T(K) (2k + í)Sok 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 G I(k0) najděte optimální jádro Koptok 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 /iopt,o,fc- 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. />-w) = -r—ž^(F^)- nľlopt,0,k i=l ^nopt,0,k' 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 15 IJ2 32 vx l)(7x2-3) 2,0165 1,0734 0,0131 6 105 í- 2 256 vx - l)(33x4- 30x2 + 5) 2,0834 1,6460 0,0125 8 315 / 2 4096 vx - 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 S^, např. kvar-tické jádro: K(x) = j|(l — x2)2I[_ltl](x), nebo jádro triweight: K(x) = ||(1 — x2)3I[_ltl](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, hREF = 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 ftopt,o,2 = 0,5122 (b) Odhad s hREF = 0,9639 (e) Odhad s hjT = 0,5314 (f) Odhad s h = 1,6460 a K e S06 Obrázek 3.15: Srovnání odhadů pro data z příkladu 3.1 48 (e) Odhad podle procedury s h = 140,46 a X 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, /iMS = 31,62, hcv = 30,55, ftIT = 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 Kq$. 49 (a) Odhad s hREF = 5,7856 (b) Odhad s hMs = 6,2480 (e) Odhad podle procedury s h = 31,73 a K e 5o8 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 1 n y 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/(, h) = ?¥p- + J^h2k(32(K)V(f^). AIVW " IÍSB(I) ' Optimální vyhlazovací parametr vzhledem k AMISE pro odhad hustoty pravděpodobnosti s jádrem řádu k je tvaru . 2fe+i »ofc y^-i °^'k 2knV(f(k))' tj. Aopt.o.fe = 0(n-1/(2fe+D), AMISE f(;hopt,0,k) = 0(n-2k/(2k+^). Metody pro odhad optimální hodnoty vyhlazovacího parametru h • metoda referenční hustoty hREF = {-^jiTj Sok(7n 2fc+1' • metoda maximálního vyhlazení metoda křížového ověřování r 2 n ~ hcv = arg min GV(h) = / f2(x, h)dx--V f-t(Xu h) heHn J nz—' i—l • iterační metoda iV(K) hiT = řešení rovnice: h - „-——-— = 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 /iopt,o,4- 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) = — x2)2 pro x e [—1,1]. Pomůcka: V(f") = 22,5. 4. Dokažte: a) a\ = x2gk(x) dx = 2k+5' 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(^-aľ) má rozptyl a2. 5. Dokažte vztah (3.15). 6. Spočítejte /ims pro • 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ů. 53 příslušnou hustotou f (x), pak platí • 0 < F (x) < 1, • F (x) je zprava spojitá, • lim^-oo F (x) = 0, lirn^oo F (x) = 1, . F(x) = ffoof(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 F, 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 " 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 9x2 1 (x-2)2 fix) = 0,5^= e"— +0,5^= e--— . 2V^ V^Ťř 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. 54 Užijeme-li substituce y = (t — Xi)/h, dostaneme F(X,h) = ^J2Í " K{y)áy=l-YjWÍX Xl To znamená, že odhad F v bodě x e M je definován takto ' x — Xi F(x,h) = -J2w h W(x) = / K(t) dt. (4.1) Zde předpokládáme, že K e Sq2,K(x) > Oprox € [—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(a;)dx < W(x) dx = 1, 3. f\w(x)K(x)dx = \, 4. ^(x)!^) dx = \ (\ - 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 ÍO x<-\, W(x) = <\(-x3+ 3x + 2) \x\<í, 1 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ě iel: EF(x,h)= [W(^-y)f(y)dy = 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 li počítáme metodou per partes a využijeme vlastnosti funkce W(t) h = h j W(t)f(x - ht) át u = W(t) v! = W'(t) = K(t) v' = f(x - hť)h v = -F(x - hť) = [-W(t)F(x - hť)]1^ + J F(x- hť)W'{ť) át = -F(x -h)+ j K(t)F(x - hť) át. (4.2) J-i Dále použijeme Taylorův rozvoj F(x - hť) = 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 - hť) át, uvažujeme-li substituce x — ht = z, dostaneme ř-00 px — h = - / f(z)áz= / f(z)áz = F(x-h). J x — h J — 00 (4.3) Vychýlení odhadu je tedy tvaru bias 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) át + F(x - h) - F(x). Odtud plyne EF(x, h)= [ K(t)F(x - ht) át. 56 A dále (z Taylorova vzorce) EF(x, h) = J K (t) [F (x) - htF'(x) + -h2t2F"ix) + o(h2) ) dt = F{x)+l-h2p2{K)F"{x) + o{h2) = F(x) + o(h). Nyní dokážeme tvar rozptylu. 1 varF(x,h) = - EW . n \ \ h =o(h) x-X -EZW x-X h Zde E2W (3^) = (EW (3^)) = (F(x) + o(h)f. Počítáme tedy jen integrál (označme jej h): = |substituce: x — y = th\ _ 1 n W2(t)f(x-ht)dt + h I f(x-ht)dt =F{x-h) První integrál počítáme metodou per partes a máme h = - \-F(x - ht)W2(t)} * + - Fix- hť)Wiť)W'iť) dt + -Fix - h) n n J n = --Fix - h) + - [ Fix - hť)Wiť)Kiť) dt + -Fix - h), n n J n použijeme nyní Taylorův rozvoj funkce F 2 W(t)K(t) (F(x) - htF'(x) + o(h)) dt -F(x) / W(t)K(t)dt--hF'(x) n J_x n tW(t)K(t) dt +o = \(vlasnost 3) =k(1~S-iW2 í2"') dx) (vlastnost 4) užitím vlastností funkce W a F'(x) = f(x) dostaneme F(x)-hf(x) ( 1 - / Wz{t)dt 'h ■o|- n Rozptyl je tedy tvaru var F(x, h) = F(x)-hf(x)[í- j W2(t)dt h -F(x)(í - F(x)) - -f(x) 1 - / W2{t) dt + o - + 0[-)--{F{x)+o{h)f 1 n I n 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: MSEŕYx, /i) = —— — —hf(x) (1 — f W2(ť)dt V 7-i , (44) + i(F"(x))2/l4/322(^)+0^+/l4N 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") = J (F" (x))2 dx < oo, K e S02, Ľm^oo h = 0 a lim„_!.00 nft = oo. Pak MISEF(-,/i) = - ÍF(x)(l-F(x))dx-ci-+c2/i4 + o( -+/í4V (4.5) n J n \n / Cl = 1 - ^ T^2(í) dí, c2 = Í/32(^)l/(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(x)(í-F(x))dx-Cl-+ c2h4 . (4.6) n J n AIV(h) AISB(h) 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 = n"173 (^)= 0(n-1/3) (4.7) a pak AMISVF(-,hopt,0,2) = -JF(x)(l-F(x))dx- — ^ n~4/3. (4.8) 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 S02 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- f W2(x)dx = 1 - f — (-x3 + 3x + 2)2dx = 0,2571, J-i J-i 16 c2 = \p2(K)V(F") = i ~-3,1818. Pak ^,0,2 = 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 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 4 Volba jádra 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 S02, např. • Epanečnikovo K(x) = |(1 — x2)I[_11](x), • kvartické K(x) = ±§(1 - x2)2I[_ltl] (x), • triweight K(x) = - x2)3I{_1A](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 Xi. Pak hcv = arg min CV(/i), h£Hn přičemž Hn = [an-1/3, Ďn-1/3] 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 můžeme aplikovat Terrelovu větu 3.3 pro k = 1. V tomto případě je 59 a tedy kde cti = j" x2gi(x) dx = V(<7Í) = ^. Odtud plyne, že 0 je odhadem 0 (viz rovnice (3.13) a (3.14)). Poznámka 4.3. Hodnota /ims 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 = [hg, /ims], kde hg je nejmenší vzdálenost mezi po sobě jdoucími body Xi,i = 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í: 1 ŕ n = 100, 0 = 1,3426, /32(K) = -, a = 1 - / W2(x) dx = 0,2571. 5 J-i 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))2 dx = -J f" {x) f {x) dx. Tudíž se budeme dále zabývat odhadem funkcionálu A = J 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,A <= S24. Pak Ä = n-1 J2 f"(^h) = -~2^3 E E K(2) l—l l—l j—1 Xň — X á 60 kde podle vztahu (3.11) je Pak nopt,2,4 — -lu T—nopt, 0,4 004 Č2 = --r2(K)íi Shrnutím předchozích úvah dostaneme proceduru pro odhad distribuční funkce F: Procedura 1. Najděte optimální vyhlazovací parametr /iopt,o,4 pro odhad hustoty s optimálním jádrem KoptflA e S04- 2. Najděte optimální vyhlazovací parametr /iopt,2,4 pro odhad druhé derivace hustoty podle vztahu (3.11) s k = 4 a optimálním jádrem e Sl24- 3. Vypočtěte odhad funkcionálu í/>i s využitím hodnoty /iopt,2,4 získané v kroku 2. 4. Vyčíslete optimální hodnotu vyhlazovacího parametru hpi = n -1/3 /_£L 1/3 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 hpi = 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 8 m3/s a vyšší se vyskytuje s pravděpodobností asi 0,3, průtok 20 m3/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 h^s = 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 =l- E w i^jr)' w{x) = /lm át- 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(-, h) = — F(x)(í-F(x))dx-c1-+ c2h n J n "--v- AIV(h) kde 4 y AISB(h) ci = 1 - y W2(t) dt, c2 = -^2{K)V{F"). Optimální vyhlazovací parametr vzhledem k AMISE pro odhad distribuční funkce je tvaru hopt,o,2 = n-1/3 , t.j. /iopt.0,2 = Oin-1^). Metody pro odhad optimální hodnoty vyhlazovacího parametru h • metoda křížového ověřování hcv = argmin^g^ CV(h) CV(h) = -J2 í (/(-oo,,.](^í)-í?-í(^/í))2dx, 71 i=i • metoda maximálního vyhlazení • plug-in metoda hPI = n-^ {^A-) ' . \-A/322(K) 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 — i2)2/^^!] (x). 2. Dokažte vlastnosti 2, 3 a 4 funkce 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) = -JLe-i(*2Hy+D2)A±_ e-H-2+fe-2)2)+IAe-i(-2+fe-4)2) JK ' 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). (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, Yi],..., [Xn, Yn]) z dvourozměrného spojitého rozdělení s hustotou f(x, y). Jádrový odhad hustoty / v bodě [x, y] G M2 je definovaný vztahem -t n i n f(xlVlU) = - Y^Kn(x - XllV - Yt) = — Y^K^ix - XllV - Y,)1 n n\H.\ '-^ V i=l 1 1 i=l (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 So2)- Existují dva typy těchto jader: • součinové jádro Kp(x,y) = K\(x) ■ K\(y), 1 Používáme zde zkrácený zápis pro dvourozměrnou hustotu normálního rozdělení, a to N(fii, fi2; uf, a\, q) 67 • sféricky symetrické jádro Ks(x,y) = cfc 1 iri (\J x2 + y2) ,ck = J J K\ (y'x2 + 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) = A(l-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 odhaduje analogická konstrukci jednorozměrného odhadu. Tedy v každém bodě [Xi, Yi] sestrojíme jádro Ku a odhad v bodě [x, y] je průměr n hodnot jader v tomto bodě - viz obrázek 5.4. 3 Statistické vlastnosti jádrových odhadů hustoty Stejně jako u jádrových odhadů jednorozměrných hustot můžeme kvalitu jádrového odhadu hustoty popsat lokálně pomocí střední kvadratické chyby: MSE/(z, y, H) = i ((K^ * f)(x, y) - (Ku * f)2(x, y)) +((KU * f)(x y) - f(x,y)f, 68 Tabulka 5.1: Třídy vyhlazovacích matic nebo globálně pomocí střední integrální kvadratické chyby MISE/(-,H)= ÍÍMSEf(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 \ diag(/i2, h2) = splňují následující předpoklady: \0 h2) (i) Nechí H = H„ je posloupnost matic takových, že {nh\h\^)-x a prvky h\ a h\ konvergují k nule pro n —> oo. (íí) Dále nechí všechny druhé parciální derivace funkce f jsou ohraničené, spojité a integrovatelné se čtvercem. (ííí) Jádro K splňuje J J xK(x, y) dx dy = J J y K (x, y) dx dy = 0, x2K(x,y)dxdy = // y2 K (x, y) dxdy = f32{K). Pak\ MISE /(■, H) = AMISE /(■, H) + o(hj + h42) + o((/ll/l2n)-1), 69 77 kde AMISE/(■, H) = ¥¥9- + l-f5l(K)(h\V(fxx) + 2h\h\V(]xy) + 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\, fi2, pro které AMISE nabývá minimální hodnoty, jsou dány vztahy: h =| _V^{fyy)V{K)_ hopt \nß^K)V^(fxx)[V(fxy) + V^fxx)V^fyy)] 1 ' 1 ' W\fxx)V{K) '2, opt nß2(K)V^(fyy) [V(fxy] + V^(fxx)V^(fyy)] (5.4) Z těchto vztahů plyne, že množina přípustných vyhlazovacích parametrů je tvaru a\n 1/6 < h\ < Ďin-1/6, o^n-1/6 < hi < b2rnT1/& pro vhodné konstanty 0 < a\ < b\ < oo, 0 < a2 < 62 < °°- Ukázkový příklad 5.3. Předpokládejme, že známe tvar hustoty f(x,y) = ^ e-3^' +v ^2 pro [x,y] G M2. Vypočítejme hodnoty optimálních vyhlazovacích parametrů pro odhad se součinovým Epanečnikovým jádrem. Podle vztahů (5.3) a (5.4) potřebujeme spočítat výrazy2 V(fxx) = JJ f2x(x, y)dxdy = JJ fxxxx(x, y)f(x, y) dx dy = 34 • (2^ V(fxy) = JJ f2y(x,y)dxdy = JJ fxxyy(x,y)f(x,y) = 33 • (2^)-\ V(fyy) = JJ fyy(x,y)dxdy = JJ fyyyy(x, y)f(x, y) dx dy = 34 • (27r)_1 Pro součinové Epanečnikovo jádro K(x, y) = — x2)(l — y2) platí V(K) = JJ K2(x,y)dxdy = 0,36, I32(K) = JJ x2K{x,y)dxdy = JJ y2K{x, y) dx dy = 0,2. Pak pro optimální vyhlazovací parametry máme h\ 6__(34 ■ (2^)-1)3/4 ■ 0,36_ opt ~ n • 0,22(34 • (2tt)-1)3/4[33 • (2tí)-1 + (34 • (2tt)-1)1/2(34 . (27r)-i)i/2] ' ,6 =_(34 ■ (2^)-1)3/4 ■ 0,36_ 2'opt n • 0,22(34 • (2tt)-1)3/4[33 • [2ti)-1 + (34 • (2tt)-1)1/2(34 . (27r)-i)i/2] ' 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 HOpt,0,2 = diag(/i2, h2,) = 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 (\X\, Y\],..., [Xn,Yn]) pochází z normálního rozdělení s hustotou a jádro K je dvourozměrnou standardizovanou normální hustotou, tj. K(x,y) = — e 2 2 . Pak podle metody referenční hustoty lze získat tento odhad vyhlazovací matice 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: /0,1292 0 \ íIref — y 0 0,9945J Na obrázku 5.6 je vykreslen odhad hustoty s pomocí Gaussova jádra a porovnání odhadu se skutečnou hustotou. 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)= //(/'(x,y,H))2dxdy--^/_í(Xí,yí,H). •* •* U i=l kde 1 " f-i(Xi,Yi, H) = —- J2 KM - Xj,Yí - Yj) 71 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 C V a MISE, který uvádí následující věta. Věta 5.2. Funkce CV je nevychýleným odhadem MISE, tj. platí ECV(H) = MISE/(-,H) - J J f2{x,y)áxáy. Důkaz. Vypočtěme střední hodnotu C V: EC\(H) = E Jí p(x,y,H)dxdy-^JTEf_l(Xl,Y,H.) l — l = E jj f2(x,y,-H)dxdy-2EKIl(X1-X2,Y1-Y2) = E II f2(x,y,ií)dxdy -2 //// KH(x - u, y - v) f (x, y) f (u, v) dxdy du dv 72 a úpravou MISE MISE/'(-,H) = i? JJ(f(x,y,H)-f(x,y))2dxdy = E J J f2(x,y,H)dxdy-2E J f(x,y,H)f(x,y)dxdy + J J f{x,y)dxdy = E JJ f2(x,y,H) dxdy-2 JJJJ KH(x - u, y - v) f (x, y) f (u, v) dxdy du dv + J J f2(x,y)dxdy. Porovnáním upravených výrazů dostaneme tvrzení. □ Optimální matici vyhlazovacích parametrů vzhledem k metodě C V označíme Hcv/ tj. Hcv = arg mm 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: /1,4532 0 tlcv — y 0 0,1431 Na obrázku 5.7 je vykreslen odhad hustoty s touto maticí. k 1, k x V x *, x : « > in, n (a) Data -3-2-1 0 1 2 3 (b) Odhad s Hc -3-2-1 0 1 2 3 (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 The 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/sowc03/tables/tablel.html 73 (b) Odhad s Href a Gaussovým jádrem (c) Odhad s Hcv a Epanečnikovám jádrem 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 (1145,2 0 \ /674,14 0 íIref — h-cv — y 0 24,9) y 0 59,35 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í: /270,5 0 \ /l797,2 0 íIref = ilcv = y 0 1516,5 y y 0 892,7 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 f(x, y, H) = — K^ií-^x -Xi,y- Y)1 i=l Dva typy jader: • součinové jádro: Kp(x, y) = Ki(x) ■ Ki(y), • sféricky symetrické jádro: Ks(x,y) = c^xKi (\Jx2 + y2),Ck = JJ K\ (\Jx2 + y2) dxdy. Asymptotická střední integrální kvadratická chyba dvourozměrného jádrového odhadu AMISE/(■, H) = + l-f32(K)(h\V(fxx) + 2h2h2V(fxy) + h2V(fvv)). Optimální vyhlazovací parametry vzhledem k AMISE _V3/4(fyy)V(K) °pt 1 n%(K)W*(fxx) [V(fxy) + y/V(fxx)V(fyy)\ , I V3'\fxx)V{K) "2,opt nP2(K)V^(fyy) [V(fxy) + VV(fxx)V(fyy)\ Metody pro odhad optimálních hodnot matice vyhlazovacích parametrů H = diag(/ii, h^) • metoda referenční hustoty /íj.REF = difl~1/&, i =1,2, • metoda křížového ověřování Hcv = argniinCV(H), CV(H) = f f (f(x, y, H))2 dxdy - - V f^X,, Y, H). H.eT> J J n A—' i—l 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 («11 «12 «21 «22 \«nl «n2 «lí «2í Mu «n > 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 F(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 = argminHex> 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 -1 2 -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 Ä"7 má nosič [—7,7]- Do funkcionálu T dosadíme jádro K1{t): T(KS) k—h K2s{t) át -ä tkKs{ť) dt -ä k—h 2v+l tkJ^K (í\dt -ä 2v+l 1 «f2(:r) dx k—ľ 1 , 5fe 2^+1 ÍSľ2(aľ) dx /c—^ aľfeíí(aľ) dx 2iv+l 1 S(2u+l)(k-u) 5{k-v){2v+ľ) T{K) Cv. 4 Výraz {K * 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 = J K (z) J 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 = j, • 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í, • 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 Sq2- Nakonec, s využitím vztahu pro druhou mocninu funkce cosx, tj. cos2 x = (1 + cos2x)/2, 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: E"=i ^(^0 - ^j) E"=i -kTl-h,h] (xo - ^j) _ ^[-fe.fe] (^0 - ^) _ I[-h,h](xo ~ xi) Ših-hM {xo - Xj) nI{_hM (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 [—/i, /i]. Priestleyův-Chaův odhad: n 1 ^ 1 1 n mPC(x,h) = - Kh(x - Xí)Yí = - ^ľl-h,h] (x - xiWi = -^^2YA-h,h](x ^ i—l i—l i—l Princip je podobný jako u Nadaraya-Watsonových odhadů (str. 15). Cv.2 mNW{x,h,K) - -——1 ,x-x A -\h- h b\- 1 K(x-Xis Eľ=i Y Eľ=i i^Ks (2^} y. V™ ±kK ( ±t^\ Ei=l ft*-^ (^fe*2"1) Cv. 3 Do vztahu (2.13) dosadíme hodnoty V(K) = 1,250 a /34(if) = -0,0476 (z tabulky 1.2) a z regresního modelu ?k+1 v2V{K){k\)2 l2,rZ ~\~ _L _ opt'°'fe 2knl32(K)V(m^)] 0,25- 1,25 -(4!)2 _ iopt.o.4 - 2 . 4 . 10Q . (_o,0476)2 • 32^ ~ ó' 'Ub U ' hopt,o,4 — 0,41. Cv. 5 Ve vztahu (2.7) použijeme Taylorův rozvoj i(x - hu) = mix) - hum'ix) H-----h 1 , / hkukmP^(x) + o(hk). 79 Využijeme vlastnosti jádra K e Sok a dostaneme h Eŕhix, h) = mix) + i-í)k—l3k(K)m^(x) + o(hk) + Oŕn"1). k\ Dále postupujeme jako při odvození vztahu (2.11). Cv. 6 Spočítejme nejdříve derivaci AMISEto(-, h) (viz vztah (2.12)) dAMKE^ = 2k dh nh2 (Ä;!)2 Jelikož hledáme minimum funkce AMISE fh(-,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,o,fe ct2V(K) 2k nh, neboli opt,0,k c72V(K) (kl) = 2k- h2J;tnkrk(K)V(m^) = 0, \\2 opt 1 J^h2kt^kp2{K)V{m nhOpt,0,k ArV(/iopt,o,fc) = 2fcAISB(/iopt,o,fe)- Kapitola 3 Otázka na str. 31 Obdélníkové jádro: Kh{z) = ^ pro z e [—h, h]. Odhad hustoty ^ n 1^1 1 n f(x,h) = -J2Kh(x-Xl) = -Y.^hh-hMÍx-Xi) = ^z\ZI[-hMÍx-Xi)- i=l Tedy v každém bodě Xj 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) x i-h Cv.l AMISE/(■,/») = + j^h2kp2(K)V(f^) = T(K) V(K) 1 2k rk(K) nhT(K) (k\)2 T(K) Pomocné výpočty: T(K) = S2k(32(K), S2k+1 = V(K) _ V(K) V(f^) V(K) T(K) 52kp2{K) 52k P2{K) = P2{K) = 1 T(K) 52k(32(K) 52kk r2fc+l °0fc _ r - — 00k x i+h 80 Celkem AMISE/(-,/i) =T{K) ( Sok i h2kV(f(kk \nhT{K) (kl)2 S2k Cv. 2 Pro danou hustotu f(x) = 20x(l — x)3 spočítáme její čtvrtou derivaci: f^(x) = —480, pak y(f^) = 4802. Dosazením do vztahu (3.5) s využitím vlastností jádra K e So4 dostaneme ' 2fc+l h €+1(fcQ2 opt'°'k 2knV{f(k)y 5o94(4!)2 1,25 (-0,0476)' r(4!)2 Opt,o,4 2 • 4 • 50 • 4802 2 • 4 • 50 • 4802 hopt,o,4 = 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) o2 = J\ x2gk(x) dx, gk(x) = Ak(l - x2)k+1 a Ak = (Jgt'iU pak a'i = Ak j x2(í - x2)k+1 dx = Ak j x- x(í - x2)k+í dx = Ak J(í-x2)k+1dx = /~ I ~.~m_~2^+l u = x u' = 1 v' = x(í - x2)k+1 v = (1 - a:2)fe+2 -1 2(k + 2) (í-x2) 2\fc+2 Ak 1 2(k + 2) (í-x2)k+2dx 1 (í-x2)k+2dx Víme, že Odtud plyne, že a tedy 2(k + 2) j_ľ í Ak+1(l - x2)k+2 dx = 1 (gk+1 je hustota). J-i"-.-' =gk+1{x) (í-x2)k+2dx A, fe+i t7k-Ak'2(k + 2)' Ak fe+i (2Jfc + 3)! ((k + í) + í)\222(-k+1'>+3 (k + 1) |222fe+3 2(k + 2) (2(k + 1) + 3)! (2fc + 3)!(fc + 2)!222fc+5 _ 1 (k + í)\222k+32(k + 2)(2k + 5)! ~ 2k + 5' 81 b) x —g ( —x dx = a ) a n \ a ,-;^g(^.r) dx-í t2—g(ť)—dt=a— í ťg(t)dt = cr J 1. M = \ Eti = 1 + 1 E?=i Vlastnosti T4^): x + 1 W2(x) dx W(x) dx ^±^d,= i 4 4 (x + 1)3 x+ 1 dx = -2 2 (x+l)s Cv. 1 Podle definice stačí integrovat polynom, který odpovídá kvartickému jádru, tj. W(x) 15 16 2 „ 1 . t--t3 + -ť — (x - -x3 + -x5 + —) = — (3x5 - ÍOx3 + Í5x 16 15' 16 Cv. 2 • Vlastnost 2: platí W'(x) = K(x), W(x) dx v! = 1 u = x v = W(x) v' = W'{x) = K{x) [xWix)]1^ - I xK{x) dx = l- W(í) - (-1) • W(-í) = 1. 82 Protože 0 < W{x) < 1, pak platí j\ W2{x) dx < j\ W(x) dx = 1. • Vlastnost 3: W{x)K{x) dx = u = W(x) v! = K(x) v' = K(x) v = W(x) = [W2(x)]1_1- j K(x)W(x)dx, =T =T W(í) = laW(-l) = 0 T = 1 - T => T W(x)K(x)dx = -. xW(x)K(x) dx u = x u' = 1 v' = K(x)W(x) v = \W2(x) V\w2(x)\\ l-W2{x) dx=l--l-J^ W2{x) dx l- (l-£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(-, h) 1 A ,3 n -= -ci - + Ac2h3 = 0. dn 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\ W2(x) dx = 0,7835. Dosadíme do vztahu (4.9) a dostaname 1/3 frM^n-W7'^-0',7835^ V7-d = 4,5089 -d-n^ 15 49 Kapitola 5 Cv. 1 Hodnota ck = J J Ki (x2 + y2) dx dy pro kvartické jádro je tt/3. Cv. 2 Zderivujeme vztah (5.2) pro AMISE a dostaneme dAMISE/(-,H) V(K) 1 0 = 0 = dhi nh\h2 ' 4 dAMISE/(-,H) V(K) 1 + -B2(K)(Ah\V(fxx) + Ahxh22V(fxyj) (6.1) + -p2(K)(4hlh2V(fxy) + ^hlV(fyy)) (6.2) d/i2 nh\h\ 4 První rovnici vynásobíme /ii, druhou rovnici vynásobíme h2 a odečteme je 0 = -AP22{K)(Ah\V{fxx) + Ah\h22V{fxy)) - \p22{K)(Ahlh22V{fxy) + Ah42V(fyy)), odtud plyne ±h\V{fxx) + Ah2h2V(fxy) = Ah2h2V(fxy) + Ah4V(fyy), K = V(fyy) h\ V(fxx) ■ Nyní dosadíme do rovnice (6.1) hf = h2- yjj^r vypočítáme h2 a pak h\. 83 Cv. 3 Rovnici (6.1) vynásobíme h\, rovnici (6.2) vynásobíme h2 a rovnice sečteme ^ffi = 4 • -Bl(K) (h\V{fxx) + 2hlhjV(fxy) + h42V(fyy)) 2-AIV 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 6-12 -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 V X V X V X V X V 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 V X V X V X V X V 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] [-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;y] [x;y] [x;y] [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 y X y X y X y X y 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. [x;y] [x;y] [x;y] [x;y] [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] [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 v = 0 k 2 -f(x2-i) 4 ±J(x2-l)(7x2-3) 6 -f|(x2-l)(33a;4-30a;2 + 5) 8 ^(x2 - l)(715x6 - lOOlx4 + 385x2 - 35) 10 ~lĚk(x2 ~ l)(4199x8 - 7956x6 + 4914x4 - 1092x2 + 63) 12 524288 l)(52003x10 124355x8 + 106590x6 39270x4 + 5775x2 231) v = 1 k Kopt,!* = 3 f x(x2 - 1) 5 -^§x(x2 -í)(9x2 -5) 7 ^x{x2 - l)(143x4 - 154x2 + 35) 9 -^x(x2 - l)(1105x6 - 1755x4 + 819x2 - 105) 11 ||2||a;(a;2 - l)(6783x8 - 14212x6 + 10098x4 - 2772x2 + 231) 13 ^2%x{x2 l)(260015x10 676039x8 + 646646x6 277134x4 + 51051x2 3003) 99 Tabulka 8.2: Optimální jádra pro v = 2 a v = 3 v = 2 k Kopt,2,k = K(2) 4 -±§(x2 - l)(5x2 - 1) 6 ^(x2-l)(77x4 - 58a;2 + 5) 8 -f|§(x2 - 1) (1755a;6 - 2249a;4 + 721a;2 - 35) 10 ^§§{x2 - 1) (3553a;8 - 6392a;6 + 3618a;4 - 672a;2 + 21) 12 - (x2 -1) (676039a;10 - 1562351a;8 + 1271974a;6 - 429726a;4 + 52899a;2 -1155) 14 T§§§š(x2 - 1) (884925a;12 - 2495270a;10 + 2653027a;8 - 1315028a;6 + 301587a;4 -26598a;2 + 429) v = 3 k Kopt,3,k 5 ?§x(x2 - l)(7x2 - 3) 7 1063495a;(a;2 l)(39a;4 38a;2 + 7) 9 a;(a;2 - l)(935a;6 - 1405a;4 + 597a;2 - 63) 11 -1Wřx(x2 - 1)(29393a;8 - 59432a;6 + 40018a;4 - 10032a;2 + 693) 13 22e2i!!x(x2 l)(382375a;10 969703a;8 + 895622a;6 364078a;4 + 61347a;2 3003) 15 - ^§§§f^x(x2 -1) (510255a;12 - 1554570a;10 + 1825625a;8 - 1034540a;6 + 288145a;4 -35178a;2 + 1287) 100 Literatura [1] Anděl, J.: Základy matematické statistiky. Matfyzpress, Praha (2005) ISBN 80-86732-40-1 [2] Forbelská, M., Koláček, J.: Pravděpodobnost a statistika I. Elektronický učební text, Masarykova unverzita, Brno (2012) http : //is .muni . cz/elportal/?id=1130308 [3] Horová, I., Koláček, ]., Zelinka, J.: Kernel Smoothing in Matlab. Theory and Practise of Kernel Smoothing. World Scientific, Singapur (2012) ISBN 978-981-4405-48-5 [4] Horová, L, Vieu, P., Zelinka, J.: Optimal Choice of Nonparametric Estimates of a Density and of its Derivatives. Statistics & Decisions 20, 355-378 (2002) [5] Horová, L, Zelinka, J.: Contribution to the bandwidth choice for kernel density estimates, Comput. Stat. 22, 31-47 (2007) [6] Köhler, M., Schindler, A., Sperlich, S.: A Review and Comparison of Bandwidth Selection Methods for Kernel Regression. Discussion Paper No. 95, Georg-August-Universität Göttingen (2011) [7] Koláček, J.: Jádrové odhady regresní funkce. Disertační práce, Masarykova univerzita, Brno (2005) http://is.muni.cz/th/19999/prif_d/ [8] Müller, H.-G.: Nonparametric regression analysis of longitudinal data. Springer, New York (1988) ISBN 978-0-387-96844-5 [9] Müller, H.-G.: Smooth optimum kernel estimators near endpoints. Biometrika 78, 521-530 (1991) [10] Scott, D.W., Gorry, G.A., Hoffman, R.G., Barboriak, J.J., Gotto, A.M.: A new approach for evaluating risk factors in coronary artery disease: a study of lipid concentrations and severity of disease in 1847 males, Circulation 62,477-484 (1980) [11] Scott, D.W.: Multivariate density estimation: Theory, practise, and visualization. Wiley, New York (1992) ISBN 0-471-54770-0 [12] Silverman, B.W.: Density estimation for statistics and data analysis. Chapman and Hall, London (1986) ISBN 0-412-24620-1 [13] Vopatová, K.: Volba vyhlazovacích parametrů pro jádrové odhady hustot. Disertační práce, Masarykova univerzita, Brno (2011) http :/ / is .muni . c z /th/ 63 985/prif_d/ [14] Wand, M.P, Jones, M.C.: Kernel Smoothing. Chapman and Hall, London (1995) ISBN 0-412-55270-1 101