Ú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 šíř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/vyvi jeny-software/ 27 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 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 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.............................. 19 5.1 Metoda křížového ověřování........................... 19 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ů m(-) regresní funkce /(■) hustota pravděpodobnosti F(-) distribuční funkce h odhad vyhlazovacího parametru h a odhad směrodatné odchylky a m odhad regresní funkce m / odhad hustoty / F odhad distribuční funkce F značí integrál , pokud není uvedeno jinak jádrová funkce (jádro) V(K) = jK2(x)dx f3k(K) = /xkK(x) dx konvoluce funkcí / a g, (/ * g)(x) = J f (ť) g (x — t) dt integrál z jádra, W{x) = f*^ K (t) dt prostor funkcí, které mají spojité derivace až do řádu k včetně na intervalu [0,1] NW Nadarayovy-Watsonovy odhady PC Priestleyovy-Chaovy odhady LL lokálně lineární odhady GM Gasserovy-Múllerovy odhady / K(-) V(K) h{K) f*9 W(-) C*[0,1] 4 MSE 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. Necht v, k jsou nezáporná celá čísla, 0 < v < k, necht 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 G [-1,1], L > 0, 2. nosič(íí) = [—1,1], tj. K = 0 vně intervalu [—1,1], 3. K splňuje momentové podmínky: a J-i xkK(x) dx =jí 0, tuto hodnotu označíme fik{K). 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 J[-i,i] (x) je indikátorová funkce intervalu [—1,1], tj. nosti. (1.1) Nyní uvedeme dva typy jader - jádra s minimálním rozptylem a optimální jádra. 6 Tabulka 1.1: Jádra 5, 02 K (x) = !l[_M](x) obdélníkové jádro K(x) = Z(l-x2)I[-1,1](a Epanečnikovo jádro " -1 0 1 K{x) = -f x(l - x2)^ 0.5 A -0.5 if(a;) = if(l-a;2)2J[_1,1](a;) kvartické jádro K(x) = -f(l-3x2)IhlA](x) 1.1 Jádra s minimálním rozptylem Předpokládejme, že K e S^fc, 0 < < k—2,vak jsou obě sudá nebo lichá1. Uvažujme funkcionál V{K) = J_1 K2{x) dx a zabývejme se problémem najít takové jádro K G Svk, pro které tento funkcionál nabývá minimální hodnoty, tj. řešíme variační úlohu min V{K) za předpokladu K G 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é a 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 [3]. Příklad 1.2. Jádra s minimálním rozptylem: »502: K{x) 5i3: K(x) »524: K(x) = -f(l-3x2)/ 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. 1Obvykle se používá pojem parita, tedy vak mají stejnou paritu. 7 Tabulka 1.2: Optimální jádra pro v = 0,1,2 k Sok Kopt,0,k 2 1,7188 -i(*2-i) 4 2,0165 §(a;2-l)(7a;2-3) 6 2,0834 -±§§(a;2-l)(33a;4-30a;2 + 5) 8 2,1021 ■§§q(x2 - l)(715a;6 - 1001a;4 + 385a;2 - 35) v = 1 k Sik Kopt,i,k = 3 1,4204 f x(x2 - 1) 5 1,7656 -±§x(x2 -l)(9x2 -5) 7 1,8931 ^x(x2 - l)(143a;4 - 154a;2 + 35) ^ = 2 Jt &2k 4 1,3925 -1Sž(x2 -l)(5x2 -1) 6 1,6964 f^(a;2-l)(77a;4-58a;2 + 5) 8 1,8269 ~ÍS(x2 - l)(1755a;6 - 2249a;4 + 721a;2 - 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)=\\J xkK(x)dx\ (j K2(x)dx) j l3k(K) V(K) který lze zkráceně psát T{K) = (\l3k{K)\2v+1V{Kfk-^f,{2k+1). 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: So2- Koptfi,2(x) = f í1 - x2)I[-i,i](x) Si3- Koptth3(x) = lfx(x2 - l)/[_iii](a;) S24- Kopt,2A(x) = -±§(x2 - l)(5a;2 - l)/[_i,i](a;) Přehled vybraných optimálních jader je uveden v tabulce 1.2. Obecný vzorec pro tvar optimálních jader lze nalézt např. v [3]. Jádra Kovt^^ a Koptt2,k se používají pro odhad první a druhé derivace hustoty (viz kapitola 3.1) a z toho důvodu pro ně zavedeme dodatečné označení respektive K^2\ 8 Užitečný při odvozování statistických vlastností odhadů bude i následující pojem. Definice 1.2. Pro jádro K e Svk definujeme kanonický faktor _ ( V(K) \ ^ Vk~\Pl{K)) ' Shrnutí Funkcionály závisející na jádře K (3k{K) = J xkK(x)dx V(K) = J K2(x)dx T{K) = {\h{K)\^V{K-f-^)^ 8vk = (^fy)2fc+1 Jádra s minimálním rozptylem minimalizují funkcionál V{K). Optimální jádra minimalizují funkcionál T(K). Podrobnější popis jader a vzorců s nimi souvisejících lze nalézt v toolboxu Matlabu. Dodatky a cvičení 1. Tvrzení: Nechť K G S^+i^+i je jádro s minimálním rozptylem a KopttVtk £ Svk je optimální jádro. Pak platí K'opt,vAx) = K(x)> x e I-1' x] Důkaz viz [3]. Jako příklad lze uvést jádro Koptfi^{x) = |(1 — x2) G Sq2 a K'opt02(x) = ^1,3(2;) = -|x e Si3. 2. Nechf K e Svk, v e N, pro S > 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 / a g, pro které platí J f2 (x) dx < 00 a J g2 (x) dx < 00. Konvoluci f * g definujeme vztahem (f*g)(x) = J f(ť)g(x-ť)dt. Vlastnosti konvoluce • f *9 = 9*f, • / * (9 * h) = (f * g) * h, • f*(g + h) = f*g + f*h. Ukažte, že pro jádrovou funkci K G £02 platí vztah fl 3 = 0, (K *K)(x)xjdx = < 0 j = l, [2í32(K) 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í, Yi),i = 1,..., 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)+ei, i = l,...,n, (2.1) kde m je neznámá regresní funkce, Xí, 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 Eei = Q, var Eí = a2, i = \,...,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 Xi,i = \,...,n, platí 0 < x\ < x2 < ' ' ' < xn < 1. Cílem regresní analýzy je nalézt vhodnou aproximaci m 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, a f 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 Atixí 1 (1 + cos 0,67rei)2 " kde body xí = i/100, i = 1,..., 100, a chyby Si, i = 1,..., 100, mají normální rozdělení N(0; 0,25). (Data jsou v tabulce 6.1.) 2 1 o -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 sub intervalu, tj. ŕh(x) É YtI[Bj] (x* i=l_ n kde I[Bj] 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: [x — h,x + h], h > 0, přesněji E Yil^c-h^+h] (xi) (2.3) E I[x-h,x+h](Xi) Obr. 2.3 ilustruje aplikaci této metody na simulovaných datech příkladu 1.1. Uvedené metody 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 ffi(x,h) = YlWi(x,h)Yi, (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 G Sok, k je sudé číslo, položme Kh(-) = ^-řr(^). Mezi nejznámější typy jádrových odhadů regresní funkce patří: 1. Nadarayovy-Watsonovy odhady (1964) n J2 Kh(x - Xi)Yi mNW(x,h) = 2=i-, J2 Kh(x - Xi) i=l 2. Priestleyovy-Chaovy odhady (1972) 1 n rhpc{x,h) = - Kh(x - Xí)Yí, n A—' i=l 3. lokálně lineární odhady (Stone 1977, Cleveland 1979) 1 {s2(x, h) - ši(x, h)(x - Xi)}Kh(x - Xí)Yí (*,/0 = -£- n ~[ Š2(x, h')šo(x, h) — ši(x, h)2 kde 1 n šr(x) = - y^(x - Xi)rKh(x - xi), r = 0,1, 2, 71 ^-* n ■ i=l 4. Gasserovy-Můllerovy odhady (1979) i{x,h) = Y,Yt j Kh(x-ť)dt, mGM( kde «o = 0, Si = (xí + Xí+i)/2, sn = 1. Tento odhad je konvolučním typem odhadu. Úmluva. Uvedené jádrové odhady budeme zapisovat ve tvaru 1{x,h)=YJW(lA\x,h)Yl, =i kde váhy W^A) a A značí příslušný typ odhadu NW, PC, LL, GM. V dalším textu budeme zkráceně psát Wi. 13 V mnoha aplikacích je užitečný zejména Nadarayův-Watsonův odhad mjvw- 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 j j—t Obrázek 2.4 ilustruje konstrukci odhadu v bodě xq, který je založen na pěti pozorováních (xi ,Y{), ..., (x5, Y5) (černé křížky). Parabola reprezentuje Epanečnikovo jádro Kh a kroužky znázorňují hodnoty vah Wi = Kh{xo — x i) j J2i=i Kh{xo — x i) pro i = 1,..., 5. Výsledný odhad regresní funkce m v bodě xo je označen hvězdičkou. 1.5i-1-1-1-1-1-1-1 x5 xO+h Obrázek 2.4: Ilustrace Nadarayova-Watsonova odhadu v bodě xo Jádrový odhad není definován pro Eľ=i Kh(x — x i) = 0. Jestliže nastane případ „0/0", pak klademe mjvw(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 | X'~X3 | > 1 pro xi Xj, a tedy hodnota jádra v těchto bodech je rovna nule, s výjimkou bodu xí, v němž dostáváme odhad mNW(xi, h) K(0) Yi. To znamená, že při malé šířce vyhlazovacího okna (h Podobně pro velké hodnoty h je výraz | Xi~X|7 stejnou hodnotu jádrové funkce K ( Xi~^:' ) ~ K(0) a dostaneme tak průměr dat 0) odhad reprodukuje data (viz obr. 2.5(a)). 0, tedy pro všechny body plánu máme mNW(xi, h) Y.K{Q)Y3 K{Q)Y.Y3 „ nK(0) Tedy velká šířka okna (h —> oo) 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í, neboť i asymptoticky optimální odhady obsahují poměrně značné „množství šumu" a to nechává prostor pro subjektivní posouzení. 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 2F *~ 0 0.2 0.4 0.6 0.8 1 Obrázek 2.6: Optimální odhad, h = 0,0816 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)a2(x) _ ,, /V(K)a2(x)~ m{x,h)-Ul.9yl ^ K >,m(x,h)+Ul_*^ ^ . kde ui_a je (1 — a/2)-kvantil standardního normálního rozdělení a odhad rozptylu v bodě x je dán vztahem n a2(x) = J2W*(x>h){Y*-Mx,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 (MSE) odhadu m v bodě x, která je obecně dána vztahem MSErífi, h) = E(m(x, h) - m{x)f. 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 Upravíme tento vztah = Em2(x, h) — 2m(x)Eřh(x, h) + m2(x) = (Em{x,h) - m{x))2 + Eíh2{x,h) - (Em(x,h))2, (2.5) N-v-' V-v-' bias2 var což znamená, že střední kvadratická chyba může být vyjádřena jako součet rozptylu odhadu var m(x, a čtverce vychýlení bias2 fh(x, h) (viz obr. 2.8). Tento rozklad rozptyl-vy•chýlení usnadňuje analýzu vlastností odhadu. 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 ffiNWr íhpcr íhLLr fhcM jsou asymptoticky ekvivalentní (viz např. [8,14]). Z tohoto důvodu budeme dále podrobněji zabývat Priestleyovými-Chaovými odhady, které budeme psát bez uvedení označení PC, tedy: ŕh(x, h) a Wi. Připomeňme, že pro Priestleyovy-Chaovy odhady je váhová funkce tvaru Wi(x, h) = Kh(x - Xi) = ±tf (^íi). Pro další výpočty budeme předpokládat: (i) Jádrová funkce K je sudou funkcí na intervalu [—1,1], tj. K G S"o2, 16 (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 > n0, kde no je pevné, (iv) m G C2 [0,1], (v) Xi = ±, i = l,...,n. Je zřejmé, že pro n —> oo platí1 Em(x,h) = -y^Kh(x-xl)m(xl)= -K(^— )m(t) dt + Orn"1). (2.6) n L—' In h V h ) i=i Jl> Symbolika O a o viz příloha 1. Nechť u = odtud dt = -ftdu a tedy s využitím Taylorova rozvoje /x/h K(u)m(x — hu) du + 0(n_1) -{l-x)/h rx/h r u2h2 i K{u)\m{x) - uhm'{x) + ——m"(x) du + o(h2) + O^"1). (2.7) (i-x)/h L 2 J Podle výše uvedených předpokladů platí h < x < 1 — h, tedy x/h —> oo a — (1 — x)/h —> —oo pro /i —> 0. Odtud, s využitím faktu, že nosičem funkce K je interval [—1,1], plyne ľ1 ľ1 h2 ľ1 Em(x,h)=m(x) i K (u) du-hm'(x) i uK(u)du-\--m" (x) u2K(u) du + o(/i2) + 0(n_1). Celkem dostaneme i j-i 2 j _! h2 Eíh{x, h) = m(x) + —p2{K)m"(x) + o(h2) + C^n"1) h2 m{x) + —l32{K)m"{x). Podobně pro rozptyl platí varm(a;, h) = E(ŕh(x, h) — Eíh(x, /i))2 íl n 1 71 = e( - ^Kh(x - XijYi--^Kh(x - Xi)m(xi) , n z—' n x i=l i=l 1 / n —E\ž^Kh{x - Xí)(Yí - m(xi)) z vlastností (2.2) plyne EKh{xi)Kh{xj)eiej = 0 pro i ^ j, tedy 1 ™ i=l Opět s využitím substituce u = a vztahu 0(n_1) = o((n/i)_1) můžeme pro n —> oo psát n2 r1 n2 r1 n2 varm(x,/i) = — / K2(u) du + oŕŕn/i)"1) = — / íí2(u) du + oŕŕn/i)"1) « — Vŕif). n/i J_1 n/i J.j^ n/i Tímto jsme dokázali následující větu o tvaru střední kvadratické chyby. 1Jedná se o přibližný výpočet integrálu - viz Dodatek na str. 27 17 Věta 3.1. Nechf jsou splněny předpoklady (i) — (iii), pak střední kvadratická chyba nabývá tvaru MSE fh(x,h) = ^V(K) + /i4/322(if)i(m"(a;))2 +o(hi + (n/i)"1). (2.! 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(7i) = í MSE(x,h)dx = AMlSE(h) + o(h'í + (nh)-1). Jo AMISE je tvaru a2 h4 AMISE(7i) = AMISE m{-,h) = —V{K) + —Pl(K)V{m"), (2.9) AI V AI S B kde V (m") = (m"(x))2 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 dAMISEW_^) + dh nh? > \ n položíme ji rovnu nule a vyjádříme h _ a2V{K) ■ppt,0,2 - np2(K-)V(mii /C.0.2 = 7yT777^777~~í7\ ■ (2-10) Poznámka 3.1. Tento výpočet vede k nalezení minima AMISE, protože platí d2 AMISE > 0. dh? Poznámka 3.2. Jestliže jádro K náleží do třídy Sok, pak AMISE je tvaru AMISE(/i) = ^V(K) + J_h2kp2(K)V(m^) a pro optimální vyhlazovací parametr h platí (2.11) kde V{mP^) = (m(fe)(a;))2 dx, podrobněji např. [7]. Nyní uvedeme důležité lemma, které ukazuje zajímavou vlastnost vyhlazovacího parametru. Lemma 3.1. Pro hovt$^ platí AIV(V*,o,fc) = 2kMSQ{hopt^k). Důkaz. Viz cvičení. □ 18 Vztah (2.12) pro optimální šířku vyhlazovacího okna ukazuje, že řád konvergence optimální šířky vyhlazovacího okna hopt,o,2 závisí na řádu jádra k, tedy pro k = 2 je 0(n-1/5). Rovnice (2.12) má pouze teoretický charakter, protože hodnota hoptt0,2 závisí na neznámých veličinách a2 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.10) do vztahu (2.9) pro AMISE, dostaneme AMISE(Vť,o,2) = ^(a2V(K)ý/5(/322(K)V(m"))1/5n-4/5. (2.13) Lze ukázat (viz cvičení), že pro jádra K G Sok je AMISE(/i) = O (n~šftr j. To znamená, že s rostoucím k se zvyšuje asymptotická rychlost konvergence. Ale není zcela jasné, zda tato zvýšená rychlost konvergence vede již k zlepšení pro konečné rozsahy výběrů, neboť ostatní veličiny se rovněž mění s k. 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ě, nebof křivka AMISE(/i) 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, 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.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 (viz [9]) nebo reflexní metodou (viz [3]). 4 Volba jádra Volba jádra není z asymptotického hlediska podstatná, jak je zřejmé z faktu (2.13). Je vhodné zvolit optimální jádro, které minimalizuje fukcioná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 £02 a Soí, lze je vybrat z tabulek pro Sok- 5 Volba vyhlazovacího parametru 5.1 Metoda křížového ověřování Jednou z nejrozšířenějších a nej použí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 19 regresní funkce (2.4), v němž vynecháme i-té pozorování: n rh-i(xi, h) = ^2 WfXxi, h)Ye, i = l,...,n. 1=1 Funkce křížového ověřování je definována takto 1 n CV{h) = -Y^(íh-l{xl,h)-Yi)2 (2.14) n ■ 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 = hCY = arg min CV(/i). h£Hn Při praktických úlohách hledáme minimum na intervalu Hn = [a^n-1/(2fe+1) 5 b^rC1!(2fe+1)]/ jehož tvar plyne ze vztahu (2.12), přičemž ak,bk jsou konstanty (0 < cik < 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 5.1. Někdy se místo chyby MISE používá průměrná střední kvadratická chyba AMSE (average mean square error) 1 n AMSE(/i) = -E^2{fh{xl,h) - m{xi)f 71 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. Pro střední hodnotu funkce CV(/i) platí 1 n ECV(h) = -EY'(m(xl,h)-m(xl))2+a2. AMSE Důkaz. Funkci křížového ověřování lze rozepsat 1 " 2 CV(/i) = - ^(m_i(xi, h) - m(xi) - (Yi - m(xi))) ^ n \ 71 1 71 - y^irh-i(xi, h) - m(xi)) - 2- V" (rh-i(xi, h) - m(xi))ei + - V^ef i=l i=l i=l = - ^(m_j(xi, h) - m(xi))--'ŽZ^wZ wiix^ h)Ye - m(xi)) + ~ ~Y] n z—' n z—' v-'—' / n ■ i=l i=l í=l i=l Střední hodnota E CV(/i) je rovna součtu tří hodnot. Předpokládejme, že m_j(xi, h) « fh(xi, h), pak první ze sčítanců je roven přímo AMSE(/i): 1 n -Ej2{™-i(xi>h) -m(xt))2 = AMSE(/i). 20 Dále víme, že Yt — m(xe) + ei, tedy můžeme psát: n n =1 1=1 =--E n -E^2£i[^2 We(xi,h)(m(xe + et) - m(x =1 i=\ n n n n =i i=i =i i=i =i i=i =i i=i = 0, Stejně jako pro druhý sčítanec, i pro třetí sčítanec využijeme vlastnosti (2.2): 1 n □ Tento výsledek znamená, že minimalizace ECV(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. 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]). 21 Obrázek 2.11: Simulovaná data (x) s jádrovým odhadem regresní funkce (/icv = 0,1158) (—) a původní funkcí (—) 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 poznámky 3.2 je známo, že AMISE a hoptfitk jsou tvaru AMISEW = £viK) + ^rk(K)V(mi% h™k = «W Ze vztahu pro hoptto,k vypočteme V(m(fe)) a dosadíme do vztahu pro AMISE, použijeme vztahy Sok = (WK)) ^ T{K) = (/W)^))57*1 T{K) = 5ltPl{K) dostaneme vyjádření AMISE AMISE(Vť,o,fc) = aT!:h+1)ÔOkT(K), ve kterém jsou parametry K,hak separovány, což umožňuje vybrat tyto parametry simultánně, Právě tento vztah je základem automatické procedury. Označme a množinu vhodných řádů k označme I(ko) = {2j,j = 0,...,[!%]}, kde [z] značí celou část čísla z. Procedura pak probíhá v pěti krocích: 22 Krok 1 Pro k G I(k0) najděte optimální jádro -řropt,o,fc £ Sok, které je dáno tabulkou 1.2 a vypočtěte kanonický faktor Sok- Krok 2 Pro /c G I(ko) a Koptt0,k G iSoifc najděte optimální vyhlazovací parametr hoptfi,k- Krok 3 Pro /c G I(k0) 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. n ~hopt,0,k) = W^X' ~hopt,0,k)Y>- i=l Příklad 6.1. Aplikace procedury na data z příkladu 1.1. Maximální řád jádra zvolme k0 = 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 hoptfitk- k Kopt,0,k Sok h L(k) 2 -1) 1,7188 0,1158 0,0648 4 32 ~ l)(7x2-3) 2,0165 0,2446 0,0575 6 105 12 256 v - l)(33a;4 - 30a;2 + 5) 2,0834 0,3575 0,0574 8 315 12 4096 v - l)(715a;6 - - lOOlx4 + 385a;2 - - 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 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 (Kopt^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. 13 12 11 10 9 8 7 6 5 1880 1900 1920 1940 1960 1980 Obrázek 2.13: Úroveň hladiny Hurónského jezera 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 (Äopť,o,2), hr, 0,8826 (^,0,10) 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 60 80 100 120 140 160 18 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 18 (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 n íh(x, h) = ^ Wi(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 biasmíx, h) « —(32(K')m"(x), vaim(x,h) ~ —V(K). 2 nh Asymptotická střední kvadratická chyba jádrového odhadu regresní funkce s jádrem řádu 2 je AMISE(/i) = ^V(K) + ^h4f3Í(K)V(m"). Optimální vyhlazovací parametr vzhledem k AMISE pro k = 2 je tvaru íl A 'AMISE - nß2(K)V(m„) s řádem konvergence n-1/5. Metoda křížového ověřování pro odhad optimálního vyhlazovacího parametru 1 n CV{h) = -"y{fh_l(xl,h)-Yi)2 => hcv = argminCVÍ/i). 71 -^ 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 ffiNW,fnLL,fhcM dokažte, že „množství" vyhlazení pomocí jádra K s vyhlazovacím parametrem h je stejné, jako „množství" vyhlazení jádrem Kg s parametrem h* = h/s,q. ŕh(x, h, K) = m(x, h*,K$). 2. Dokažte vztah pro střední hodnotu odhadu (2.6). 3. Dokažte vztah (2.13) pro obecné k, tj. že platí AMISE(/lopt,0lfc) = n"2fc/(2fc+1) (a2V(K))2k/(2k+1) ^2(K)(k\)-2V(m^) (2fc)dftr Návod: Předpokládejme, že funkce m je spojitá a má spojité derivace až do řádu k včetně, tj. m G C* [0,1]. Užitím Taylorova (až do řádu k) rozvoje upravíme vztah pro střední hodnotu odhadu podobně jako ve vztahu (2.7). Užijte poznámky 3.2. 4. Dokažte vztah z lemmatu 3.1. 26 / V(K) \ l/(2fc+l) 5. Dokažte, že příspěvek jádra A'(sofe, Sok = ( ^2 (Jf 1 , ^ £ Sofc, k oběma částem chyby AMISE je stejný, tj. platí rovnost íofe / r^nk Sok \J Sok 6. Aplikujte metodu křížového ověřování a automatickou proceduru na simulovaná i reálná data. Dodatek Výpočet integrálu, U = ^, i = 1...., n, ni G(ť) dt n-l rti+1 = J2 G(t)dt s využitím Taylorova rozvoje funkce G(t) n-l = J2 / (G{ti+i) + (* - + o(l)) dť n-l n-l rti+1 = J2 G(ti+i)(ti+í -U)+ ^2 (t- tl+1)G'(U+1) dt + o(l) i=0 " ^ ' i=0 ťi 1 11 n_1 (i - i \2 1 n 1 n—1 = -EG(ť*) + 2^EG'(ť*+i) + o(1)- i=l i=0 Za předpokladu |Gf/(íi+1) | < M, pak platí r1 i n / G(í)dí = -5]G(íi) + 0(n-1). ■/o n „--i 27