Ú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í „hhovoř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 šířce vyhlazovacího okna, 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. Příslušný toolbox pro program Matlab je dostupný na adrese: https://www.math.muni.cz/veda-a-vyzkum/vyvijeny-software/ 2 7 4-matlab-toolbox.html V příloze (kapitola 6) 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 Vopatové, Ph.D., za pomoc při sazbě tohoto textu a za příspěvek ke kapitole 5 a podkapitolá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 10 1 Motivace........................................... 10 2 Základní typy neparametrických odhadů........................ 11 3 Statistické vlastnosti jádrových odhadů......................... 15 4 Volba jádra.......................................... 19 5 Volba vyhlazovacího parametru.............................. 20 5.1 Metoda křížového ověřování........................... 20 6 Automatická procedura .................................. 22 7 Aplikace na reálná data................................... 24 3 Jádrové odhady hustoty 28 1 Motivace........................................... 28 2 Základní typy neparametrických odhadů........................ 28 3 Statistické vlastnosti jádrových odhadů hustoty .................... 29 3.1 Odhad derivace hustoty.............................. 33 4 Volba jádra.......................................... 34 5 Volba vyhlazovacího parametru.............................. 34 5.1 Metoda referenční hustoty............................. 34 5.2 Metoda maximálního vyhlazení......................... 35 5.3 Metoda křížového ověřování........................... 37 6 Automatická procedura .................................. 39 7 Aplikace na reálná data................................... 41 4 Jádrové odhady distribuční funkce 45 1 Motivace........................................... 45 2 Základní typy neparametrických odhadů........................ 45 3 Statistické vlastnosti odhadu ............................... 47 4 Volba jádra.......................................... 49 5 Volba vyhlazovacího parametru.............................. 49 5.1 Metody křížového ověřování........................... 49 5.2 Princip maximálního vyhlazení.......................... 50 5.3 Plug-in metoda................................... 50 6 Aplikace na reálná data................................... 51 2 5 Jádrové odhady dvourozměrných hustot 55 1 Motivace........................................... 55 2 Základní typy odhadů................................... 56 3 Statistické vlastnosti jádrových odhadů hustoty .................... 57 4 Volba jádra.......................................... 59 5 Volba vyhlazovacího parametru.............................. 59 5.1 Metoda referenční hustoty............................. 59 5.2 Metoda křížového ověřování........................... 60 6 Aplikace na reálná data................................... 61 6 Přílohy 64 1 Symbolika O, o ....................................... 64 2 Datové soubory....................................... 65 3 Seznam použitého značení h vyhlazovací parametr H matice vyhlazovacích parametrů to(-) regresní funkce /(•) hustota pravděpodobnosti F(-) distribuční funkce h odhad vyhlazovacího parametru a odhad směrodatné odchylky a m odhad regresní funkce m f odhad hustoty / F odhad distribuční funkce F J značí integrál , pokud není uvedeno jinak K(-) jádrová funkce (jádro) V(K) V(K) = / K2(x)dx pk{K) pk{K) = / xkK{x)áx f * g konvoluce funkcí / a g, (/ * g){x) — J f(t)g(x — ť) dí W(-) intergál z jádra, W(x) = f*^ K(t) dí NW Nadarayovy-Watsonovy odhady PC Priestleyovy-Chaovy odhady LL lokálně lineární odhady GM Gasserovy-Múllerovy odhady 4 M SE střední kvadratická chyba (mean square error) MISE střední intergá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 intergated 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 5 Kapitola 1 Jádrové funkce a jejich vlastnosti Výstupy z výukové jednotky Student • bude znát základní třídy jádrových funkcí, jejich vlastnosti a metody jejich konstrukce. 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 vlast- Definice 1.1. Nechť i/, 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,ye [-1,1],L>0, 2. nosič(.?f) = [-1,1], tj. K = 0 vně intervalu [-1,1], 3. K splňuje momentové podmínky: Taková funkce K se nazývá jádro řádu k a třída všech těchto funkcí se označuje Svk- Příklad 1.1. V tabulce 1.1 jsou uvedeny příklady několika jader společně s jejich grafy. Přitom funkce I[-i,í\ (x) je indikátorová funkce intervalu [—1,1], tj. nosti. (1.1) a xkK{x) dx ^ 0, tuto hodnotu označíme f3k(K). Nyní uvedeme dva typy jader - jádra s minimálním rozptylem a optimální jádra. 6 Tabulka 1.1: Jádra S 02 obdélníkové jádro K(x) = l(l-x2)I{_1A](x) Epanečnikovo jádro K(x) = &il-X2)2I[_ll]ix) kvartické jádro 1.1 Jádra s minimálním rozptylem Předpokládejme, že K e Svk, () < v < k — 2, v ak jsou obě sudá nebo lichá1. Uvažujme funkcionál V(K) — K2{x) dx a zabývejme se problémem najít takové jádro K e Suk, pro které tento funkcionál nabývá minimální hodnoty, tj. řešíme variační úlohu min V(K) za předpokladu K e Svk- Řešení této úlohy se nazývají jádra s minimálním rozptylem, což jsou polynomy stupně k — 2 na intervalu [—1,1]. Tyto polynomy jsou sudé funkce pro k sudé a liché funkce pro k liché. Mají k — 2 různých kořenů v intervalu (—1,1). Obecný vztah pro jádra s minimálním rozptylem lze nalézt ve [4]. Příklad 1.2. Jádra s minimálním rozptylem: S02: K(x) 24: K(x) 2-4-1,1 (x) K(x)^-f(í-3x2)I[_1A](x) Poznámka 1.1. Jádra s minimálním rozptylem mají skoky v koncových bodech intervalu [—1,1], což negativně ovlivňuje hladkost výsledného odhadu. obvykle se používá pojem parita, tedy v a k mají stejnou paritu. 7 Tabulka 1.2: Optimální jádra pro ľ — 0,1,2 v = 0 k C>0fc Kopt,0,k 2 1,7188 -|(x2-l) 4 2,0165 i|(x2-l)(7a;2-3) 6 2,0834 -i2|(x2-l)(33x4 - 30a;2 + 5) 8 2,1021 j^(x2 - l)(715x6 - lOOlx4 + 385x2 - 35) v= 1 k Sik Kopt,i,k 3 1,4204 ^x(x2 - 1) 5 1,7656 -^x(x2-l)(9a;2-5) 7 1,8931 ^x(x2 - l)(143x4 - 154x2 + 35) i/ = 2 $2k Kopt,2,k 4 1,3925 -±§(x2 - l)(5x2 - 1) 6 1,6964 ^(x2-l)(77x4 - 58a;2 + 5) 8 1,8269 -||§(a;2 - l)(1755x6 - 2249x4 + 721x2 - 35) 1.2 Optimální jádra Při vyšetřování statistických vlastností se setkáváme s následujícím funkcionálem T(K) xkK(x) dx 2v+l K2(x) dx k—u\ 2fe+l V(K) který lze zkráceně psát T(K) = {\pk{K)\2v+1V{K)^-^)2'^^. 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ů. Příklad 1.3. Optimální jádra: * e [-1,1] Důkaz viz [4]. Jako příklad lze uvést jádro Koptfl^{x) — |(1 — x2) e So 2 a -KTÓpt 0 2(^) = Ki,3(x) = -|x e Si3. 2. Nechf e S„fc, v e N, pro 5 > 0 položíme Jádra K, Kg nazýváme ekvivalentní. Dokažte: Funkcionál T(K) je invariantní vzhledem k této transformaci, tj. T{K) = T{KS). 3. Nechť jsou dány funkce /ag, pro které platí J f2 (x) dx < 00 a J g2 (x) dx < 00. Konvoluci f * g definujeme vztahem (f*g)(x) = J f(t)g(x-t)dt. Vlastnosti konvoluce • f*g = g*f, • f * (g *h) = (f * g)*h, • f*(g + h) = f*g + f*h. Ukažte, že pro jádrovou funkci K e S02 platí vztah 1 3 — 0, (K * K)(x)x3 dx = { 0 j = 1, 2/32(^) j = 2. 9 Kapitola 2 Jádrové odhady regresní funkce Výstupy z výukové jednotky Student • bude znát základní typy jádrových odhadů regresní funkce a jejích derivací. • bude schopen analyzovat statistické vlastnosti odhadů. • bude mít přehled o metodách pro volbu vyhlazovacího parametru. • se seznámí s automatickou procedurou pro simultánní volbu vyhlazovacího parametru, jádra a jeho řádu. • bude schopen analyzovat daný soubor dat a aplikovat uvedenou proceduru na tento soubor. • bude schopen použít příslušný toolbox v Matlabu a zkonstruovat odhad regresní funkce. 1 Motivace Předpokládejme, že pro pevné nebo náhodné hodnoty nezávisle proměnné X máme k dispozici naměřené hodnoty závisle proměnné Y. Chceme-li tato data analyzovat, musíme nalézt vhodný funkční vztah mezi těmito proměnnými. Jestliže dvojice bodů (xí,Y),í — l,... ,n, znázorníme graficky, pak pouhý pohled na takový 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 Yi — m(xi) + Si, i — l,...,n, (2.1) kde m je neznámá regresní funkce, xir i — 1,..., n, jsou body plánu a e„ 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 Esí — 0, var Si — u2, i — l,...,n. (2-2) Poznámka 1.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 Xi,..., Xn jsou náhodné veličiny se stejnou hustotou /, 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. 10 Bez újmy na obecnosti budeme v dalším předpokládat, že pro body x i, í — 1,..., n, platí 0 < x\ < x2 < • • • < xn < 1. Cílem regresní analýzy je nalézt vhodnou aproximaci to neznámé funkce to. 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 1.1. Obr. 2.1 ilustruje na simulovaných datech nevhodnost aplikace parametrického přístupu. V tomto případě byla data generována podle vztahu sin Attxí Y —_-__|- e ■ * (1 + cos Ofiirxi)2 kde body Xi — l/100,z — 1,..., 100, a chyby eí, í — 1,..., 100, mají normální rozdělení ÍV(0; 0,25). (Data jsou v tabulce 6.1.) 2 1 0 -1 -2 0 0.2 0.4 0.6 0.8 1 Obrázek 2.1: Simulovaná data (x) s regresní přímkou (—) a původní funkcí (—) 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.1 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í funkcí), 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í 11 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 subintervalu, tj. ŕh{x, h) E ^/[i^OO i=l_ n ' i=l kde /[g j je indikátorová funkce subintervalu B j. Výsledek aplikace regresogramu na simulovaná data příkladu 1.1 je znázorněn na obr. 2.2. Vidíme, že tento odhad „vhodně" vystihuje tvar funkce, ale výsledný odhad je příliš hrubý. 0.2 0.4 0.6 0.Í Obrázek 2.2: Regresogram (—) pro simulovaná data (x) z příkladu 1.1 s původní funkcí (—) 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. Přesněji rh{x, h) E Yil[x-h,x+h}{%i) i=l_ n E I[x-h,x+h]{xi) (2.3) Obr. 2.3 ilustruje aplikaci této metody na simulovaných datech příkladu 1.1. Uvedené metody 0.6 O.f Obrázek 2.3: Klouzavý průměr (—) pro simulovaná data z příkladu 1.1 12 patří mezi nejjednodušší neparametrické vyhlazovací metody. Jádrové odhady lze považovat za zobecnění těchto metod. Připomeňme zde základní myšlenku vyhlazování tak, jak ji formuloval R. Eubank v r. 1988: Jestliže předpokládáme, že m je hladká funkce, pak pozorování v bodech xi blízko bodu x obsahují informace o hodnotě m v bodě x. Bylo by tedy vhodné užít lokálních průměrů dat blízko bodu x, abychom získali odhad m(x). Obecně lze jádrové odhady regresní funkce m v bodě x definovat takto n m(x,h)=J2W>(x>h)Y» (2-4) i=l kde funkce Wi,i — 1,..., n, se nazývají váhy, nezávisí na Y, 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 W, závisí na jádrové funkci K. Nechť K e Sofc, k je sudé číslo, položme Kh(-) — x^(tí)' Mezi nejznámější typy jádrových odhadů regresní funkce patří: 1. Nadarayovy-Watsonovy odhady (1964) n E Kh{x- xx)Y% mNW(x, h) = 2=±-; Y, Kh(x - Xi) i=l 2. Priestleyovy-Chaovy odhady (1972) 1 " -,{x, h) = - V" Kh(x - xl)Yl, 71 * * mpc{ n ■ i=l 3. lokálně lineární odhady (Stone 1977, Cleveland 1979) 1_ {^(x, h) - ši{x, h)(x - Xj)}Kh(x - Xj)Yt i=l kde 0, h) = - V ■ rí ^—' n 4^ š2{x, h)šo(x, h) — ši(x, h)2 1 " šr(x) = - y^(x - xl)rKh(x - Xi), i=l 4. Gasserovy-Múllerovy odhady (1979) n s* r(x,h) = '22Yi J Kh(x-t)dt, mgm \ kde So — 0, Si — (xí + Xj+i)/2, sn — 1. Tento odhad je konvolučním typem odhadu. Úmluva. Uvedené jádrové odhady budeme zapisovat ve tvaru n mj(x,h) = J2WíA)(x>h)Y» i=l kde váhy WJ;A) a A značí příslušný typ odhadu NW, PC, LL, CM. V dalším textu budeme zkráceně psát Wj. 13 V mnoha aplikacích je užitečný zejména Nadarayův-Watsonův odhad trmw- Popíšeme nyní jeho konstrukci a budeme ilustrovat vliv vyhlazovacího parametru na kvalitu odhadu. Pro daný bod xq, h < xq < 1 — h, jsou váhy Nadarayova-Watsonova odhadu dány vztahem Wi(x0,h)= Kh(£°-Xi) J2WÁ*o,h) = l. Obrázek 2.4 ilustruje konstrukci odhadu v bodě x0, který je založen na pěti pozorováních (xi, Yi), (x5,15) (černé křížky). Modrá parabola reprezentuje Epanečnikovo jádro Kh a modré kroužky znázorňují hodnoty vah Wi — Kh(xo — x i) j Y^l=\ Kh(%o — xi) pro i — 1,... ,5. Výsledný odhad regresní funkce fh v bodě x0 je označen červenou hvězdičkou. 1.51—1-1-1—1-1-1-1 Obrázek 2.4: Ilustrace Nadarayova-Watsonova odhadu v bodě x0 Jádrový odhad není definován pro £™=1 Kh(% — x i) — 0. Jestliže nastane případ „0/0", pak klademe mNw(x, h) — 0. Omezíme se nyní na odhady funkce m v bodech plánu Xi, í — 1,..., n. Pak pro h —>• 0 platí mNW{x1, h) ->• = li To znamená, že při malé šířce vyhlazovacího okna odhad reprodukuje data (viz obr. 2.5(a)). Dále, pro h —>• 00 platí n n ^(o)E^- x n mNW(Xi,h) -+ ^ = ngJ(o) = « i=i Tedy velká šířka okna vede k přehlazení, a to k průměru dat (viz obr. 2.5(b)). Na obrázku 2.6 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í, nebof i asymptoticky optimální odhady obsahují poměrně značné „množství šumu" a to nechává prostor pro subjektivní posouzení. Poznámka 2.1. Intervaly spolehlivosti pro hodnotu regresní funkce m 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 m v bodě x. Jsou definovány takto ,V(K)d2(x) _ V(K)Č2(x) 14 O 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (a) Podhlazený odhad, h = 0,02 (b) Přehlazený odhad, h = 0,4 Obrázek 2.5: Podhlazený a přehlazený odhad regresní funkce z příkladu 1.1 2 F x x 0 0.2 0.4 0.6 0.8 1 Obrázek 2.6: Optimální odhad, h = 0,0816 kde Mi-f je (1 — a/2)-kvantil standardního normálního rozdělení a odhad rozptylu v bodě x je dán vztahem n d2(x) ^^2wi(x,h)(Yl-fh(x,h))2. i=l Ukázka intervalu spolehlivosti pro a — 0,05 je na obrázku 2.7. 3 Statistické vlastnosti jádrových odhadů Kvalitu jádrového odhadu lze lokálně popsat pomocí střední kvadratické chyby. Budeme se zabývat asymptotickým tvarem této chyby, nebof pro neparametrické odhady, na rozdíl od odhadů parametrických, neexistuje nevychýlený odhad, tj. takový odhad že, Em — m pro skoro všechna i e M [2j. Věta 3.1. Nectí K e Sok,EY2 < oo a nectí posloupnost vyhlazovacích parametrů h — hn (n — 1,2,...) splňuje podmínky: hn —>• 0, nhn —>• oo pro n —>• oo. Pak v každém bodě spojitosti funkce m platí n fh(x, h) — Wi(x, h)Yi m(x), (2.5) i=l kde —^ značí konvergenci podle pravděpodobnosti (viz např. [1 ]). Uvedené odhady m jsou tedy konzistentními odhady m. 15 O 0.2 0.4 0.6 0.8 1 Obrázek 2.7: Interval spolehlivosti pro data z příkladu 1.1 při a — 0,05 Střední kvadratická chyba MSE odhadu m v bodě x je obecně dána vztahem MSEm(x, h) — E{rh{x, h) — m(x))2. Upravíme tento vztah — Em2(x, h) — 2m(x)Em(x, h) + m2(x) — (Em(x, h) — m{x)Ý + Erh2{x, h) — (Em(x, /i))2, (2-6) bias2 var což znamená, že střední kvadratická chyba může být vyjádřena jako součet rozptylu odhadu var m(x, h) a čtverce vychýlení bias2 m(x, h) (viz obr. 2.8). Tento rozklad rozptyl-vychýlení usnadňuje analýzu vlastností odhadu, kterou se budeme zabývat. 0 0.2 0.4h 0.6 h 0.8 1 opt Obrázek 2.8: MSE (—) jako součet rozptylu ( ) a vychýlení (—) Všechny uvedené jádrové odhady regresní funkce m^wr ffipci ™gm jsou asymptoticky ekvivalentní. 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 W^. Připomeňme, že pro Priestleyovy-Chaovy odhady je váhová funkce tvaru Wi{x, h) = Kh{x - x,) = Ik(^í) . 16 Pro další výpočty budeme předpokládat: (i) Jádrová funkce K je sudou funkcí na intervalu [—1,1], tj. K e So 2- (ii) Vyhlazovací parametr h — h(n) je posloupností 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é. Je zřejmé, že pro n —>• oo platí1 1 ™ Z"1 1 2ľ í Efh(x,h) =-^2Kh(x - Xi)m(xi) = I -K(^——^m(t) dt + Oin^1). (2.7) n i=l ^° Symbolika O a o viz kapitola 1. Nechf u — r£j^, odtud dí — — hdu a tedy s využitím Taylorova rozvoje /x/h K(u)m(x — hu) du + 0(n_1) -{\-x)/h ľx/h r ?/2/?2 i = / m(x) - uW(x) + ——m"(x) dw + o(/i2) + O^1). (2.8) J-(l-x)/h L 2 J Podle výše uvedených předpokladů platí h < x < 1 — h, tedy i//i^ooa- (1 — x) / /i —>• — oo pro /i —>• 0. Odtud, s využitím faktu, že nosičem funkce je interval [—1,1], plyne E^,h) = ™W / *<„)d. - / ■*(„) iu + / „>*(„) d, + „A') + 0(„-) r1 r1 h? r1 = m{x) J K(u)du-j uK{u)du+ —m"(x)J u2K{u) du + o{h2) + C^n"1). Celkem dostaneme h2 Eíh(x, h) = m(x) + —p2(K)m"(x) + o(/i2) + C^n"1) /i2 « m(x) + —j32{K)m"(x). Podobně pro rozptyl platí varm(x, h) — E(íh{x, h) — Erh{x, /i))2 /1 " 1 " v i=l =1 1 / " = — EÍ ^2Kh(x - xl)(Yl -m(xi)) z vlastností (2.2) plyne EKh{xi)Kh{xj)eiSj — 0 pro i ^ j, tedy 1 " = — Y.KUx-x^Ee2 i=l =^^)=ä(í>(^)*+o<»-'> 1 Jedná se o přibližný výpočet integrálu. 17 Opět s využitím substituce u — a vztahu 0(n 1) — o((nh) 1) můžeme pro n —> oo psát 2 r 2 pl 2 varm(i,/i) = / K2(u) du + o^nh)-1) = / K2(u) du + o^nh)-1) « ?—V(K). nh J nh J_1 nh Tímto jsme dokázali následující větu o tvaru střední kvadratické chyby. Věta 3.2. Nectí jsou splněny předpoklady (i) — (iii), pak střední kvadratická chyba nabývá tvaru MSEto(x, h) = ^V(K) + h4p2(K)^(m"(x))2 +o(h4 + (nh)-1). (2.9) 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í intergá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í MISE(/i) = / MSE(x, h) dx = AMISE(/i) + o(h4 + (nh)-1). Jo AMISE je tvaru a2 h4 AMISE(/i) = AMISEm(-, h) = ~^{K) + -jPl(K)V(m"), (2.10) AIV AISB kde V(m") = J1 (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ším cílem je minimalizovat AMISE(/i), tzn. najít takovou hodnotu vyhlazovacího parametru h, pro kterou AMISE nabývá minimální hodnoty, a tedy odhad bude nejlepší ve smyslu AMISE. Užijeme metody matematické analýzy a spočítáme derivaci dAMISE(/i) _ a2V(K) dh nh2 položíme ji rovnu nule a vyjádříme h + hifá(K)V(m"), o2V(K) Vt,o,2- np2{K)v{m„y ^A1> Poznámka 3.1. Tento výpočet vede k nalezení minima AMISE, protože platí d2 AMISE ^2—>°- Poznámka 3.2. Jestliže jádro K náleží do třídy Sok, pak AMISE je tvaru AMISE(/0 = ^V(K) + J^h2k(32(K)V(m^), a pro optimální vyhlazovací parametr h platí (2.12) fe+1 *>V(K)(k!)> °^°'k 2knP2(K)V(m^)y { ó> kde V(m^) — (nrSk\x))2 dx, podrobněji např. [8]. Nyní uvedeme důležité lemma, které ukazuje zajímavou vlastnost vyhlazovacího parametru. 18 Lemma 3.1. Pro /iopt,o,fe platí AIV(/iopt,o,fc) = 2A: AISBrA^.o.k)- Důkaz. Viz cvičení. □ Vztah (2.13) pro optimální šířku vyhlazovacího okna ukazuje, že řád konvergence optimální šířky vyhlazovacího okna /iOpt,0,2 závisí na řádu jádra k, tedy pro k — 2 je 0(n_1/5). Rovnice (2.13) má pouze teoretický charakter, protože hodnota /iOpt,0,2 závisí na neznámých veličinách cr2 a m"(x), a tedy není užitečná pro praktické účely. Na druhou stranu umožní vyjádřit asymptotickou rychlost konvergence AMISE(/i). Dosadíme-li (2.13) do vztahu (2.10) pro AMISE, dostaneme AMISE(/lopt,o,2) = n-4/5((r2l/(^))4/5(/322(^)(A;!)-2l/(m"))1/5. (2.14) Odtud plyne, že řád konvergence pro jádro řádu k je n~2fe+!, tedy s rostoucím k se zvyšuje asymptotická rychlost konvergence. Ale není zcela jasné, zda tato zvýšená rychlost konvergence vede již k zlepšení pro konečné rozsahy výběrů, neboť ostatní veličiny se rovněž mění s k. O těchto problémech pojednáme podrobněji v dalším odstavci. 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. 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 je plošší. Poznámka 3.3. Vyšetřování kvality odhadu obvykle probíhá za předpokladu, že pracujeme s vnitřními body intervalu [0,1]. V hraničních oblastech, tj. v intervalech [0, h) U (1 — h, í], je kvalita odhadu ovlivněna negativně skutečností, že jádro K zde nesplňuje momentové podmínky (1.1). To je způsobeno tím, že blízko krajních bodů nosič jádra K zasahuje do oblasti, kde nejsou žádná data, což zhoršuje odhad - viz obr. 2.9. Hraniční efekty jsou také patrné na obrázcích 2.6 a 2.12, 1.5,-1-,-1 X Obrázek 2.9: Hraniční efekt zejména u pravého okraje intervalu. Problém okrajových efektů lze překonat např. použitím hraničních jader [9] nebo reflexní metodou [4]. 4 Volba jádra Volba jádra není z asymptotického hlediska podstatná, jak je zřejmé z faktu (2.14). Je vhodné zvolit optimální jádro, které minimalizuje fukcionál T(K), neboť tato jádra jsou spojitá na M a odhadovaná regresní funkce tak „zdědí" hladkost jádra. Vhodná jsou jádra třídy So 2 a S04, lze je vybrat z tabulek pro Sofc. 19 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 ŕh-i{xi,h) — wÁxi, h)Yt, i — í,...,n. i=i Funkce křížového ověřování je definována takto 1 ™ CVt^-^ím.ifii./il-ľi)2 (2.15) i=l a odhadem optimální hodnoty vyhlazovacího parametru je bod, v němž nastává minimum této funkce, tj. hopt,o,k = /icv = arg min CV(h). Při praktických úlohách hledáme minimum na intervalu Hn — [afcn-1/(2fe+1); 6fcn-1/(2fe+1)]/ jehož tvar plyne ze vztahu (2.13), přičemž cifc, b]~ jsou konstanty, 0 < cik < bk < oo. Pro ekvidistantní body plánu byl na základě zkušeností odvozen interval [^,2]. Poznámka 5.1. Někdy se místo chyby MISE používá průměrná střední kvadratická chyba AMSE (average mean square enoi) 1 ™ AMSE — -Ej2(m(xt,h) - m(xt)ý n i=l Využívá se zejména v případech, kdy je nevyhovující použít numerické integrování související s chybou MISE. Věta 5.1. Střední hodnota funkce CV(h) je 1 ™ ECV(h) = -£y(m(i„/i)-ray)2V. n ^-^ i=l AMSE Důkaz. Funkci křížového ověřování lze rozepsat 1 ™ 2 CV(/i) — — ^(m_i(xi, /i) - m(xj) - (Yt - m(xt))) ^ n \ n \ n = - y^(m-i(xi, h) - m(xl)) - 2- V"(fh-i(xj, /i) - m(xí))eí + - V" e2 i—l i—l i—l — - '^2(m-l(xl, h) - m(xi))--£í (5ľ We(xl,h)Yi - m(xl)j + - y^e; n —' n —' \*—' / n i=l i=l i=l i=l 20 Pak střední hodnota E CV(h) je rovna součtu tří hodnot 1 ™ ))2 = AMSE, n n i=l i=\ n Druhý a třetí vztah plyne z vlastností (2.2). □ Tento výsledek znamená, že minimalizace EGV(h) odpovídá minimalizaci AMSE. Jestliže tedy předpokládáme, že minimum CV(/i) je blízko minima E CV(h), pak tato minimalizace dává dobrou volbu vyhlazovacího parametru - viz ilustrace na obr. 2.10. Obrázek 2.10: Porovnání minima AMSE (červenou) a minima funkce křížového ověřování CV (modrou) pro simulovaná data Příklad 5.1. Použijeme metodu křížového ověřování pro nalezení vyhlazovacího parametru pro data z příkladu 1.1. Při použití Epanečnikova jádra získáme vyhlazovací parametr hcv — 0,1158. Na obrázku 2.11 je zobrazen odhad regresní funkce s tímto parametrem. Obrázek 2.11: Simulovaná data (x) s jádrovým odhadem regresní funkce (hcv a původní funkcí (—) 0,1158) (-) 21 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ř. [4, 7]). 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 [4]), která odhadne všechny parametry tak, aby byla minimalizována AMISE. Procedura byla původně odvozena pro odhad hustoty pravděpodobnosti ([5]), ale lze ji aplikovat i pro odhad regresní funkce. Uvedeme zde její zjednodušenou verzi. Podle poznámky 3.2 je známo, že AMISE a /iOpt,0,fe jsou tvaru n2 1 n2h2k+1ih\\2 AMISE(„, = Z-VW + *»& = Dosadíme-li vyhlazovací parametr /iOpt,0,fe do chyby AMISE a použijeme vztahy (^g)^ T(K) = {pk(K)V\K)y dostaneme vyjádření AMISE AMISE^o,,) = ^h+1)SokT(K), v němž je vidět separace parametrů K, h a k, což umožňuje vybrat tyto parametry simultánně, Právě tento vztah je základem automatické procedury. Označme 2nkhoptfl^k a množinu vhodných řádů k označme /(fco) = {2j,j = 0,...,[^]}, kde [z] značí celou část čísla z. Procedura pak probíhá v pěti krocích: Krok 1 Pro k G I(ko) najděte optimální jádro Koptfl,k <= Sok, které je dáno tabulkou 1.2 a vypočtěte kanonický faktor S^k- Krok 2 Pro k G I(ko) a Koptto,k <= Sofc najděte optimální vyhlazovací parametr /iOpt,0,fe- Krok 3 Pro k 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. Krok 4 Vypočtěte optimální hodnotu řádu k, které minimalizuje funkcionál L(k). Krok 5 Použijte parametry z předchozích kroků k získání optimálního jádrového odhadu regresní funkce, tj. ™(x> ^opf,0,fc) = W*(X> hpt,0,k)Y>- i=l 22 Příklad 6.1. Aplikace procedury na data z příkladu 1.1. Maximální řád jádra zvolme ko — 8, tedy množina možných řádů jader je 1(8) — {0, 2,4, 6, 8}. Pro tyto řády spočítejme hodnoty z kroků 1-3, v kroku 2 jsme použili metodu křížového ověřování pro nalezení optimálního vyhlazovacího parametru hoptflJ;. k Kopt,0,k ^Ofe h L(k) 2 -1) 1,7188 0,1158 0,0648 4 15 IJ2 32 vx l)(7x2-3) 2,0165 0,2446 0,0575 6 105 f 2 256 v - l)(33x4 - 30x2 + 5) 2,0834 0,3575 0,0574 8 315 12 4096 v - l)(715x6 - - lOOlx4 + 385x2 - - 35) 2,1021 0,4543 0,0592 Z tabulky vidíme, že optimální řád jádra je k — 6. Výsledný odhad je uveden na obrázku 2.12. 2- * x x X -2[__,_,_,_,_: 0 0.2 0.4 0.6 0.8 1 Obrázek 2.12: Simulovaná data (x) s jádrovým odhadem regresní funkce při použití procedury ( ) a skutečnou funkcí (—) společně s 95% intervalem spolehlivosti 23 7 Aplikace na reálná data První 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 6.5 a na obrázku 2.13. 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 (tfopt.0,2), Voc = 0,1525 (^,0,12). Výsledné odhady na obrázku 2.14 ukazují, že metoda křížového ověřování spíše podhlazuje. 13r 12-11 -10-9-8-7-6-5- \ X x * * V x x mm* x *« XX x x * x x *x* Obrázek 2.13: Úroveň hladiny Hurónského jezera 1880 1900 1920 1940 1960 1980 1880 1900 1920 1940 1960 1980 (a) CV (b) procedura Obrázek 2.14: 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) 24 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 6.6 a na obrázku 2.15. Odhady vyhlazovacích parametrů podle metody křížového ověřování a pomocí automatické procedury: hCv = 0,1865 (Kopt,o,2), 0,8826 (K0pt,o,w)- Výsledné odhady regresní funkce jsou na obrázku 2.16. Je vidět, že CV metoda dává spíš pod-hlazený odhad. Na druhou stranu, odhad pomocí procedury se může zdát již přehlazený. X XX X X XX 20- , 15- 60 80 100 120 140 160 180 Obrázek 2.15: Axiální délka krystalů ledu 60 80 100 120 140 160 180 (a) CV 60 80 100 120 140 160 180 (b) procedura Obrázek 2.16: Odhadnuté regresní funkce - na ose x je vynesen čas ve vteřinách a na ose y délka krystalu v mikrometrech 25 Shrnutí Odhad regresní funkce m(x) v bodě x je tvaru i(x,h)='52Wi(x,h)Yi. i=l Wi závisí na K a h. Vychýlení (bias) a rozptyl (var) odhadu s jádrem řádu 2 jsou h2 a2 biasřn(x,/i) ~ —j32{K)m"{x), va,vfh(x,h) ~ —V(K). 2 nh Asymptotická střední kvadratická chyba jádrového odhadu regresní funkce s jádrem řádu 2 AMISE(/i) = ^rV(K) + l-h^fi2JK)V{m"). nh 4 Optimální vyhlazovací parametr vzhledem k AMISE pro k — 2 je tvaru 5 _ v2V{K){k\)2 n a 'AMISE - 4n/32(Ä-)y(m„) s řádem konvergence n-1/5. Metoda křížového ověřování pro odhad optimálního vyhlazovacího parametru 1 ™ CV(h) = -^(m_i(xi,/i) -y,)2 /icv = 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. Dodatky a cvičení 1. Pro odhady ííínw, ít^ll, íncM dokažte, že „množství" vyhlazení pomocí jádra K s vyhlazovacím parametrem h je stejné, jako „množství" vyhlazení jádrem Kg s parametrem h* = VM- ŕh{x, h, K) — m{x, h*, Kg). 2. Dokažte vztah pro střední hodnotu odhadu (2.7). 3. Dokažte vztah (2.14) pro obecné k, tj. že platí , / \ 2fe/(2fe+l) / , N \ l/(2fc+l) AMISE(/lopt,0,fe)=n-2fe/(2fe+1)((T2l/(^)) \p2{K){k\)-2V{m^)) Návod: Předpokládejme, že funkce m je spojitá a má spojité derivace až do řádu k včetně, tj. m e Ck[0,1]. Užitím Taylorova (až do řádu k) rozvoje upravíme vztah pro střední hodnotu odhadu podobně jako ve vztahu (2.8). Užijte poznámky 3.2. 4. Dokažte vztah z lemmatu 3.1. 26 (t r/ jy-\ \ -L/ V^^T-1- } -šhjt) ) ' K e S°k' k oběma částem chyby AMISE je stejný, tj. platí rovnost KL(t)dt = / *fc^í0fc(*)dí ■ -Sok \J Sok ) 6. Aplikujte metodu křížového ověřování a automatickou proceduru na simulovaná i reálná data. 27 Kapitola 3 Jádrové odhady hustoty Výstupy z výukové jednotky Student • bude znát tvar jádrových odhadů hustoty pravděpodobnosti. • bude schopen analyzovat statistické vlastnosti těchto odhadů. • se seznámí s metodami pro volbu vyhlazovacího parametru. • porozumí automatické proceduře pro simultánní volbu parametrů vyhlazování. • zvládne použití toolboxu v Matlabu a dokáže pro daný soubor dat zkonstruovat jádrový odhad hustoty a jejích derivací. 1 Motivace Hustota pravděpodobnosti je základním pojmem ve statistice [1, 3]. Odhadem hustoty rozumíme rekonstrukci hustoty z množiny naměřených dat. Tato rekonstrukce může poskytnout důležité informace o dané množině dat. Předpokládejme, že máme k dispozici nezávislé náhodné proměnné Xi,..., Xn, které mají tutéž spojitou hustotu /. Můžeme předpokládat, že neznámá hustota patří do třídy hustot, které závisejí na nějakých parametrech. Pro odhad hledané hustoty je tedy třeba odhadnout tyto parametry. Tento přístup se nazývá parametrický. My se zaměříme na neparametrický přístup, ve kterém se předpokládá pouze jistá hladkost odhadované hustoty (tj. dostatečný počet spojitých derivací). 2 Základní typy neparametrických odhadů Nejstarším neparametrickým odhadem hustoty je histogram [12, 11, 14]. Histogram zobrazuje relativní četnosti třídicích intervalů jako plochy obdélníků sestrojených nad těmito intervaly. Pak definujeme odhad hustoty četnosti f(x, h) — —(počet Xi ve stejném intervalu jako x), kde h značí šířku třídicích intervalů (obvykle se volí stejná šířka pro všechny intervaly). Nevýhody histogramu: • Histogram je citlivý na počet tříd a jejich šířku. 28 • Histogram je schodovitá funkce, ale přitom předpokládáme, že neznámá hustota je spojitá. Příklad 2.1. Mějme dán datový soubor generovaných ze směsi dvou normálních hustot N(0; 0,25) a N(2; 0,25), který má rozsah n = 100. 1 2 1 0-2)2 f (x) = 0,3^== e~fe +0,7 _ e ~ (Data jsou v tabulce 6.2.) Na obr. 3.1 je patrné, že histogram nevystihuje hustotu pravděpodobnosti dat. -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 (a) 7 intervalů o šířce h = 0,55 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 (b) 13 intervalů o šířce h = 0,29 Obrázek 3.1: Histogramy s různými počty třídicích intervalů Výše uvedené problémy lze odstranit použitím jádrových odhadů. Jádrový odhad hustoty / v bodě x G M je definovaný vztahem (Rosenblatt (1956), Parzen (1962)) 1 n Y 1 n (3.1) K G Sofc a. h je vyhlazovací parametr nebo také šířka vyhlazovacího okna. Jádrový odhad hustoty závisí na třech parametrech: jádře, které hraje roli váhové funkce, vyhlazovacím parametru, který řídí hladkost odhadu, a na řádu jádra, který odpovídá předpokládanému počtu derivací neznámé hustoty. Popíšeme konstrukci jádrového odhadu. V každém bodě Xi sestrojíme jádro Kh a odhad v bodě x je průměr n hodnot jader v tomto bodě - viz obrázek 3.2(a). Nyní uvedeme ještě vztah pro jádrový odhad v-té derivace hustoty. Budeme předpokládat, žeO < v < k — 2 a k, v jsou stejné parity. Pak 1 n Y (3.2) Konstrukce jádrového odhadu derivace je stejná jako konstrukce odhadu hustoty - obr. 3.2(b) - pouze místo jádra K G sq/c používáme jádro třídy Svk- 3 Statistické vlastnosti jádrových odhadů hustoty Stejně jako u jádrových odhadů regresní funkce lze kvalitu jádrového odhadu hustoty popsat lokálně pomocí střední kvadratické chyby. 29 0.8 (a) Hustota -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 (b) Derivace hustoty Obrázek 3.2: Konstrukce jádrového odhadu hustoty a její derivace Věta 3.1. MSE f(x,h) = E(f(x,h) - f(x)f = -((^ * f)(x) (Kh * /)2(x)) + ((Kh * /)(*) - /(x)V Důkaz. Spočítejme střední hodnotu odhadu f(x,h) 1 " /" Ä) = E- Kh(x - Xi) = ~X)= Kh(x - y)f(y) dy = * f)(x). i=0 Vychýlení (bias) pak bude mít tvar bias/(x, h) — Ef(x, h) — f(x) — (Kh * f)(x) — f(x). Dále upravíme vztah pro rozptyl 1 " 1 " 1 var f(x, h) — var — N K^(x — Xi) — — var N K^(x — Xi) — — var K^(x — X) i=l i=l = -EKl(x -X)- - (EKh(x - X))2 n n 1 1' K2h(x-y)f(y)dy--((Kh*f)(x))2 n .1 n 1 ((K2h*f)(x)-(Kh*f)2(x)). □ Důsledek. MISE = i ( /(if2 * /)(*) dx - /(A^ * /)2(x) dx) + /((Kh * /)(x) - /(x))2 dx Podobně jako u odhadu regresní funkce můžeme použít globální pohled na kvalitu odhadu, a to pomocí střední integrální kvadratické chyby (MISE) a jejího asymptotického tvaru (AMISE). Věta 3.2. Nectí funkce f má spojité derivace až do řádu k0 (tj. / e Ck°) pro 0 < k < k0, K e Sok a j(/(fe) (x)) dx < oo, dále předpokládejme h —>• 0 a nh —>• oo pro n —>• oo. PflA: pZarí MISE/(■,/!) = MISE(/i) = + ^/^(jQytfW) + 0(h2k + (nh)-1), kdeV(f^) = J(/(fe)(x))2dx. 30 Důkaz. Nejprve vypočteme střední hodnotu Ef{x,h) = (Kh*f)(x) = j Kh(x-y)f(y)dy = J (^-) 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^k\x)hkzk + o(hk) K(z) [f(x) - f'(x)hz + ■■■+ (-^f(-k\x)hkzk + o(hk)] dz = m + {~l)kf^){x)hkh{K) + o(hk). Nyní dokážeme vztah pro rozptyl. Víme, že platí vztahy 1 var/(*, h) = ^-((K2 * f)(x) - (Kh * f)2(x)) tedy Ef(x, h) = f(x) + ( VLf^(xhkpk(K) + o(hk) = f(x) + o(í), v&rf(x,h) = ± J ^K2(^-)f(y)dy-^(f(x) + o(l))2 = K2(z) f(x - hz) dz - i(/(a;) + o(l))2 =f(x)+o(l) = \( K2(z)(f(x) + o(l)) dz + -(f(x) + o(l))2 nh J n = IM. f K2(z)dz + o((nh)-1). nh . □ Důsledek. Nechí h —>• 0, nh —>• oo pro n —>• oo, paA: / je konzistentním odhadem /, ŕ/. E1/ —>• / a var / —>• 0. Vztah mezi chybami MSE, MISE a AMISE je podobný jako u odhadu regresní funkce, tedy platí MISEO) = J MSE(x, h) dx = AMISE(/i) + o(h2k + (nh)-1), kde AMISE je tvaru AMISE/(-, h) = AMISE(/i) = + J^h2k (322(K)V(f^). (3.3) V dalších částech textu budeme využívat označení jednotlivých částí chyby AMISE, která je součtem asymptotického tvaru integrálu rozptylu AIV (asymptotic integrated variance) a asymptotického tvaru integrálu druhé mocniny vychýlení AISB (asymptotic integrated squared bias): AlV(h) = TO AiSB(h) = J^h2kp22(K)V(f^), tedy AMISE (h) = AIV(/i) + AISB(/i). 31 Užitím vztahů T{K) = (V{K)k pk{K))2,{2k+1) a 5^+1 = pro K e Sok lze AMISE zapsat ve tvaru AMISEW=T(/0(^ + ^p). 0.4) Odtud je zřejmé, že vyhlazovací parametr, pro nějž AMISE nabývá minimální hodnoty, je dán vztahem . 2fe+i _ "ofc y^-i opt'°'k 2knV(f^)' tj. hoMk = 0(n-V(2fc+D). Vypočtěme hodnotu AMISE při dosazení optimálního parametru /iOpt,0 fe- AMISE( = V(^)n-2fc/(2fc+i) (2fc(^2+;(2fc+1), (3.5) tj. AMISE(/loptj0,fc) = 0(n-2fe/(2fe+D). I v tomto případě, podobně jako u odhadhu regresní funkce, platí vztah mezi AIV(h) a AISB(/i): AW(hopt^k) = 2fcAISB(/ioptj0,fc)- (3.6) Nyní uvedeme zajímavou vlastnost vyhlazovacího parametru. Poznámka 3.1. Nechť Ä" e So2- Pak optimální hodnota vyhlazovacího parametru je h5 - 02 opt'0'2 nV(f")' Počítejme derivace AMISE (3.3) Řešením rovnice d3 AMISE / dh3 — Oje V(K) č,5 02 n(32(K)V(f") nV(f") opt'0'2' tj. hopt02 také realizuje minimum d2AMISE / dh2. Obecně lze ukázat, že d2AMISE(/(-,/i)) d/i2 , _ 2fc-2 N — 0{n 2k+1 h—hoptro,k a to znamená, že pro jádra vyšších řádů je minimum AMISE plošší a tedy volba h blízká optimální hodnotě /i0pt,o,fe nevede k velkému růstu AMISE (obr. 3.3). Vztah pro optimální hodnotu vyhlazovacího parametru poskytuje informaci, že asymptoticky je /i ~ n-1/(2fe+1). Ale vztah má pouze teoretický charakter, protože optimální parametr závisí na neznámé hustotě /. Poznámka 3.2. Z předchozích úvah je zřejmé, že množina přípustných hodnot vyhlazovacích parametrů je dána vztahem H„ = [afen-Wi),&fen-Wi)], kde ak, bk jsou konstanty, 0 < ak < bk < oo. 32 Obrázek 3.3: AMISE pro jádra vyšších řádů s vyznačenými minimálními hodnotami pro jádra řádu 2,4, 6 3.1 Odhad derivace hustoty Pojednáme nyní stručně o statistických vlastnostech jádrových odhadů derivace hustoty. Připomeňme, že jádrový odhad derivace hustoty je dán vztahem (3.2), tj. 1 n Y i=l Předpokládejme nyní, že platí 0 < v < k — 2, linin^oo h — 0, linin^oo nh2v+1 — oo, / e Ck° (k < ko) a V(f^) — j(f(k\x))2dx < oo. Pak lze ukázat, že asymptotická střední kvadratická chyba AMISE /»(■, h) je tvaru AMISE/»(,/*) = + ^h2^Pl{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. Optimální hodnota vyhlazovacího parametru je dána vztahem k+1 S^(2v + l)(klf V {KM) flopt,u,k 2n(k - v)V(f(k)) ' vk PliKWy Tento vzorec umožňuje výpočet optimálního vyhlazovacího parametru pro pomocí /iOpt,0,fc a /i0pt,i,fe- Předpokládejme nejdříve, že ľ a k jsou sudá čísla. Pak hopt,,,k _ f(2v+l)k\1/i2k+1) 5vk hopt,o,k V k — v J Sok Pro vak lichá platí (3.7) iort,v,k = /(2^+l)(fc-l)\1/(2fc+1) 'iopt,i,k v - í/) / 5ifc' Speciálně pro i/ — 2, A; = 4 dostáváme velmi užitečný vztah ^opt,2,4 — ÍO1/9-^/iopt,o,4, (3.9) 004 33 pncemz KoptM*) = §(^2 - l)(7x2 - 3)> ^04 = 2,0165, 105 ^(2)(x) = ^opt,2,4(a:) = -^-(1 - x2)(5x2 - 1), ô24 = 1,3925. 4 Volba jádra Volba jádra není z asymptotického hlediska podstatná, jak je zřejmé z faktu (3.5). Je vhodné zvolit optimální jádro, které minimalizuje fukcionál T(K), neboť tato jádra jsou spojitá na M a hladkost jádra také „zdědí" odhadovaná hustota. 5 Volba vyhlazovacího parametru Volba vyhlazovacího parametru pro jádrový odhad hustoty je, stejně jako u regrese, zásadním problémem. I když tomuto problému byla a je věnována značná pozornost, doposud neexistuje univerzální přístup k řešení tohoto problému. Nejjednodušší metoda je „okometrická". Je účelné „nakreslit" několik křivek s různými vyhlazovacími parametry dříve než užijeme nějakou automatickou proceduru. Je třeba zdůraznit, že z hlediska analýzy všechny volby vyhlazovacího parametru vedou k užitečnému odhadu hustoty. Velká šířka okna charakterizuje globální strukturu hustoty a naopak malá šířka odhaluje lokální strukturu, která může nebo nemusí být přítomná v přesné hustotě. Tuto myšlenku ilustruje obrázek 3.4, na němž jsou zobrazeny odhady pro simulovaná data (n — 100) s hodnotami vyhlazovacího parametru z intervalu [0,05,1]. Jednotlivé odhady hustoty příslušející těmto hodnotám jsou znázorněny tenkými čarami. Silná křivka znázorňuje odhad s optimální hodnotou h — 0,4217. Třída těchto odhadů ukazuje široký rozsah vyhlazení od pod-hlazení až k přehlazení. 2 Obrázek 3.4: Volba vyhlazovacího parametru 5.1 Metoda referenční hustoty Nejčastěji se pro odhad neznámé veličiny V(f^) (viz rovnice (3.3)) používá parametrické třídy hustot. Jednou z možností je použít standardní normální hustotu / s rozptylem a2, tj. předpokládáme, že 1 /(*) = -/Ť^e 2'2 ■ \/2lT(7 34 V tomto případě je odhad optimálního vyhlazovacího parametru tvaru pro K e Sofc /22fc(fc!)y^\^_ ^ref = —,n,(,, sokem , (2k)\k ) Je třeba ještě odhadnout směrodatnou odchylku a. To lze dvěma způsoby: Vn — 1 z—' / n z—' (3.10) (3.11) (3.12) kde $ 1 je standardní normální kvantilová funkce a číslo A[3„/4], respektive X[„/4], je horní, respektive dolní výběrový kvartil. Je vhodné volit mm{dsD,&iQii}- Poznámka 5.1. Pokud za jádro K zvolíme Gaussovo jádro (k — 2) K(x) 2tt pak dostaneme jednoduchý vztah (3.13) Příklad 5.1. Použijme odhad vyhlazovacího parametru pro data z příkladu 2.1 metodou referenční hustoty. Pro Epanečnikovo jádro, které je řádu k — 2, se vztah (3.10) zjednoduší na tvar ísv^\1/5x ~ -1/5 /iref = ( J s02cm 1/5. Dále odhadneme směrodatnou odchylku: cťsd — 1,0325, viqr — 1,3344, tedy a — 1,0325. Po dosazení počtu prvků n — 100 a parametru S02 — 1,7188 získáme hodnotu vyhlazovacího parametru pro odhad hustoty /iref — 0,9639. Na obrázku 3.5 je vykreslen odhad hustoty s tímto parametrem. Obrázek 3.5: Odhad hustoty s /iref — 0,9639, odhad (—), původní funkce (—) 5.2 Metoda maximálního vyhlazení Princip maximálního vyhlazení (maximal smoothing) - MS (nebo přehlazení) znamená, že vybereme největší stupeň přehlazení kompatibilní s odhadovanou hustotou. Získáme tak horní 35 hranici pro odhad optimální šířky vyhlazovacího okna. Tato hodnota pak může sloužit jako počáteční aproximace pro některé z dalších metod. Princip spočívá v tom, že hledáme hustotu, pro kterou V(f^) nabývá minimální hodnoty, a tedy vztah pro /iOpt,0,fe nabývá maximální hodnoty. Věta 5.1 (Terrell 1990). Mezí všemi hustotami f s nosičem [—1,1] má hustota rozdělení Beta(k+2, k + 2) ( (2fc+3)! (í _ 2)k+l gk(X) = \ (fe+l)!22fe+3 ^ X > |X| - 1' 10 jinak, nejmenší hodnotu integrálu (/(fe) (x))2 dx. Lze ukázat, že platí L al = /-i x29k(x) dx = 2^5. 2. Pro r > 0, f (rg(k\rx))2 dx — r2k+1 J (g^ (x))2 dx pro každou hustotu, pro kterou integrál existuje. 3. Jestliže hustota g má rozptyl a2, pak hustota ^-g{^-x) má rozptyl a2. (Podrobněji např. [6].) Jestliže / je neznámá hustota s rozptylem a2 a gk je hustota rozdělení Beta(/j + 2, k + 2), pro kterou je V(g^) minimální, pak hopt,0,k < í>0fc Hodnotu Xj) = EK^ X^ 1 ' i=i j=i Odtud Kh(x -y) f (y)J (x) dxdy — E \ f(x,h)f(x)dx E f (x,h) ECV(h)=E / f(x,h)dx-2E / f(x,h)f(x)dx □ Odhad /ioptío,fc je dán vztahem/icv — argmin^g^ C V (h). Odtud plyne, že C V (h) + j f2 (x) dx je pro každé h nevychýleným odhadem MISE(/i). Protože J f2 (x) dx nezávisí na h, minimalizace EGV(h) odpovídá minimalizaci MISE. Jestliže předpokládáme, že minCV(/i) ~ min E GV(h), dostaneme dobrou aproximaci optimální hodnoty h. Obrázek 3.7: Porovnání minima MISE (červenou) a minima funkce křížového ověřování CV (modrou) pro simulovaná data Poznámka 5.2. Předpokládejme, že k — 2. Pak vychýlení odhadu může být velké, jestliže /" nabývá velkých hodnot, tj. křivost hustoty je velká. Při vyhlazování se tato objevuje ve „vrcholech", kde je vychýlení záporné, nebo v „údolích", kde je vychýlení kladné. Odhad má tendenci „vyhladit" tyto jevy, jak je patrné z obrázku 3.8. Příklad 5.4. Jádrový odhad hustoty dat z příkladu 2.1 je zobrazen na obr. 3.9. Pro rekonstrukci bylo použito Epanečnikovo jádro a vyhlazovací parametr určený metodou křížového ověřování hCY = 0,5628. Shrneme-li na závěr doposud vypočítané hodnoty vyhlazovacích parametrů pro simulovaná data z příkladu 2.1, můžeme vizuálně porovnat jednotlivé odhady - viz obrázek 3.10. Vt,o,2 = 0,5122 hREF = 0,9639 hMS = 1,0409 hCy = 0,5628 38 Obrázek 3.8: Zahlazení vrcholů a údolí při odhadu hustoty směsi normálních rozdělení 0.6r Obrázek 3.9: Odhad hustoty s hcv — 0,5628, odhad (—), původní funkce (—) 6 Automatická procedura Obdobně jako v případě regresní funkce můžeme nalézt podobnou formuli pro AMISE/(-, h), ve které budou jednotlivé parametry K, h, k separovány, což nám umožní navrhnout proceduru pro simultánní volbu těchto parametrů. Vyjdeme ze vztahu AMISE(/iopt,o,fc) = T{K) nh. ôok , hfMkV(f^)\ s^k(k\y )■ opt,0,k Ze vztahu pro hopt^,k vypočteme V(f^) V U I - oí.„r.2fe+l Ohnh2k+1 Alínnoptfl,k a tuto hodnotu dosadíme do předchozího vztahu: AMISE(/iopt,o,fc) = T{K) (2k+í)60k 2nkh, opt ,0,/c Tento vztah je základem automatické procedury, označme jej L(k) ve shodě s označením u regresní funkce. Podobně množinu vhodných řádů k označme /(fco) = {2j,j = 0,...,[^]} 39 (a) Odhad s hopt,o,2 (b) Odhad s aref o 1 2 (c) Odhad s/ims o 1 2 (d) Odhad s hcv Obrázek 3.10: Srovnání odhadů pro data z příkladu 2.1 Krok 1 Pro k G I{ko) najděte optimální jádro Koptto,k <= Sok, které je dáno tabulkou 1.2, k němu příslušný kanonický faktor Sok- Krok 2 Pro k e /(fco) a Koptto,k <= Sofc najděte optimální vyhlazovací parametr hoptto,k- Krok 3 Pro k e I(ko) vypočtěte hodnotu výběrového kriteria L(k) s využitím hodnot získaných v krocích 1 a 2. Krok 4 Vypočtěte optimální hodnotu řádu k, které minimalizuje funkcionál L(k). Krok 5 Použijte parametry z předchozích kroků ke konstrukci optimálního jádrového odhadu hustoty, tj. ' x - X, ^^opí,0,/c i—l opí,0,/c Příklad 6.1. Aplikace procedury na data z příkladu 2.1. Maximální řád jádra zvolme ko — 8, tedy množina možných řádů jader je 1(8) — {0, 2, 4, 6, 8}. Pro tyto řády spočítejme hodnoty z kroků 1-3. V toolboxu Matlabu, který je doprovodným materiálem těchto skript, je jako implicitní metoda pro odhad vyhlazovacího parametru automatickou procedurou použita iterační metoda (podrobněji např. [4]). Proto při výpočtu optimálních parametrů použijeme mezivýpočty z tohoto toolboxu. 40 k K opt,o,k Sok h L(K) 2 -|(x2-l) 1,7188 0,5314 0,0141 4 i§(x2 - l)(7x2 - 3) 2,0165 1,0734 0,0131 6 -i2|(x2-l)(33a;4-30a;2 + 5) 2,0834 1,6460 0,0125 8 J^(x2-l)(715x6-100Lr4 + 385:r2 - 35) 2,1021 2,1367 0,0126 Z tabulky vidíme, že optimální řád jádra je k — 6. Výsledný odhad je uveden na obrázku 3.11. 0.6 0.1 0.5 0.4 0.3 0.2 4 Obrázek 3.11: Simulovaná data (x) s jádrovým odhadem hustoty při použití procedury (—) a původní funkcí (—) Při bližším pohledu na odhadnutou hustotu je patrný vliv použití optimálního jádra vyššího řádu. Jádra vyšších řádů mohou nabývat záporných hodnot a tím ovlivnit výslednou odhadnutou funkci - viz obrázek 3.12. V takovém případě je vhodné použít jinou metodu pro nalezení vyhlazovacího parametru, případně použít jiné jádro. Lze doporučit jádra třídy S02, např. kvar-tické jádro: K{x) — 1|(1 — x2)2I[_ltl](x), nebo jádro triweight: K{x) — ||(1 — x2)3I[_ltl](x). Datový soubor obsahuje morfologická měření padesáti exemplářů od obojího pohlaví a obou barevných forem (oranžová a modrá) krabů rodu Leptograpsus.1 Pro odhad hustoty nám postačuje jeden druh měření, vybrali jsme délku podél středové osy krunýře, která byla měřena v milimetrech. Data jsou uvedena v tabulce 6.7. Užitím výše uvedených metod pro odhad vyhlazovacího parametru jsme (při použití Epaneč-nikova jádra) dostali následující hodnoty: U automatické procedury je v toolboxu implicitně nastavena iterační metoda pro odhad vyhlazovacího parametru, proto jsme tuto volbu ponechali i zde, af má čtenář možnost porovnání při vlastních výpočtech. Při použití procedury vyjde vyhlazovací parametr roven hpľoc — 31,7329 s optimálním jádrem Kopt,o,8- Výsledné odhady hustoty na jsou zachyceny na obrázku 3.13. Celý datový soubor je dostupný v programu R. 7 Aplikace na reálná data ^ref — 5,7856, hMs = 6,2480, hCv = 8,1317. 41 (a) Jádro A"optio,6 (b) Odhad hustoty Obrázek 3.12: Jádro třídy Soe a k němu příslušný odhad hustoty při použití procedury (—) a původní funkcí (—) 10 15 20 25 30 35 40 45 50 55 10 15 20 25 30 35 40 45 50 55 (c) CV (d) procedura Obrázek 3.13: Grafy odhadnutých hustot pro délku krunýře 42 Shrnutí Odhad hustoty pravděpodobnosti / v bodě x je tvaru ' x — Xi nh i=l Asymptotická střední kvadratická chyba jádrového odhadu hustoty pravděpodobnosti s jádrem řádu k je součtem asymptotického tvaru rozptylu (AIV) a druhé mocniny vychýlení (AISB) AMISE(/0 = + J^h2kfá(K)V(fW). AIV " AISB ' Optimální vyhlazovací parametr vzhledem k AMISE pro odhad hustoty pravděpodobnosti s jádrem řádu k je tvaru . 2fe+i °ofc y^-i °Pt'°'fe 2knV(f(k))1 tj. hoMk = 0(n-V(2fc+D). Metody pro odhad optimální hodnoty vyhlazovacího parametru /i • metoda referenční hustoty /22fe/j!3v/^\5^(. __i_ %EF = (vT^) Sok(7n 2k+1' • metoda maximálního vyhlazení hMS = *n-W°+»bk, • metoda křížového ověřování r 2 ™ - hCv = arg min CV(/i) = / /2(x, /i) dx--V" f-l(Xl, h) i=i Automatická procedura pro simultánní volbu optimálního jádra, vyhlazovacího parametru a řádu jádra je dostupná v toolboxu Matlabu. Dodatky a cvičení 1. Dokažte vztah (3.4) pro tvar chyby AMISE. 2. Dokažte vztah (3.13). 3. Spočítejte (3.10). 4. Spočítejte hyis Pr° • Epanečnikovo jádro K(x) — |(1 — x2), je-li = (i/zy ~ „ *0 _ 5/7 čf(iľ) - (1/7)2 kvartické jádro K{x) = ±§(1 - x2)2, je-li = = 35. 43 5. Aplikujte metody pro odhad vyhlazovacího parametru a automatickou proceduru na simulovaná i reálná data. 44 Kapitola 4 Jádrové odhady distribuční funkce Výstupy z výukové jednotky Student • bude znát základní typy jádrových odhadů distribuční funkce a jejich statistické vlastnosti. • získá přehled o metodách pro volbu vyhlazovacího parametru. • bude schopen navrhnout a implementovat proceduru pro zpracování reálných dat. • se naučí používat příslušný toolbox v Matlabu a dokáže zkonstruovat jádrový odhad distribuční funkce pro daná reálná data. 1 Motivace Distribuční funkce popisuje rozložení pravděpodobnosti náhodné veličiny (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. Nejužívanějším neparametrickým odhadem distribuční funkce F 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ů Nechť Xi,..., Xn jsou nezávislé náhodné proměnné, které mají tutéž spojitou hustotu / a distribuční funkci F. Nejjednodušší neparametrický dohad distribuční funkce F je empirická distribuční funkce Fn definovaná v bodě x vztahem 1 " F„(x) = - ^/(_00í2.](A"í). i=l Tento odhad má sice dobré statistické vlastnosti, ale je to schodovitá funkce (viz obr. 4.1, a proto se budeme zabývat postupy, které umožní zkonstruovat „hladký" odhad distribuční funkce F. 45 Příklad 2.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 3 9a:2 1 O-2)2 fix) = 0,5—7= e" V +0,5^ e-1-^ 2V27T (Data jsou v tabulce 6.3.) Z obrázku 4.1 je patrné, že schodovitá funkce nevystihuje plně charakter distribuční funkce. Obrázek 4.1: Empirická distribuční funkce (červeně) a skutečná distribuční funkce (modře) pro data z příkladu 2.1 Nejznámější postup spočívá v integraci jádrového odhadu hustoty, t.j. 't-X. Užijeme-li substituce y — (t — Xi)/h, dostaneme /X -t n r-X fit, /i) dí = — £ / K -co n" „■_-, J — co dt. F(x,h) = ±J2Í h Kiy)dy^-J2wÍX X> To znamená, že odhad F v bodě x e M je definován takto ' x - X, Fix,h) = WW i=l W{x) = / K(t) dt. (4.1) Zde předpokládáme, že K e So 2/ K(x) > 0 pro x e [—1,1]. Níže jsou shrnuty základní vlastnosti funkce W: 1. W(x) — 0 pro x e (—00, —1] a W(x) — 1 pro x e [1, 00), 2- /-1 W2{x) dx < J\ W{x) dx = 1, 3. /_\ W{x)K{x) dx = i, 4. J\ xWix)Kix) dx=\(l- J\ W2ix) dx). 46 -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 Obrázek 4.2: Epanečnikovo jádro K (vlevo) a k němu příslušná funkce W (vpravo) Příklad 2.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. Pro data z příkladu 2.1 je jádrový odhad distribuční funkce zachycen na obrázku 4.3. Obrázek 4.3: Jádrový odhad distribuční funkce s parametrem h — 1,5 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)f = (EF(x, h) - F{x)f + E(F(x, h)f - (EF(x, h)f . 47 Spočítejme nejdříve hodnotu EF(x, h) v bodě ieR: EF(x,h) = Jw(^^jf(y)dy W(t)f(x-ht)dt + h / W(t)f(x-ht)dt. -oo 7l Podrobnější výpočet je uveden v práci [10]. Zde uvedeme jen vybrané kroky výpočtu. Označme první integrál E a druhý I2 ■ Spočítejme první z nich. S využitím vlastností funkce W dostaneme E = -F (x - h) + J K (t) F (x - ht) dt. Dále použijeme Tayloruv rozvoj (4.2) h2t2 F (x - ht) = F (x) - htF'(x) + —F" [x) + o(h2), tedy Je zřejmé, že pro integrál I2 platí E = -F(x -h) + F (x) + ^F"\x)h2 p2{K) + o(h2). h = F(x-h). (4.3) Vychýlení odhadu je tedy tvaru Has F {x, h) = ^F"(x)h2fi2(K) + o(h2). Poznámka 3.1. Vztahy (4.2) a (4.3) dávají zajímavý vztah pro vychýlení EF(x,h) - F(x) = J K(t)F(x- ht)dt- F(x). Výpočet tvaru pro rozptyl je komplikovanější, a proto ho nebudeme uvádět. Platí varF(x,h) = -F(x)(l-F(x))--hf(x)(l- f W2(t) dt] + o (— Důkaz tvaru rozptylu lze najít např. v [4]. Výše uvedené výsledky můžeme nyní zformulovat v následující větě: Věta 3.1. Nectí F e C2, h —>• 0, nh —>• 00 pro n —>• 00. PflA: MSEFfo/i) = 3F(j:)(1-F(i))-3/i/(i) ( 1- / W2(t) dt ) + -{F" {x))2hA fi22{K)+o ( 3 +/ť (4.4) n n \ / 4 \ n Globální pohled na kvalitu odhadu lze získat prostřednictvím střední integrální kvadratické chyby (MISE). Věta 3.2. Nectí F e C2, V(F") = j(F"(x))2 dx < 00, K e S02, lim„^oo h = 0 a lim^oo nh = 00. Pak MISEF(-,/i) = - /f(x)(l-f(x))d:r-c1- + C2/i4 + o( - +/i4 ) , (4.5) n 7 n \n J kde ci = 1 - y i w2(t) dt, C2 = i/3|(i<:)y(F"). 48 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: 1 ľ h AMISEF(-,h — AMISE(/i) — - / F(x)(í - F(x)) dx - Cl- + c2h4 . (4.6) n J n AISB AIV Nyní už lze standardními metodami matematické analýzy nalézt takovou hodnotu h, pro kterou AMISE(h) nabývá minimální hodnoty. Je snadné ukázat, že hm2 = n-1^ (JLJ = 0(n-i/3) (4.7) a pak AMISE(/lopt,o,2) = - J F{x){\ - F{x)) dx - — (^) n~4/3. (4.8) Poznámka 3.2. Optimální hodnota vyhlazovacího parametru pro odhad distribuční funkce je řádu n-1/3, zatímco pro odhad hustoty s jádrem K e S02 je vyhlazovací parametr řádu n-1/5. 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) — i|(l - x2)2I[_1 ^(x), • triweight K(x) = - x2)3/[_ia] (x). 5 Volba vyhlazovacího parametru 5.1 Metody křížového ověřování Metody křížového ověřování patří k nejužívanějším metodám pro volbu vyhlazovacího parametru. Zde uvedeme pouze metodu navrženou A. Bowmanem (1998). Funkce křížového ověřování je v tomto případě tvaru CVW = Í Ž / (k-oo,x] (Xi) - F-i(*, h)) 2 dx, kde F-i(x, h) je jádrový odhad distribuční funkce s vynecháním bodu Xi. Pak /icv — arg min CV(/i), h£Hn přičemž Hn — [an-1/3, brT1^] pro vhodná 0 < a < b < 00. 49 5.2 Princip maximálního vyhlazení Myšlenka této metody je stejná jako pro odhad hustoty. Užijeme-li faktu, že J(F"(x))2dx = J(f(x))2dx, můžeme aplikovat Terrelovu větu 5.1 pro k — 1. V tomto případě je 15 #1 (z) = Tg (1 - z2)2i"[-i,i] (x), a tedy kde (Ti — J x2gi(x) dx — j, V(g[) — Odtud plyne, že a je odhadem a (viz rovnice (3.11) a (3.12)). Hodnota Iims 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 h t je nejmenší vzdálenost mezi po sobě jdoucími body Xi, í — 1,..., n. Příklad 5.1. Pro data z příkladu 2.1 zvolme Epanečnikovo jádro. Pak hodnoty potřebné pro odhad vyhlazovacího parametru metodou maximálního vyhlazení jsou následující: i ŕ n =100, £ = 1,3426, /32(K) = -, a = 1 - / W2(x) dx = 0,2571. 5 J-i Pak platí fíMS — 1; 1037 a na obrázku 4.4 je zobrazen odhad distribuční funkce. -2 -1 0 1 2 3 4 5 Obrázek 4.4: Odhad distribuční funkce s /ims — 1,1037, odhad (—), původní funkce (—) 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))2dx = - J f'(x)f(x)dx. 50 Tudíž se budeme dále zabývat odhadem funkcionálu V>i= / 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.7) použijeme k odhadu druhé derivace s jádrem — Kopt,2,4 € S24. Pak n n n / Y _ Y a = n-1 y, /"(*>h) = -~2^3 E E Ki2)' i=l i=l kde podle vztahu (3.9) je a tedy «opt,2,4 — 1U T—«opt,0,4, ^04 1 "2, 4' Shrnutím předchozích úvah dostaneme proceduru pro odhad distribuční funkce F: Krok 1 Najděte optimální vyhlazovací parametr hopt,o,4 Pro odhad hustoty s optimálním jádrem Kopt,0,4 <= i s využitím hodnot /iOpt,0,4 a hoptt2,4 získaných v krocích 1 a 2. Krok 4 Vyčíslete optimální hodnotu vyhlazovacího parametru / \1/3 API = n-v3 ' cl » Krok 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 5.2. S použitím funkce toolboxu zjistíme, že pro data z příkladu 2.1 je vyhlazovací parametr určený plug-in metodou roven hpi — 0,5717. Na obrázku 4.5 je odhad distribuční funkce společně se skutečnou distribuční funkcí. 6 Aplikace na reálná data 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 dihydrofolát reduktázy byla měřena pomocí asociační konstanty. Data jsou v tabulce 6.8 a jsou dostupná na osobních stránkách Dennise D. Boose1, 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.6 jsou uvedeny odhady distribuční funkce s těmito parametry a také je zde pro srovnání uvedena empirická distribuční funkce. 1 http://www4.stat.ncsu.edu/~boos/var.select/pyrimidine.html 51 Obrázek 4.5: Odhad distribuční funkce s hpi — 0,5717, odhad (—), původ] Shrnutí Odhad distribuční funkce F(x) v bodě x je tvaru n*) = \ E w (^ir1)' w{x) = /lm áf- Asymptotická střední kvadratická chyba jádrového odhadu distribuční funkce je součtem asymptotického tvaru rozptylu (AIV) a druhé mocniny vychýlení (AISB) 1 ľ h AMISE(/i) = - / F(x)(l - F(x))dx-Cl-+c2h4, n J n ^-^^ "-v-' AISB AIV kde ci = 1 - J i W2(t) dt, c2 = -^1{K)V{F"). Optimální vyhlazovací parametr vzhledem k AMISE pro odhad distribuční funkce je tvaru 1/3 t.j. /iopt.o.2 = C^n-1/3). Metody pro odhad optimální hodnoty vyhlazovacího parametru h • metoda křížového ověřování hcv — argmin/jg^ CV(/i) CVW=^E / (l(-oo,X](Xi) - F-i(X>h))2 i—l • metoda maximálního vyhlazení -"/3 (lák)1,3 ^ • plug-in metoda \-a/32(K)J Dodatky a cvičení 1. Odvoďte tvar funkce pro kvartické jádro K(x) — y|(l — x2)2/[_i!i](x). 2. Dokažte vlastnosti 2,3 a 4 funkce W. 3. Dokažte vztahy (4.7) a (4.8). 4. Odvoďte tvar vyhlazovacího parametru podle metody maximálního vyhlazení pro Epa-nečnikovo jádro. 53 5. Odvoďte tvar vyhlazovacího parametru podle plug-in metody pro Epanečnikovo a pro kvartické jádro. 6. Aplikujte metodu maximálního vyhlazení a plug-in metodu na simulovaná i reálná data. 54 Kapitola 5 Jádrové odhady dvourozměrných hustot Výstupy z výukové jednotky Student • bude znát součinová a sférická dvourozměrná jádra pro odhady dvourozměrných hustot. • porozumí principu vyhlazování dvourozměrných hustot. • pochopí 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. 1 Motivace V této kapitole se budeme zabývat rozšířením jádrových odhadů pro jednorozměrné hustoty na odhad vícerozměrných hustot. 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) v kapitole 3, a zaměříme se zejména na odhad dvourozměrné hustoty. 3,-.---g-.-1 3 g I-,-,-,-1 _gl-,-,-1— 2 1 0 1 2-2-1 0 1 Obrázek 5.1: Náhodný výběr a jeho jádrový odhad 55 Poznámka 1.1. Jádrové odhady dvourozměrných hustot se obvykle znázorňují pomocí vrstevnic, které umožňují snazší náhled na odhadnutou funkci. 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 2.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) f(x v) = lJLe-l(*2+(y+tr)+l±_ e-H-2+fe-2)2)+IAe-!(-2+fe-4)2) Jy 3 2tt 3 2tt 3 2tt (Data jsou v tabulce 6.4.) Z obrázku 5.2 je patrné, že histogram nepostihuje charakteristické rysy hustoty pravděpodobnosti dat. (a) 7 x 7 obdélníků (b) 5 x 10 obdélníků Obrázek 5.2: Histogramy s různými počty třídicích obdélníků 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 i n i n f(x, y; h) = - Ka(x -Xi,y- Yt), = — ]T ^(h^ -XllV- y,)T) (5.1) 71 i=i n' ' i=i přičemž h je matice vyhlazovacích parametrů a K je dvourozměrné jádro. Jádro K je dvourozměrná funkce, kterou můžeme získat pomocí jednorozměrného symetrického jádra K\ (K\ e £02)- Existují dva typy těchto jader: • součinové jádro Kp(x, y) — K\(x) ■ K\(y), • sféricky symetrické jádro Ks(x, y) — ck~1Kí (\Ac2 + y2), ck — J J K1 (\Ac2 + y2) dx dy. Příklad2.2. Epanečnikovo jádro, které je v jednorozměrném případě tvaru K(x) — x2)I[_i ^ (x), má následující dvourozměrné varianty Kp(x1,x2) = -^(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 součinového Epanečnikova jádra v závislosti na třídě matic. Tabulka 5.1: Třídy vyhlazovacích matic S V T Budeme se zabývat jádrovými odhady s diagonální vyhlazovací maticí. Jádrový odhad s maticí třídy S dává ve všech směrech stejnou míru vyhlazení, což neponechává příliš mnoho prostoru pro zachycení variabilty dat. Na druhou stranu při použití matice třídy T je potřeba odhadnout větší počet parametrů, což znamená vyšší výpočetní náročnost. Konstrukce jádrového odhadu je analogická konstrukci jednorozměrného odhadu. Tedy v každém bodě [Xi, 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 f(x, y; H) = - * /)(*, V) ~ {Ku * f?{x, y)) + ((Ku * f)(x, y) - f(x, y)f, var bias 57 nebo globálně pomocí střední integrální kvadratické chyby MISE/(x,y;H) = J J MSE f (x, y; H) dx dy. Poznámka 3.1. Podobně jako v jednorozměrném případě se definuje konvoluce funkcí dvou proměnných. Nechť jsou dány funkce / a g, pro které platí f f f2 (x, y) dx dy < oo a f f g2 (x, y) dx dy < oo. Konvoluci f * g definujeme vztahem (f*g)(x,y)= // f{t,u)g(x -t,y -u)dtdu. 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ěta 3.1. Předpokládejme, že funkce f, jádro K a matice vyhlazovacích parametrů H — diag(/ii, h2)2 splňují následující předpoklady. (i) Nectí H — H„ je posloupnost matic takových, že (n \ H |) ~1 a prvky matice H konvergují k nule pro n —> oo. (íí) Dále nectí všechny druhé parciální derivace funkce / jsou ohraničené, spojité a integrovatelné se čtvercem. (ííí) Jádro K splňuje xK(x,y)dxdy — J J yK(x, y) dx dy — 0 x2K(x, y) dxdy — jjy2 K (x, y) dxdy — Í32(K). Pak platí kde MISE(H) = AMISE(H) + o(hj + h22) + o((/ii/i2n)"1), AMISE(H) = AMISE/(■, H) ^ +^2(K)(hiV(fxx) + 2h2h2V(fxfy) + h2V(fyy)), (5.2) přičemž označení je ve shodě s předchozími kapitolami, tj. V(g) — JJ g2(x, y) dx dy. 1' hx 0 2Užíváme zde zkráceného zápisu: diag(fti, h^) ' 0 h2 58 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\, h^, pro které AMISE(/ii, h^) nabývá minimální hodnoty, jsou dány vztahy: h = (_V^(fyy)V(K) \1/6 'OPt \nfô(K)V^(fXX)[V(fXfy) + y/V(fxx)V(fyy)\ J ' ( V3/4(fxx)V(K) h2'°Pt = \nPl{K)VV^fyy)[V{fxfy) + ^V(fxx)V(fyy)\ ) ■ (5'4) Z těchto vztahů plyne, že množina přípustných vyhlazovacích parametrů je tvaru ain-1/6 < hi < birT1^, a^n-1/6 < hi < b2"nTxl^ pro vhodné konstanty 0 < a\ < b\ < oo, 0 < a,2 < &2 < oo. 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 4.1. V literatuře se také využívá Gaussovo jádro K(x, y) — (2tt)^1 e-^' +v ^2, které se zdá být výhodnějším při studiu asymptotických vlastností jádrového odhadu. Na druhou stranu má nevýhodu, že jeho nosičem je celá reálná osa, což způsobuje „nedokonalost" při odhadech hustot s omezeným definičním oborem. 5 Volba vyhlazovacího parametru 5.1 Metoda referenční hustoty Předpokládejme, že náhodný výběr ([Xi, Y\],..., [Xn,Yn]) pochází z normálního rozdělení s hustotou f(x,y) = --e "i -i Z7T(7l(72 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 tyto odhady vyhlazovacích parametrů /ij.ref = CTjn-1/6, 2 = 1,2. Tento vztah je také znám jako Scottovo pravidlo ([11]). Příklad 5.1. Pro simulovaná data z příkladu 2.1 vychází matice vyhlazovacích parametrů podle metody referenční hustoty s využitím součinového Epanečnikova jádra takto: /0,3595 0 Href — \ 0 0,9972^ Na obrázku 5.5 je vykreslen odhad hustoty s touto maticí a porovnání odhadu se skutečnou hustotou. Poznámka 5.1. V toolboxu, který doplňuje tato skripta, se uvádí druhá mocnina matice vyhlazovacích parametrů, tedy H2. Z obrázku 5.5 je patrné, že jádrový odhad je podhlazený, zejména ve směru osy x. Je to způsobeno odhadem směrodatné odchylky, která je základem metody referenční hustoty. 59 -2 0 2 -2 0 2 (a) Skutečná hustota (b) Data (c) Odhad s Href (d) Porovnání odhadu (—) se skutečnou hustotou (—) Obrázek 5.5: Jádrový odhad dvourozměrné hustoty - referenční hustota 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, li] 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). n ■ i=l kde 1 " f-i(Xi, li, H) = —- Ka(Xi - Xj,Yí - Yj) i=i Někdy se metoda C V nazývá nevychýlená metoda křížového ověřování (unbiased cross-validation), důvodem je jednoduchý vztah mezi CV a MISE, který uvádí následující věta. Věta 5.1. Funkce CV je nevychýleným odhadem MISE, tj. platí ECV(H) = MISE(H) - J J f(x,y)dxdy. Důkaz. Vypočtěme střední hodnotu CV: £CV(H) = £ (if2{x,y1H)áxáy--YJEf_l{XllYlllÍ) •* •* n i=l = E J J /2(x,y,H) áxdy-2EKn{X1-x2,Y1-Y2) = E f2(x,y,H)dxdy - 2 //// Ku(x - u, y - v)f(x, y)f(u, v) dxdy du dv 60 a úpravou MISE MISE(H) = e J J (f(x,y,H)-f(x,y)) dxdy = eJJ f2(x,y,H)dxdy-2E J f {x,y,K) f {x,y) dxdy + J J f {x,y) dxdy = E íl f2 (x, y,H) dxdy - 2 jjjj kH(x - u, y - v) f (x, y) f (u, v) dxdy dudv f2(x, y) dx dy. Porovnaním upravených výrazů dostaneme tvrzení. □ Optimální matici vyhlazovacích parametrů vzhledem k metodě C V označíme H^v/ tj. Hcv — arg min CV(H). Hex> Příklad 5.2. Použijeme-li součinové Epanečnikovo jádro, pak pro simulovaná data z příkladu 2.1 dostanem matici vyhlazovacích parametrů určenou podle metody křížového ověřování v následujícím tvaru: '1,2055 0 0 0,3783, Na obrázku 5.6 je vykreslen odhad hustoty s touto maticí. H cv -202 (a) Skutečná hustota X Xw X *x «*„*»«-* •3<« * x x x * x -2 0 2 (b) Data (c) Odhad s hc (d) Porovnání odhadu (—) se skutečnou hustotou (—) Obrázek 5.6: Jádrový odhad dvourozměrné hustoty - metoda křížového ověřování Odhad hustoty na obrázku 5.6 se jeví podhlazený, zejména ve směru osy y. Metoda křížového ověřování nedává dobré výsledky pro odhad hustoty vícerozměrných dat, částečně je to způsobeno problémy při minimalizaci funkce křížového ověřování. 6 Aplikace na reálná data Tento 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. Data jsou shrnuta v tabulkách 6.9 a 6.10. 61 Obrázek 5.7: 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) Matice vyhlazovacích parametrů určené podle metody referenční hustoty a metody křížového ověřování jsou následující: í16,45 0 \ /42,39 0 \ íIref — tlcv — \ 0 38,94y y 0 29,88y Na obrázku 5.7 jsou znázorněna data a vrstevnice jádrového odhadu s Epanečnikovým součinovým jádrem. 62 Shrnutí Odhad dvourozměrné hustoty pravděpodobnosti f(x, y) v bodě [x, y] je tvaru 1 n n' ' i=i Dva typy jader: • součinové jádro: Kp(x,y) — Ki(x) ■ K\(y), sféricky symetrické jádro: Ks{x, y) — ck xK\{^Jx2 + y2), Ck — jj K\{^J x2 + y2) dx dy. Asymptotická střední integrální kvadratická chyba dvourozměrného jádrového odhadu AMISE(H) = + -^l{K)(h\V{fxx) + 2hjhlV(fxfy) + h22V{fyy)). Optimální vyhlazovací parametry vzhledem k AMISE \ 1/6 _V3/\fyy)V(K) 'OPt \nft(K)W*(fxx) [V(fXfy) + y/V(fxx)V(fyy)\ h\ opt — \np2(K)v-^(;xx)[V(;x;y) + ^v(Ux)v(hv)\ ) 1/6 _ /_V3/\fxx)V(K)_ l2'°Pt " \n^(K)V^(fyy)[V(fxfy) + ^V{fxx)V{fyy)\ Metody pro odhad optimálních hodnot matice vyhlazovacích parametrů H — diag(/ii, h2) • metoda referenční hustoty • metoda křížového ověřování HCV = argminCV(H), CV(H) = í í (f(x, y, H))2 dxdy - - V f^Xi, Yit H). i=i Dodatky a cvičení 1. Určete součinové a sféricky symetrické dvourozměrné jádro odvozené z kvartického jádra A» = if(l-x2)2. 2. Odvoďte vztahy (5.3) a (5.4) pro optimální vyhlazovací parametry. 3. Aplikujte metodu referenční hustoty a metodu křížového ověřování na simulovaná i reálná data. 63 Kapitola 6 Přílohy 1 Symbolika O, o Symboly O a o se často používají pro vyjádření chyb matematických výrazů Nechť / je funkce definovaná v okolí bodu a. Symbol f (x) — O (g (x)) pro x —>• a znaci, ze Podobně symbol znaci, ze nm sup —-—r- < oo. x^a g{x) f (x) = o(g(x)) pro x \m\ lim 0. 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 = 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(<7(a;)) + 0(ff(x oíffO)) + o(g(x 0(g(x)) + o(g(x 0{g{x)) ■ 0{h{x o(g(x)) ■ o(h(x O (g (x)) ■ o(h(x o(g(x = 0{g{x)) = o{g{x)) = 0(í/(x)) = 0(ff(x) • h(x)) = o(ffO) • h(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) — O (g (x)). Opačně to ovšem neplatí. 64 Datové soubory Tabulka 6.1: Hodnoty simulovaných dat z příkladu 1.1 - regrese X y X y x y x y x y 0,01 0,3002 0,02 0,9792 0,03 -1,0372 0,04 0,5519 0,05 0,3070 0,06 -0,4816 0,07 -0,0225 0,08 0,3848 0,09 2,0187 0,10 1,6268 0,11 -0,4240 0,12 1,7734 0,13 0,6198 0,14 0,2228 0,15 0,6049 0,16 0,1343 0,17 0,1602 0,18 0,9489 0,19 0,8871 0,20 0,8664 0,21 0,4661 0,22 -0,5034 0,23 0,4270 0,24 0,8499 0,25 0,2444 0,26 0,4820 0,27 0,2926 0,28 -0,2577 0,29 0,0068 0,30 -0,5664 0,31 0,2407 0,32 -0,8052 0,33 -0,7914 0,34 -0,6835 0,35 -1,7689 0,36 0,4086 0,37 -0,1572 0,38 -0,7018 0,39 0,3614 0,40 -1,1739 0,41 -0,3584 0,42 -0,4120 0,43 -0,1105 0,44 -0,0875 0,45 -0,6454 0,46 -0,1926 0,47 -0,2206 0,48 0,2188 0,49 0,4979 0,50 0,5546 0,51 -0,3811 0,52 0,1413 0,53 -0,4521 0,54 -0,3496 0,55 0,2547 0,56 1,0735 0,57 -0,0313 0,58 0,5820 0,59 0,3219 0,60 1,0265 0,61 -0,0495 0,62 0,5318 0,63 0,8049 0,64 1,0842 0,65 1,3028 0,66 0,5615 0,67 -0,2485 0,68 0,0955 0,69 -0,1043 0,70 1,5522 0,71 0,0104 0,72 0,6246 0,73 0,0784 0,74 0,5351 0,75 -0,3824 0,76 -0,7979 0,77 -0,9098 0,78 -0,0599 0,79 -0,5005 0,80 -0,6184 0,81 0,0816 0,82 -0,5874 0,83 -0,7349 0,84 -0,1342 0,85 -1,4160 0,86 -0,7407 0,87 -0,7340 0,88 -1,3214 0,89 -1,1229 0,90 -1,8261 0,91 -1,8089 0,92 -1,1517 0,93 -0,7882 0,94 0,2239 0,95 -1,2950 0,96 -0,7328 0,97 -0,7041 0,98 -1,4370 0,99 -0,4688 1,00 -0,8973 65 Tabulka 6.2: Hodnoty simulovaných dat z příkladu 2.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 Tabulka 6.3: Hodnoty simulovaných dat z příkladu 2.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 66 Tabulka 6.4: Hodnoty simulovaných dat z příkladu 2.1 - dvourozměrná hustota [x;v] [-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] 67 Tabulka 6.5: Hodnoty reálných dat z kapitolky 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. 68 Tabulka 6.6: Hodnoty reálných dat z kapitolky 7 - krystaly ledu x V x V x V x V x V 50 19 60 20 60 21 70 17 70 22 80 25 80 28 90 21 90 25 90 31 95 25 100 29 100 30 100 33 105 32 105 35 110 28 110 30 110 30 115 30 115 31 115 36 120 25 120 28 120 36 125 28 130 31 130 32 135 25 135 34 140 26 140 33 145 31 150 33 150 36 155 33 155 41 160 30 160 37 160 40 165 32 170 35 180 38 Tabulka 6.7: Hodnoty reálných dat z kapitolky 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 69 Tabulka 6.8: Hodnoty reálných dat z kapitolky 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 70 Tabulka 6.9: Hodnoty reálných dat z kapitolky [x;y] [x;y] [x;y] 6 - koncentrace lipidů - část 1. [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] 71 Tabulka 6.10: Hodnoty reálných dat z kapitolky 6 - koncentrace lipidů - část 2. [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] 72 Literatura [1] Anděl, J.: Základy matematické statistiky. Matfyzpress, Praha (2005) ISBN 80-86732-40-1 [2] Collomb, G.: Estimation non paramétrique de la regression par la méthode du noyau. Disertační práce, Univerzita Paula Sabatiera, Toulouse (1976) [3] 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 [4] 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 [5] Horová, I., Vieu, P., Zelinka, J.: Optimal Choice of Nonparametric Estimates of a Density and of its Derivatives. Statistics & Decisions 20, 355-378 (2002) [6] Horová, I., Zelinka, J.: Contribution to the bandwidth choice for kernel density estimates, Comput. Stat. 22, 31-47 (2007) [7] 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) [8] Koláček, J.: Jádrové odhady regresní funkce. Disertační práce, Masarykova univerzita, Brno (2005) http://is.muni.cz/th/19999/prif_d/ [9] Müller, H.-G.: Smooth optimum kernel estimators near endpoints. Biometrika 78, 521-530 (1991) [10] Nováková, J.: Jádrové odhady distribuční funkce. Diplomová práce, Masarykova univerzita, Brno (2009)http://is.muni.cz/th/151323/prif_m/ [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] Terrell, G.R.: The maximal smoothing principle in density estimation. J. Am. Stat. Assoc. 85, 470-477 (1990) [14] Wand, M.P, Jones, M.C.: Kernel Smoothing. Chapman and Hall, London (1995) ISBN 0-412-55270-1 73