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 " Fn{x) = - ^/(.^(Xi). 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. 49 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 0-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, h)dt=—Y,l 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). 50 -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<-\, = { \(-x3 + 3x + 2) \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 . 51 Spočítejme nejdříve hodnotu EF(x, h) v bodě ieR: EF(x,h) = Jw(^^jf(y)dy /l f- oo W(t)f(x-ht)dt + h / W(t)f(x-ht)dt. -oo Jl Předpokládejme dále, že F e C2. Označme první integrál li a druhý I2- Integrál ii počítáme metodou per partes a využijeme vlastnosti funkce W(t) h = h í W(t)f(x - ht) dt J — OO u = VK(í) w' = W(í) = if(í) í/ = f(x - ht)h v = -F(x - /ií) = [-F(x - ht)}1^ + J F(x- hť)W'(t) dt = -F(x - h) + J K (t) F (x - ht) dt. (4.2) Dále použijeme Tayloruv rozvoj h2t2 F (x - ht) = F (x) - htF'(x) + —F" {x) + o(h2), tedy Ji = -F(x - /t) + F(x) + ^F"(x)h2fo{K) + o{h2). Počítejme nyní integrál I2: /oo W(t)f(x - ht) dt. Uvažujeme-li substituce x — ht — z, dostaneme r- —OO r-X — h h = -l f(z)dz = f(z)dz = F(x-h). J x —h J—in f ty Vychýlení odhadu je tedy tvaru biasF(x,h) = ^F"(x)h2p2(K) + o(h2). Poznámka 3.1. Vztahy (4.2) a (??) dávají zajímavý vztah pro vychýlení EF(x,h) - F(x) = / if(í)F(x-/ií)dí-F(x). J-i Odtud plyne EF(x, h) = y K (t) F (x - ht) dt a také (z Taylorova vzorce) EF(x,h) =F(x) + o(h). 52 Nyní dokážeme tvar rozptylu. varFfo/i) = - ( EW2 (^-r^] - E2W x-X h Zde E2W (^jr) = {EW (2^rL)Y = (F(x) + o{h)f. Počítáme tedy pouze integrál I3: — |substituce: x — y — th\ ( \ 1 f- OO W2(t)f(x-ht)dt + h / f{x-hť)át oo Jl v =F{x-h) / První integrál počítáme metodou per partes a máme 1 2 ľ 1 J3 = - [_F(x - ^í)W2(í)l * + - / F(x - /ií)VK(í) W(í) dí + -F(x - h) tri <- -I 1 irt I tri --F(x -h) + - f F(x- ht)W(t)K(t) dí + -F(x - h) n n n použijeme nyní Taylorův rozvoj 2 W(t)K(t) (F(x) - htF'(x) + o(h)) dí -F(x) / W(t)K(t) dí - -hF'(x) / tW(t)K(t) dí + o - '/i o I - n užitím vlastností funkce W a F'(aľ) — f (x) dostaneme F (x) - h f (x) [í- W2{ť)dt Rozptyl je tedy tvaru var F{x, h) — — F (x) - h f (x) ( 1 - / W2(t)dt i + o{-)-(F(x)+o(h)Y = -F(x)(í - F(x)) --f (x) (í - f W2(t) dí ) + o (-n n V J-i J \n Výše uvedené výsledky můžeme nyní zformulovat v následující větě: Věta 3.1. Nechí F e C2, h ->• 0, nh ->• oo pro n ->• oo. PflA: MSEF>,/0 = -F(a:)(l-F(a;))-Í/i/(a;) fl- / T^2(í) dt)+^-(F"(x))2h^2 (K)+o (- + \ŕ n n \ J_i J 4 \n (4.3) Globální pohled na kvalitu odhadu lze získat prostřednictvím střední integrální kvadratické chyby (MISE). 53 Věta 3.2. Nechí F e C2, V(F") = / {F" [x])2 dx < oo, K e S02, lim„^oo h—Oa lim^oo nh = oo. Pak MISEF(-,/i) = - />f1(x)(l-f1(a;))da;-ci- + C2/i4 + o( - + /i4) , (4.4) n J n \n J kde ci = i - y' w2(í) dí, c2 = i/3|(i<:)y(F"). Naším cílem je nalézt takovou hodnotu vyhlazovacího parametru, pro kterou bude MISE nabývat minimální hodnoty. Ale uvedený tvar MISE není pro takovou analýzu vhodný, a proto (stejně jako při odhadu hustoty a regresní funkce) budeme uvažovat asymptotickou střední integrální kvadratickou chybu AMISE, která v tomto případě je tvaru: 1 ľ h AMISE F (■, h) = AMISE (h) = - / F (x) (í - F (x)) dx - Cl- + c2h4 . (4.5) n J n "-v-' 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 /iopt.o.2 = n"1/3 (j^ = Oln-W) (4 6) a pak _L^4/3„-4/3 AMISE(/lopt,o,2) = \J F{x){\ - F{x))dx-4j^ (j)' (4.7) Poznámka 3.2. Optimálni hodnota vyhlazovacího parametru pro odhad distribuční funkce je řádu n-1/3, zatímco pro odhad hustoty s jádrem K G So 2 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[_\^(x), • kvartické K(x) — ±|(1 - x2)2I[_ltl](x), • triweight K(x) = §§(1 - x2)3/[_ia] (x). 5 Volba vyhlazovacího parametru 5.1 Metody křížového ověřování Metody křížového ověřování patří k nejužívanějším metodám pro volbu vyhlazovacího parametru. Zde uvedeme pouze metodu navrženou A. Bowmanem (1998). Funkce křížového ověřování je v tomto případě tvaru CV^ = \ Ž / {h-ooA - F-i(x, /i)) 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. 54 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) = yg (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, p2{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. 55 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^,4 G S24. Pak n n n / V V Ä = n"1 Y, /"(*, Ä) = n-2/*-3 E E ( 2 — 1 2 — 1 J—1 kde podle vztahu (3.9) je Pak _lnl/9524, lopt,2,4 — 1U T—"opt, 0,4-004 1 £ = —jim 4' Shrnutím předchozích úvah dostaneme proceduru pro odhad distribuční funkce F: Krok 1 Najděte optimální vyhlazovací parametr /iOpt,0,4 pro odhad hustoty s optimálním jádrem Kopt,0,4 G »$04. Krok 2 Najděte optimální vyhlazovací parametr hopt2,4 pro odhad druhé derivace hustoty podle vztahu (3.9) s & = 4 a optimálním jádrem e S24. Krok 3 Vypočtěte odhad funkcionálu i/>i s využitím hodnoty /i0pt,2,4 získané v kroku 2. Krok 4 Vyčíslete optimální hodnotu vyhlazovacího parametru / \ 1/3 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 56 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*) = \ £ 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)(í - 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/jgij^ CV(/i) i—l • metoda maximálního vyhlazení -"/3 (lák)1,3 ^ • plug-in metoda 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.6) a (4.7). 4. Odvoďte tvar vyhlazovacího parametru podle metody maximálního vyhlazení pro Epa-nečnikovo jádro. 58 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. 59