Skalární skór Zdeněk Fabián Ústav informatiky AV ČR, Praha
Abstrakt. Po stručném přehledu základních inferenčních funkcí matematické statistiky zavedeme jednu novou a porovnáme její přednosti a nedostatky s ostatními.
1
Úvod
Výsledky pozorování náhodných veličin lze reprezentovat pomocí ukazatelů a čísel, které shrnují nějaké obecnější tendence. K tomu slouží statistické charakteristiky dat typu ’středu’ a ’poloměru rozptýlenosti’, šikmosti či špičatosti, a statistické analýzy, např. korelace či regrese, kterými se hledají souvislosti mezi náhodnými veličinami. Charakteristiky populace jsou získávány z pozorovaných dat obvykle jako řešení rovnic typu Ψ(data; θ) = 0 kde θ je vektor parametrů zvoleného statistického modelu a Ψ statistická inferenční funkce. V tomto textu popisuji přednosti a nedostatky nejužívanějších inferenčních funkcí a představuji jednu novou, uvedenou do statistické literatury pod jménem core funkce v r. 2001 (viz [3]), kterou dnes nazývám skalárním skórem. Skalární skór je užitečný sice univerzálně, ale zdá se, že nejvíce pro rozdělení, jejichž hustoty neklesají dostatečně rychle nule, t.zv. rozdělení s těžkými chvosty, která nemusí mít střední hodnotu ani rozptyl.
2
Základní pojmy a značení
R značí reálnou osu a X ⊆ R otevřený (konečný či nekonečný) interval. P značí pravděpodobnost. Spojitá náhodná veličina X má rozdělení F , popsané distribuční funkcí F (x) = P (X < x) na nosiči X , jestliže je její pravděpodobnostní hustota (dále jen hustota) f (x) = dF (x)/dx na X kladná a na R − X rovna nule. Místo F s nosičem X někdy říkáme F na X . 1
Budeme se zabývat úlohami, kdy je možné předpokládat, že pozorovaná data Xn = (x1 , ..., xn ) jsou náhodným výběrem z F , t.j. jsou realizací spojitých náhodných veličin X1 , ..., Xn , nezávislých a stejně rozdělených podle (neznámého) rozdělení F zkoumané populace (a nabývajícich hodnot v X ). Dat je obvykle méně, než by bylo potřeba k určení hustoty populace např. aproximací histogramu a je nutno použít jiných metod. Parametrický přístup spočívá v předpokladu, že F je člen Fθ0 předpokládaného modelu, rodiny {Fθ : θ ∈ Θ}, kde Θ ⊆ Rm , s hustotami f (x; θ), které jsou spojité a mají v x i θ potřebné derivace. Často užívanou modelovou rodinou je normální rozdělení (Tabulka 1) s nosičem R a s parametrem θ = (µ, s). Jako µ ∈ R se obecně označuje parametr polohy, t.j. x-ová souřadnice maxima hustoty (mód), a jako s ∈ (0, ∞) značíme parametr měřítka, popisující ’rozptýlenost’ hodnot kolem ’středu’ µ. Hodnota estimátoru Tn (X) = T (X1 , ..., Xn ) parametru θ0 , určená z náhodného výběru Xn , odhad θˆn , kolísá od výběru k výběru: Tn (X) je náhodná veličina, která má určité rozdělení, které vypovídá o přesnosti odhadu. Toto rozdělení se obvykle nedá odvodit pro pevné n a je fajn, když určíme alespoň asymptotické, limitní rozdělení pro n → ∞. Známe-li je, lze pak doufat (a případně podepřít výsledky počítačových simulací metodou zvanou bootstrap), že přibližně platí i pro taková n, která máme k dispozici. Pak můžeme konstruovat intervaly spolehlivosti, určovat míru rizika že se na základě daného odhadu rozhodneme chybně atd. V dalším nebudeme rozlišovat mezi estimátorem (náhodnou veličinou) a její konkrétní hodnotou pro určitý náhodný výběr (odhadem). Estimátor je konzistentní když θˆn → θ0 pro velká n (pokládá se za nutnou vlastnost). Konzistentní odhad θˆn je asymptoticky normální, AN (θ0 , v 2 /n), když√pro velká n má θˆn − θ0 přibližně normální rozdělení s parametrem θ = (0, v/ n). Odhady (kromě výběrového průměru) značíme v dalším stříškou a index n vynecháváme.
3
Identická inferenční funkce
Předností identické inferenční funkce Ψ(x) = x je samozřejmě její jednoduchost. Charakteristikami rozdělení jsou momenty Z k µk = EΨ (X) = xk f (x) dx (1) X
2
a centrální momenty νk = E(X − m)k . Pokud integrály (1) konvergují (momenty existují), populace má střední hodnotu m = µ1 a rozptyl σ 2 = ν2 . Veličina ν3 /σ 3 charakterizuje šikmost a ν4 /σ 4 špičatost rozdělení. Předpokládejme, že střední hodnota a rozptyl populace s rozdělením F existují. Data pak lze charakterizovat konečnou aproximací µ1 , výběrovým průměrem n 1X x¯ = xi (2) n i=1 ¯ = 1 Pn Xi ), a výběrovou (připomeňme, že se jedná o náhodnou veličinu X i=1 n středněkvadratickou odchylkou σ ˆ , kde n
1 X σ ˆ = (xi − x¯)2 n − 1 i=1 2
(3)
je výběrový rozptyl. x¯ a σ ˆ 2 jsou konzistentními odhady veličin m a σ 2 , a pokud je σ známé, platí podle centrální limitní věty že x¯ je AN (µ, σ 2 /n). Tento hluboký, dobře známý výsledek teorie pravděpodobnosti má nepříjemný důsledek: x¯ se zhusta považuje za typickou pozorovanou hodnotu, což asi není ono v případě šikmých rozdělení s nosičem (0, ∞). Jednoduchou mírou závislosti (korelace) náhodných veličin X, Y je první smíšený moment, jehož výběrovou podobou je Pearsonův korelační koeficient n
rˆXY =
1 X (xi − x¯) (yi − y¯) n i=1 σ ˆX m ˆY
(4)
kde σ ˆX , σ ˆY jsou příslušné výběrové středněkvadratické odchylky. Nejjednodušší inferenční funkce, Ψ(x) = x, má samozřejmě své mouchy. Generuje jednoduché charakteristiky rozdělení a dat, které nepříliš dobře charakterizují data z šikmých rozdělení a které jsou zcela nepoužitelné pro data z rozdělení s těžkými chvosty, pro které integrály (1) nekonvergují a výběrové průměry nemají asymptoticky normální rozdělení. Charakteristiky (2), (3) i (4) zásadním způsobem ovlivňují odlehlá data: výběrové charakteristiky nejsou robustní.
4
Skór pro parametr polohy
Podle principu maximální věrohodnosti je nejvhodnějším kandidátem pro odQn ˆ had parametru θ takové θ, které maximalizuje věrohodnost i=1 f (xi ; θ), t.j. 3
pravděpodobnost současného výskytu x1 , ..., xn v předpokládaném modelu. Protože f (x) a log f (x) jdou ’paralelně’, máme θˆ = arg max θ
n X
log f (xi ; θ).
(5)
i=1
Pro rodinu s nosičem R a hustotou f (x−µ) dostáváme maximálně věrohodný (ML) odhad parametru polohy z rovnice n X
Uµ (xi − µ ˆM L ) = 0,
(6)
i=1
kde Uµ (x − µ) =
∂ log f (x − µ) ∂µ
(7)
je t.zv. skór pro parametr polohy. Podle (6) je odhad µ ˆM L (módu) konzis2 tentní a za nepříliš omezujících podmínek je AN (µ, σµ /n), kde σµ2 = 1/Iµ a kde Iµ = EUµ2 je Fisherova míra informace o parametru µ. Veličina σµ2 /n je pro dané n Cramér-Raova dolní mez rozptylu odhadu: µ ˆM L je vydatný (eficientní). Jasná zpráva, jako skalární inferenční funkci můžeme uvažovat funkci S(x − µ) = Uµ (x − µ). V obecném případě rozdělení s nosičem R a hustotou f (x; µ, θ2 , ..., θm ) je S(x; θ) = Uµ (x; µ, θ2 , ..., θm ). Několik rozdělení s parametrem polohy uvádíme v Tabulce 1. Z důvodu, který bude jasný později, značíme hustotu g(y; µ, s) a položíme y−µ u= . s Z tabulky je patrné, že identická inferenční funkce je vlastně skór pro parametr polohy normálního rozdělení. Skór je v tomto případě neomezený pro y → ±∞. Skór Uµ Gumbelova rozdělení je neomezený pro y → ∞, logistické rozdělení má skór omezený a Cauchyho dokonce redescentní, t.j. klesající pro velká |y| k nule. Tabulka 1. Několik rozdělení s parametrem polohy. 4
Rozdělení
g(y; µ, s)
Uµ (y; µ, s)
− 12 u2
normální Gumbelovo logistické Cauchyho
√1 e 2πs 1 u −eu e e s 1 eu s (1+eu )2 1 1 πs 1+u2
1 u s u
1 (e − 1) s 1 eu −1 s eu +1 1 2u s 1+u2
Poslední dvě rozdělení tabulky mají těžké chvosty, ale jejich skóry pro parametr polohy jsou omezené. V povaze modelu tedy je generovat odlehlá data, která si ale sám ’reguluje’ ve smyslu že µ ˆM L k nim není příliš citlivý, odhad je robustní: přidáme-li k výběru Xn obrovské xn+1 , µ ˆn+1 ≈ µ ˆn . Poznamenejme zde že Cauchyho rozdělení sice nemá střední hodnotu, ale má mód y ∗ = µ, který lze považovat za alternativní charakteristiku ’středu’ (centrální tendence) rozdělení, a právě toto my učiníme pro všechna rozdělení s nosičem R (hustota rozdělení na R musí mít maximum, a když je má vícenásobné, lze ji chápat jako hustotu směsi). Rozdělení, jejichž skór pro parametr polohy je u některého z ’konců’ nosiče neomezený, jsou v té části naopak k odlehlým hodnotám v datech citlivá. Zde předpokládáme, že odlehlá data pocházejí z kontaminovaného modelu Fc = (1 − ²)F + ²H,
² ¿ 1.
kde H je rozdělení s velkým rozptylem. Pro takové situace doporučuje robustní statistika modelu Fθ do jisté míry nedbat a použít jako inferenční funkci nějakou spojitou a omezenou ’psí-funkci’ ψ(x) s vlastností Eψ = 0, a namísto (6) řešit rovnici n X
ψ(xi − µ ˆM ) = 0.
(8)
i=1
Odhad µ ˆM je zobecnění ML odhadu, t.zv. M-odhad. Je konzistentní a asympto2 ticky normální, AN (µ, σM /n), kde 2 = σM
Eψ 2 (x) . [Eψ 0 (x)]2
(9)
2 2 Poznamenejme, že poměr σM L /σM je asymptotická vydatnost M-odhadu parametru polohy. ² je vždy menší než 1 (Cramer-Raova mez). Volbou ψ se volí
5
určitý kompromis mezi vydatností a robustností odhadu. V případě kontaminovaného normálního rozdělení se jako psí funkce často používá Huberova funkce pro x > b b x pro |x| < b ψ(x) = −b pro x < −b, kde b je určeno z požadovaného kompromisu mezi robustností a vydatností.
5
Skórová funkce
Parametrický model má obvykle parametr θ s více (2-3) složkami. Derivováním (5) podle θj dostáváme soustavu m rovnic pro m složek θ, n X
Uθj (xi ; θˆM L ) = 0,
j = 1, ..., m
(10)
i=1
kde Uθj (x; θ) =
∂ log f (x; θ) ∂θj
(11)
jsou parciální skóry pro jednotlivé složky. Maximálně věrohodné odhady θˆM L jsou konzistentní, asymptoticky normální a vydatné. Platí, že (θˆM Lj − θ0 ) je AN (0, Jjj /n),
j = 1, ..., m,
kde J = {I −1 (θ0 )}|i,j=1,...,m a I(θ) = {EUθi Uθj (θ)}|i,j=1,...,m je Fisherova informační matice. Vektorová skórová funkce Ψ(x; θ) = [Uθ1 (x; θ), ..., Uθm (x; θ)]
(12)
je inferenční funkcí klasické statistiky a θˆM L tím nejlepším odhadem z hlediska minimálníko rozptylu. Má to však několik háčků. a) Namísto výsledku ve tvaru f (x; θˆM L ) je často pro další analýzy dobré mít několik čísel charakterizujících data. Parametry různých modelů ale mají
6
různou strukturu a různý význam. Ideální by bylo charakterizovat data nějakou typickou hodnotou a ’poloměrem rozptýlenosti’. Nabízejí se ovšem momenty: střední hodnota a středněkvadratická odchylka jako funkce odhadnutého θˆM L . Tento postup se však nepoužívá; momenty jsou často vyjádřeny pomocí speciálních funkcí a momenty rozdělení s těžkými chvosty neexistují. b) Skóry pro různá θj mají různou citlivost k odlehlým hodnotám. Například rovnice (10) pro model s parametry polohy a měřítka s hustotami s−1 f ((x − µ)/s) mají tvar µ ¶ n X xi − µ Uµ = 0 (13) s i=1 ¶ µ n X 1 xi − µ xi − µ Uµ = 1. (14) n i=1 s s Inferenční funkce rovnice (14) je tedy Ψ(ξ) = ξUµ (ξ), což je neomezená funkce, takže simultánní odhady (ˆ µM L , sˆM L ) nejsou robustní ani když je Uµ omezená. Robustní statistika navrhuje v tomto případě použít M-odhady (µM , sM ) určené z rovnic µ ¶ n X xi − µ ψ = 0 s i=1 µ ¶ n 1X xi − µ χ = 1 n i=1 s kde ψ(x) je omezená funkce a χ(x) = dψ(x)/dx. Rovnice je nutno řešit iteračním postupem s vhodnými počátečními robustními odhady parametrů, obvykle se používá µ0 = med(x), což je medián (teoreticky bod x˜, kde P (X < x˜) = P (X > x˜) a prakticky prostřední pozorovaný bod), a s0 = MAD(x) = med{x − med(x)} (pro normálně rozdělená data je s = 1.4785 ∗ MAD(x)). Robustní statistické metody (např. [11,13]) eliminují vliv odlehlých pozorování. Volba psí funkce však částečně eliminuje informaci, která vedla k uvažovanému modelu, který může mít i jiné parametry než jsou parametry polohy a měřítka. Pro modely s jinými parametry nebo více parametry jsou však robustní postupy známy spšee jen výjimečně. c) Robustní verze korelačního koeficientu (4) je µ ¶ µ ¶ n 1X xi − m ˆX yi − m ˆY rˆXY = ψX ψY , (15) n i=1 σ ˆX σ ˆY 7
kde ΨX , ΨY jsou vhodné psí funkce a odhady nšjaké robustní odhady středních hodnot a rozptylu (medián, MAD). Vektorová skórová funkce klasické statistiky je však je pro tento a jiné podobné účely příliš komplikovaná. V praxi se tedy používají korelační koeficienty (4) a v případě možného výskytu odlehlých hodnot (15) nebo koeficient konstruovaný z uspořádaného výběru (Spearmanův). Žádný z nich nevyužívá předpokládané modely náhodných veličin X a Y .
6
Skalární skórová funkce
V této kapitole popíšeme postup, kterým lze zavést skalární inferenční funkci S(x), resp. S(x; θ), která koresponduje s hlavními rysy rozdělení) pro libovolné spojité, parametrické či neparametické rozdělení. Uvidíme, že ji lze zvolit tak, aby pro některá rozdělení byla identická s parciálním skórem pro určitý (nejdůležitější) parametr rozdělení. Nazveme ji skalární skórovou funkcí nebo skalárním skórem. Jejími numerickými charakteristikami budou skórové momenty: zobecněné momenty tvaru Z k Mk (θ) = ES (x; θ) = S k (x; θ)f (x; θ) dx. (16) X
Namísto řešení soustavy (6) získáme odhad θˆSM parametru θ = (θ1 , ..., θm ) jako řešení rovnic n
1X k S (xi ; θ) = Mk (θ) n i=1
k = 1, ..., m.
(17)
Z důvodů popsaných v kap. 4 je přirozené zvolit za skalární skór náhodné veličiny Y s rozdělením G na R s hustotou g funkci SG (y) = −
g 0 (y) g(y)
(18)
vyjadřující relativní změnu hustoty. Řešením y ∗ rovnice SG (y) = 0, a v parametrickém případě SG (y; θ) = 0), je mód rozdělení. Funkce (18) bohužel není vhodným popisem rozdělení s nosičem X 6= R (viz třeba exponenciální rozdělení na X = (0, ∞) s hustotou g(y) = e−y , pro které SG (y) ≡ 1 nebo rovnoměrné rozdělení, kde SG (x) ≡ 0). Hustota f (x) 8
nemusí mít na X maximum (může mít třeba supremum v x = 0). Následující konstrukce byla publikována v [3]: X s rozdělením F na X = 6 R považujeme za transformovanou náhodnou −1 veličinu X = η (Y ), kde η : X → R je vhodné spojité rostoucí zobrazení a Y = η(X) je ’prototyp’ s rozdělením G na R. Rozdělení náhodné veličiny X je F (x) = G(η(x)) a má hustotu f (x) = g(η(x))η 0 (x),
(19)
kde η 0 (x) = dη(x)/dx je Jacobián zobrazení. Za signifikantní funkci rozdělení F považujeme transformovaný skalární skór prototypu, funkci TF (x) = SG (η(x)).
(20)
Z (18) a (19) se snadno odvodí, že (20) lze vyjádřit pomocí hustoty f rozdělení F jako µ ¶ 1 d 1 TF (x) = − f (x) . (21) f (x) dx η 0 (x) Pro porovnání vlastností různých rozdělení na X je třeba zvolit pro daný nosič jedno určité zobrazení, nejlépe takové, aby (21) byla vyjádřena jednoduchými vzorci pro podstatnou část prakticky užívaných rozdělení. Protože velké množství modelových hustot je vyjádřeno pomocí exponenciálních funkcí, je pro tento účel vhodné Johnsonovo zobrazení [12] ½ log(x − a) pro X = (a, ∞) η(x) = (22) x log 1 − pro X = (0, 1) x (jeho zobecnění pro obecný interval je zřejmé). Po dosazení (22) do (21) dostaneme explicitní tvar funkce TF pro různé nosiče X ve tvaru 0 (x) X =R − ff (x) f 0 (x) −1 − (x − a) f (x) X = (a, ∞) TF (x) = (23) f 0 (x) −1 + 2x − x(1 − x) f (x) X = (0, 1). TF (x) a její parametrický tvar TF (x; θ) budeme nazývat t-skór. Těžištěm x∗ rozdělení F nazveme řešení rovnice TF (x) = 0. 9
(24)
Protože podle (20) 0 = TF (x∗ ) = SG (η(x∗ )) = −
g 0 (y ∗ ) , g(y ∗ )
je těžiště rozdělení F obrazem módu svého prototypu G. Takto jsme tedy vyřešili dodnes diskutovaný problém, totiž co je nejvhodnější ’centrální charakteristikou’ (typickou hodnotou) rozdělení na nosiči (0, ∞), zda střední hodnota (nemusí existovat) či mód (nemusí existovat) nebo medián (matematicky špatně uchopitelná veličina). Odpověď zní těžiště: transformovaný mód prototypu (ten, jak jsme viděli v kap. 4, existuje). Uvažujme nyní parametrická rozdělení Fθ na X 6= R. Může se stát, že x∗ je přímo jedním z parametrů. Je-li µ parametrem polohy prototypu Gθ˜ s vektorovým parametrem θ˜ = (µ, θ2 , ..., θm ), parametrem rozdělení Fθ (x) = Gθ˜(η(x)) je θ = (τ, θ2 , ..., θm ), kde τ = η −1 (µ)
(25)
je obrazem parametru polohy prototypu. Označme třídu těchto rozdělení třeba DX . V Tabulce 2 jsou uvedeny hustoty a t-skóry transformovaných rozdělení z třídy D(0,∞) s prototypy z Tabulky 1. τ (25) je transformovaný parametr polohy, který se obvykle považuje za parametr měřítka; my jej však interpretujeme jako parametr polohy rozdělení s nosičem (0, ∞). 1/x je Jakobián transformace, c = 1/s a výraz u z Tabulky 1 se transformuje na log y − log µ ³ x ´c v= = . s τ Tabulka 2. Transformovaná rozdělení z třídy D(0,∞) . Rozdělení lognormální Weibullovo log-logistické log-Cauchy
f (x; τ, c) − 12
2
log v √c e x 2π c ve−v x v c x (v+1)2 1 c πx 1+log2 v
TF (x; τ, c) c log v c(v − 1) c v−1 v+1 2c log v 1+log2 v
Ve třídě rozdělení DX platí vztah, dokázaný v [3], že totiž Uτ (x; θ) = η 0 (τ )TF (x; θ). 10
(26)
Parciální skór pro τ , definovaný pomocí derivace podle parametru, se nám podařilo rozložit na součin Jakobiánu transformace a členu určeného derivací podle proměnné. Ovšem ne každý prototyp musí mít parametr polohy a většina rozdělení na (0, ∞) není z D(00 ∞) (jsou to většinou rozdělení s prototypy, které se neuvažují a neužívají, takže není na první pohled patrné, že je lze chápat jako transformovaná). Rozhodujícím krokem bylo si uvědomit, že τ funguje v (26) nejen jako parametr, ale také jako těžiště rozdělení. Rozdělení, která nejsou v DX , sice nemají parametr, který by byl obrazem módu prototypu, mají však těžiště. V práci [4] je definována funkce S(x; θ) = η 0 (x∗ )TF (x; θ)
(27)
kterou nazýváme skalárním skórem a považujeme ji za obecnou skalární inferenční funkci. Je utvořena podle vzoru (26), proto lze o ní předpokládat, že bude mít podobný význam jako má (26) pro rozdělení z DX . Funkce sice není tak důkladným popisem rozdělení jako vektorová skórová funkce, ale zachycuje jeho hlavní rysy a je to snadno manipulovatelná funkce, a to je ono. Pro rozdělení, která nejsou ve třídě DX , byla (27) (stejně jako TF ) dosud neznámou funkcí. Veličinu ES 2 lze pak chápat jako Fisherovu informaci o těžišti (pro rozdělení z třídy DX je to podle (26) skutečná Fisherova míra informace Iτ (θ) o τ ). Jako charakteristiku variability rozdělení F pak definujeme její převrácenou hodnotu, 1 ω2 = , (28) ES 2 které budeme říkat s-variance (ze score variance). Podobně jako x∗ , ω 2 existuje i pro rozdělení s těžkými chvosty (poznamenejme, že podmínka 0 < Fisherova míra informace < ∞ je obvyklou podmínkou regularity rozdělení). Skalární skór normálního rozdělení s parametry µ a s je S(x; µ, s) = x−µ a protože ES 2 = s2 , platí ω 2 = s2 . Pro normální rozdělení tedy nic s2 nového. V Tabulce 3, kde B značí beta funkci a Γ gamma funkci, uvádíme kromě t-skóru střední hodnotu, těžiště a s-varianci několika dalších rozdělení s nosičem X 6= R, která nejsou z DX . Skalární skór rozdělení s nosičem (0, ∞) je S(x; θ) = TF (x; θ)/x∗ , Paretovo rozdělení má nosič (1, ∞) a skalární skór S(x; θ) = TF (x; θ)/(x∗ − 1) = c(1 − x∗ /x) s použitím (27) a (22). Pro gamma rozdělení s lineárním t-skórem se střední hodnota a těžiště shodují. Zbývající tři rozdělení jsou rozdělení s těžkým chvostem v +∞; 11
jejich skalární skór je omezený (ale všiměte si, že skór Fréchetova rozdělení je neomezený v nule). Pro hodnotu parametru menší nebo rovnou jedné střední hodnoty neexistují, ale jsou nesmyslné velké i v oboru, kde existují, když se Tabulka 3. Těžiště a s-variance některých rozdělení. rozdělení gamma Paretovo beta 2.druhu Fréchetovo
f (x)
γ α α−1 −γx x e Γ(α) c+1
c/x
xp−1
1 B(p,q) (x+1)p+q c τ c −( τx )c ( )e x x
TF (x) γx − α c − c+1 x c[1
qx−p x+1 ¡ ¢c − xτ ]
m α/γ
x∗ α/γ
ω2 α/γ 2
c c−1 p q−1
c+1 c p q
c+2 c3 p(p+q+1) q3 2 2
τ Γ(1 − 1c )
τ
τ /c
hodnota parametru blíží k jedné. Na obr.1 jsou znázorněny hustoty a normované skalární skóry S(x)/ES 2 Weibullova (Tabulka 2) a Paretova (Tabulka 3) rozdělení pro několik hodnot parametru c. Hustoty Weibull
Hustoty Pareto
c=0.75
c=0.5
c=2.5
c=2.5
5
10
15
3
Skalar score Weibull
5
7
Scalar score Pareto
c=2.5 c=2.5
c=0.5 c=0.75
5
10
15
3
5
7
Obr. 1. Skalární skóry dvou rozdělení, která nejsou v D(0,∞) Weibullovo rozdělení s těžištěm x∗ = 5 a hustotami rychle klesajícími k nule má pro x → +∞ neomezené skalární skóry. Paretovo rozdělení je jednoduché rozdělení s těžkým chvostem a tedy omezenými skóry v +∞. 12
7 7.1
Použití skalární skórové funkce Zobecněné momentové odhady
Skalární skóry jsou často dány jednoduchými vzorci, ’pasují’ totiž ke svému rozdělení, takže skórové momenty (16) jsou pro běžná rozdělení elementárními (ale někdy ne zcela jednoduchými) funkcemi parametrů. Odhady parametrů metodou skórových momentů (SM odhady, rovnice (17) kde S je dáno vztahem (27)) jsou M-odhady, takže jsou konzistentní a asymptoticky normální. V jednodušších případech lze skalární skór vyjádřit ve tvaru S(x; x∗ ) kde x∗ = x∗ (θ). Protože M1 = 0, odhad xˆ∗SM těžiště x∗ se určí z rovnice n X
S(xi ; xˆ∗SM ) = 0
(29)
i=1
a má asymptotický rozptyl (9), kde ψ(x; x∗ ) = S(x; x∗ ). Tak např. pro gamma rozdělení (Tabulka 4) je S(x; α, γ) =
γ (γx − α) = γ(x/x∗ − 1), α
takže z (29) plyne, že x∗SM = x¯ je AN (α/γ, α/nγ 2 ). Pozoruhodné je, že v rovnicích (16) se proměnná vyskytuje pouze uvnitř skalárního skóru. Pokud je tedy S omezená funkce, jsou odhady všech parametrů robustní. V této skutečnosti spočívá hlavní přednost metody skórových momentů oproti metodě maximální věrohodnosti. V případě, že pracujeme s modelem, ve kterém lze předpokládat výskyt odlehlých pozorování, tato pozorování jen málo zkreslí výsledky. Cenou za tuto vlastnost je menší vydatnost nových odhadů. V obr. 2 jsou vyneseny průměry ML a SM odhadů těžiště xˆ∗ na základě 5000 výběrů z Paretova rozdělení s ω = 1, t.j. (c = 1.52), v závislosti na velikosti kontaminace d. Rozdělení jsme kontaminovali tak, že k desetině hodnot výběru X100 byla přičtena konstanta d. Pro představu uvádíme průměrné hodnoty kvantilů uspořádaných výběrů při ω = 1 : q25 = 1.23, q50 √= 1.70, q75 = 2.91, q90 = 4.65. Průměrná standardní odchylka odhadů, σ/ n, byla 0.063 pro ML a 0.08 pro SM. Vydatnost SM odhadů je tedy zhruba 0.8, ale s rostoucím d roste vychýlení ML odhadů, kdežto xˆ∗SM se ustálí na určité zvýšené hladině.
13
Pareto: x*
* ML o SM
5
10
15
20
25
d
Obr. 2. Odhady těžiště pro kontaminované Pareto Pro rozdělení s neomezeným skalárním skórem (je to většina prakticky používaných rozdělení) lze jednoduše použít postupů robustní statistiky, např. zobecnit Huberovu funkci (kap. 3). Buď X = (a, b) 6= R, a < u < v < b a nechť parametrické rozdělění (rodina) Fθ má nosič X a na něm neomezený skór S(x; θ). Definujme zobecněnou Huberovou funkci vztahem v pro x > v S(x; θ) pro u < x < v ΨF (x) = u pro x < u. Robustní odhady parametrů ’ve shodě s uvažovaným modelem’ pak obecně dostaneme řešením soustavy (17) ve tvaru n
1X k ˜ k (θ) = 0 Ψ (xi ; θ) − M n i=1 F
k = 1, ..., m
˜ k (θ) je k-tý moment funkce Ψk (x; θ). Kromě nejednodušších případů, kde M F praktická aplikace metody zatím naráží na obtížný výpočet těchto momentů.
7.2
Charakteristiky rozdělení a datových souborů
Skórové momenty Mk = ES k představují nové charakteristiky rozdělení. M1 = 0 a kořen rovnice S(x) = 0 je těžiště (’střed’), ω 2 = 1/M2 je charakteristikou variability (ω je ’poloměr’), M3 jistým způsobem charakterizuje šikmost a 1/M4 špičatost rozdělení. Namísto výběrové střední hodnoty a výběrového rozptylu můžeme za charakteristiky výběru z rozdělení považovat odhady těžiště a s-variance. ˆ aω ˆ zkonstruujeme jejich odhady z (třeba Položíme-li xˆ∗ = x∗ (θ) ˆ 2 = ω 2 (θ), 14
maximálně věrohodných) odhadů parametrů. Rodina, do které patří rozdělení populace F , je obvykle známa jen přibližně, rozumné je proto vyzkoušet více modelů. Odhady parametrů různě parametrizovaných rodin se porovnávají jen obtížně, je však snadné porovnávat odhady těžiště a s-variance (a jejich přesnost). Z tohoto hlediska nemusí být odhady parametrů konečným cílem zpracování, ale prostředkem k sestrojení výběrového těžiště a výběrové s-variance. Podobně jako pro rozdělení gamma v kap 7.1, z rovnice (29) lze nalézt analytické výrazy pro xˆ∗ i pro některá další rozdělení, viz Tabulka 4. Výběrovým těžištěm gamma (a tedy i exponenciálního) rozdělení je aritmetický průměr, lognormálního geometrický průměr a Paretova harmonický průměr. Výběrové těžiště Lomaxova rozdělení (beta 2.druhu pro p = 1) je dáno jakýmsi originálním ’průměrem’. V tabulce uvádíme i výběrové s-variance. Pro rozdělení s jedním parametrem je ω ˆ 2 ovšem funkcí xˆ∗ , pro rozdělení se dvěma parametry se po připojení rovnice pro druhý skórový moment našly analytické vzorce pro gamma a lognormální rozdělení. Obecně je třeba x ˆ∗ and ω ˆ2 hledat iteračními postupy. Tabulka 4. Výběrové těžiště a výběrová skórová variance některých rozdělení. rozdělení gamma lognormální Paretovo Lomaxovo
7.3
xˆ∗ x¯ Q x¯G = n1P xi x¯H = n/P 1/xi x¯L =
P
1 xi +1 xi xi +1
ω ˆ2 P 1 (xi − x¯)2 nP 1 log2 xi /¯ xG n (2¯ xH − 1)(¯ xH − 1)2 x¯2L (2¯ xL + 1)
Vzdálenosti
Z předchozího textu kapitoly je patrné, že za relevantní charakteristiky rozdělení lze považovat skórové momenty a za charakteristiky dat výběrové skórové momenty, obecně konstruované z odhadnutých parametrů, které jsou určeny z rovnic, ve kterých vystupují pozorovaná data Xn pouze prostřednictvím skalárního skóru. Je tedy možné namísto náhodných veličin X1 , ..., Xn s identickým rozdělením F uvažovat náhodné veličiny S(X1 ), ..., S(Xn ), kde S je skalární skór rozdělení F , a namísto pozorovaných dat uvažovat ’latentní data’ S(x1 ; θ), ..., S(xn ; θ), 15
(30)
kde S(x, θ) je skalární skór příslušný uvažovanému modelu Fθ 9podobným způsobem se můžeme dívat na data v modelech se skórem pro parametr polohy, a podle (26) i v modelech z třídy DX )). Z tohoto pohledu je přirozená vzdálenost (diference) bodů x1 , x2 ve výběrovém prostoru X rozdělení F se skalárním skórem S dána vztahem D(x1 , x2 ) = |S(x2 ) − S(x1 )|
(31)
a ’latentní’ vzdálenost ve výběrovém (datovém) prostoru Dθ (x1 , x2 ) = |S(x2 ; θ) − S(x1 ; θ)|. Po odhadnutí θ lze diferenci Dθˆ(ˆ x∗ , x∗ ) využít ke konstrukci intervalů spolehlivosti pro odhad těžiště, viz [6].
7.4
Informace a neurčitost
Domníváme se, že veličina ES 2 (Fisherova informace) představuje informaci rozdělení. ES 2 je pro štíhlá rozdělení velká a malá pro placatá, což opravňuje volbu s-variance jako ω 2 = 1/ES 2 . Funkce neurčitosti náhodné veličiny X je definována v [7] jako S 2 (x) U (x) = . (32) (ES 2 )2 Střední neurčitost rozdělení je pak EU = ω 2 . Na obr. 3 jsou znázorněny funkce neurčitosti tří rozdělení s nosičem (0, ∞). 4 3.5
U(x)
3 2.5 2 1.5 1 0.5 0
0
2
4
6
8
x
Obr. 3. Funkce neurčitosti rozdělení gamma (plná čára), lognormalného (čárkovaná) and log-logistického (beta 2. druhu s p = q = 1, tečkovaná).
16
7.5
Korelace a regrese
Máme-li ’data’ typu (30) pro dvě náhodné veličiny, můžeme je použít pro studium jejich vzájemného vztahu. Skórový kovarianční koeficient náhodných veličin X, Y s rozděleními FX , FY a skalárními skóry SX , SY definujeme jako CovS (X, Y ) = ESX SY . Skórový korelační koeficient je pak p R(X, Y ) = CovS (X, Y )/ M2 (X)M2 (Y ). (33) a jeho výběrový protějšek je Pn SX (xi , θˆX )SY (xi ; θˆY ) ˆ RXY = qP i=1 , Pn n 2 2 ˆ ˆ S (x , θ ) S (y , θ ) i X i Y i=1 X i=1 Y
(34)
kde θˆX a θˆY jsou odhady parametrů (marginálních) rozdělení X a Y . V horní části obr. 4 je znázorněna závislost průměrného Pearsonova, Spearmanova a skórového korelačního koeficientu z 5000 výběrů délky n = 100 náhodných veličin X a Y , kde X a Z jsou generovány z Paretova rozKorel. koeficient, Pareto
* ML o SM x Spearman
0.5
1
1.5
ω
2
1.5
ω
2
Chyba odhadu
0.5
1
Obr. 4. Korelační koeficienty pro kontaminované výběry z Paretova rozdělení 17
dělení a Y = 0.132X + 0.868Z, na rostoucí variabilitě rozdělení popsané veličinou ω. Teoretický korelační koeficient je rXY = 0.15. Z horního obˆ XY roste (v tomto případě) podobně jako Spearmanův rázku je patrné, že R koeficient. Pearsonův rˆXY sice vypadá na horním obrázku méně ovlivněn generovanými velkými hodnotami při rostoucím ω, ale ze spodního obrázku, kde je vynesena průměrná středněkvadratická odchylka odhadu, je jasné, že v jednotlivých případech rˆXY ’poskakuje’ zcela libovolně. Předpokládejme teď, že náhodnou veličinu Y lze vyjádřit jako Y = α0 + α1 X + ε, kde ε je náhodná proměnná s rozdělením Fε . Pro odhad koeficientů α0 a α1 požadujme minimální Fisherovu informaci reziduí, t.j. n
1X 2 S (εi ) = min ., n i=1 ε
(35)
Zde εi = yi − (α0 + α1 xi ) jsou rezidua a Sε je skalární skór rozdělení Fε . Ve dvojicích (xi , yi ) vyznačených kroužky na obr. 5 je xi realizací rovnoměrně rozdělené X a yi = −2.3 + 1.2xi + (ε − ε∗ ), kde ε je generována jako náhodná veličina s nesymetrickým Lomaxovým rozdělením s těžištěm ε∗ . Předpokládáme-li tento model reziduí, určíme během iteračního procesu jeho parametry a výsledkem ’skórové regrese’ je přímka, která se od přímky proložené standardní metodou nejmenších čtverců i od přímky proložené robustní metodou s obvyklým předpokladem symetrie rozdělení chyb (kód robustfit v Matlabu) docela výrazně liší. Přímka proložená podle (35) bere nesymetrii v úvahu. Lineární regrese
..... klasická −−−− Huberova ___ beta−prime 5
10
15
Obr. 5. Porovnání tří typů lineární regrese 18
Literatura [1] Antoch J., Vorlíčková D. (1992). Vybrané metody statistické analýzy dat. Academia. [2] Casella, G., Berger, R.L. (2002). Statistical inference. Duxbury. [3] Fabián, Z. (2001). Induced cores and their use in robust parametric estimation. Comm. Statist. Theory Methods 30, 537-556. [4] Fabián, Z. (2007). Estimation of simple characteristics of samples from skewed and heavy-tailed distribution. In Skiadas, C. (ed.) Recent Advances in Stochastic Modeling and Data Analysis, Singapore, World Scientific, 43–50. [5] Fabián, Z. (2009). The t-information and its use in multivariate problems and time series analysis. J. Statist. Planning and Inference 139, 3773-3778. [6] Fabián, Z. (2009). Confidence intervals for a new characteristic of central tendency of distributions. Comm. Statist. Theory Methods 38, 1804 - 1814. [7] Fabián, Z. (2009). O rozděleních s těžkými chvosty. Informační bulletin ČSS 20, 3. [8] Fabián, Z. (2010). Score moment estimators. Proc. of COMPSTAT 2010. [9] Fabián, Z. (2010). Characteristics of data from skewed distributions. In V. Snášel, V. Voženílek (eds.): Artificial inteligence in GIS, 11-22. [10] Fabián, Z. (2011). A new statistical tool: Scalar score function. Computer Technology and Application 2, 109-119. [11]. Huber, P. J., Ronchetti, E. M. (2003). Robust Statistic. The Approach Based on Influence Functions, Wiley, New York. [12] Johnson, N.L. (1949). Systems of frequency curves generated by methods of translations. Biometrika 36, 149-176. [13] Jurečková, J. (2001). Robustní statistické metody. Nakl. Karolinum. [14] Stehlík, M., Potocký, R., Waldl, H., Fabián, Z. (2010). On the favorable estimation for fitting heavy tailed data. Comput. Stat. 25, 485-503.
19
