DETEKCE LINEÁRNÍHO TRENDU V ROZPTYLU NORMÁLNÍHO ROZDĚLENÍ

ROBUST’2004

c JČMF 2004

DETEKCE LINEÁRNÍHO TRENDU V ROZPTYLU NORMÁLNÍHO ROZDĚLENÍ Luboš Prchal Klíčová slova: Detekce změny v rozptylu, regrese v L1 a L2 normě, radioaktivní záření. Abstrakt: Tento příspěvek je věnován detekci a odhadu dvou neznámých bodů změny, mezi nimiž se rozptyl nezávislých normálně rozdělených náhodných veličin lineárně mění z jedné konstantní úrovně na druhou. V první části je navržena vhodná testová statistika, v druhé části se pak věnujeme porovnání L1 a L2 odhadů bodů změny a parametrů rozptylu. Postup je ilustrován na reálné analýze variability vertikálních profilů radioaktivního záření.

1

Úvod

Studium problematiky detekce a odhadu měnícího se rozptylu normálně rozdělených náhodných veličin bylo motivováno praktickou potřebou analyzovat variabilitu měření vertikálních profilů radioaktivního záření. Tato data statistické veřejnosti představil na konferenci Robust’98 Hlubinka [2]. Připomeňme, že data se měří pomocí meteorologických balónů vypouštěných ze stanici v Praze-Libuši a stoupajících do výšky kolem 35-ti km, přičemž výsledkem měření jsou dvojice (xi , yi ), i = 1, . . . , n, představující průměrnou intenzitu záření yi v nadmořské výšce xi . Typický vertikální profil gama radiace je znázorněn na obrázku 1.

Obrázek 1: Typický průběh průměrného počtu gama částic v závislosti na nadmořské výšce. Hlubinka [2] podrobně diskutuje jak parametrický tak neparametrický přístup k modelování „trenduÿ pomocí regresního modelu Y = m(X)+η, kde X a Y jsou náhodné veličiny odpovídající nadmořské výšce, resp. intenzitě radiace, m(·) představuje průměrnou radiaci a η náhodnou složku měření.

316

Luboš Prchal

Nedostatky navrhovaného parametrického modelu založeného na derivaci tzv. Richardsovy růstové křivky jsou pak odstraněny jeho rozšířením podrobně popsaným v práci [3]. V tomto příspěvku se zaměříme na evidentně se měnící variabilitu měření radioaktivního záření. Odhad rozptylu radiace v závislosti na výšce
2 σi2 = var [Yi | X = xi ] založíme na čtvercích reziduí ̺2i = Yi − m(X b i ) . Jejich průběh proložený jádrovým odhadem
Pn x−xi i=1 yi K h
gK (x) = Pn (1) x−xi i=1 K h

s normálním jádrem a vyhlazovacím parametrem h zvoleným pomocí křížového ověřování (cross validation) je znázorněn na obrázku 2. Poznamenejme, že budeme-li dále mluvit o jádrové regresi, pak budeme mít na mysli neparametrický jádrový odhad (1).

Obrázek 2: Typický průběh čtverců reziduí ̺2i proložený jádrovou regresí s normálním jádrem a vyhlazovacím parametrem h = 2, 74. Pro parametrický popis chování reziduí ̺2i se jako vhodný jeví lineární model ve tvaru ̺ = Dβ + ε, kde ̺ = (̺21 , . . . , ̺2n )′ , β = (σ 2 , δ 2 )′ je vektor
neznámých reálných parametrů, σ, δ > 0, a regresní matice D = 1 d n×2 je dána vektorem samých jedniček 1n×1 a vektorem d = (d1 , . . . , dn )′ definovaným předpisem di = 0, xi − xs = , xt − xs = 1,

i = 1, . . . , s, i = s + 1, . . . , t, i = t + 1, . . . , n.

Uvědomme si, že β obecně závisí na neznámých bodech změny xs , xt , a dodejme, že s pomocí uvedeného regresního modelu odhadneme podmíněný b kde d′ představuje i−tý řádek matice D a β b je rozptyl σi2 jako σ bi2 = d′i β, i „vhodnýÿ odhad neznámých parametrů.

Detekce lineárního trendu v rozptylu normálního rozdělení

2

317

Testování modelu

Podívejme se nyní, jak otestovat hypotézu o konstantním rozptylu proti alternativě, že existují dva body změny, v nichž se charakter rozptylu mění v duchu výše popsaného regresního modelu. Než se dostaneme k samotnému testování, dodejme, že v této části budeme předpokládat normální rozdělení a nezávislost jednotlivých měření Yi . Dále předpokládejme, že „sousední měření mají stejný rozptylÿ, přesněji, že 2 2 σ2j−1 = σ2j , j = 1, . . . , n e, kde n e = ⌊n/2⌋ je spodní celá část n/2. Uvědomme si, že z tohoto předpokladu přirozeně vyplývá, že body změny s a t jsou sudá čísla. Pracujeme tedy s náhodnými veličinami
2 Y2j−1 ∼ N m(x2j−1 ), σ2j−1 , j = 1, . . . , n e,
2 Y2j ∼ N m(x2j ), σ2j−1 , j = 1, . . . , n e, a chceme testovat hypotézu o konstantnosti jejich rozptylu, tj. H1 : σi2 = σ 2 ,

proti alternativě

∀xi ∈ {x1 , x2 , . . . , xn },

A1 : ∃ xs , xt ∈ {x1 , x2 , . . . , xn }, σi2

2

=σ ,

xi − xs 2 = σ2 + δ , xt − xs = σ2 + δ2 ,

xs < xt , xi ∈ {x1 , . . . , xs }, xi ∈ {xs+1 , . . . , xt }, xi ∈ {xt+1 , . . . , xn }.

Uvažme, že díky normálnímu rozdělení veličin Yi mají náhodné veličiny
2
2 Y2j−1 − m(x2j−1 ) + Y2j − m(x2j ) , j = 1, . . . , n e, Vj = 2

2 exponenciální rozdělení s parametry ϑj = σ2j−1 . Nahraďme proto nelineární regresní model m(x) jeho odhadem m(x) b a neznámý podmíněný roz2 ptyl σ2j−1 jeho parametrickým odhadem založeným na lineárním modelu 2 σ b2j−1 = β ′ d2j−1 , j = 1, . . . , n e. Definujme dále náhodné veličiny


2
2 ̺22j−1 + ̺22j Y2j−1 − m(x b 2j−1 ) + Y2j − m(x b 2j ) Wj = = , j = 1, . . . , n e, 2 2

a předpokládejme, že i veličiny Wj mají díky normálnímu rozdělení Yi a odhadu m(x) metodou nejmenších čtverců exponenciální rozdělení, tentokrát 2 s parametry θj = σ b2j−1 = β ′ d2j−1 . Pomocí právě popsané transformace jsme nejen přešli od normálně rozdělených Yi k výběru Wj ∼ Exp(θj ), j = 1, . . . , n e, ale současně naši hypotézu H1 , resp. alternativu A1 , můžeme ekvivalentně přepsat jako hypotézu o konstantní hodnotě parametru exponenciálního rozdělení

318

Luboš Prchal H2 : θj = σ 2 ,

proti alternativě A2 : ∃ se, e t ∈ {1, 2, . . . , n e},

j = 1, . . . , n e,

2

θj = σ ,

= σ2 +

x2j−1 − x2es 2 δ , x2et − x2es

= σ2 + δ2 ,

1 ≤ se < e t≤n e, j = 1, . . . , se,

j = se + 1, . . . , e t,

j=e t + 1, . . . , n e,

kde se = ⌊s/2⌋ a e t = ⌊t/2⌋. Odvození testové statistiky T pro úlohu testování hypotézy H2 proti alternativě A2 vychází z [4] a [1] a je podrobně popsáno v [3]. Na tomto místě jen uveďme, že ji můžeme vyjádřit vztahem Pne j=2 γj Wj T = , Pne Γ · j=1 Wj

přičemž normovací konstanta Γ a konstanty γj mají v tomto konkrétním případě tvar Γ=

n e

1X γj n e j=1

a

γj =

n e j−1

(j − 1)(j − 2) X X x2j−1 − x2es + . 2 x2et − x2es e s=1 t=j e

Lze ukázat, že testová statistika T má za platnosti hypotézy H2 a při n e → ∞ asymptoticky normální rozdělení, a tudíž T −ET U= √ var T

má asymptoticky normované normální rozdělení N(0, 1), přičemž E T = 1 a rozptyl var T lze vyjádřit vztahem   Pne Pne 2 2 γ n e 1 n e   j=2 j j=2 γj var T = . 1 + P 2  − 1 = P 2 − n e+1 n e+1 n e+1 n e n e j=2 γj j=2 γj

Ze simulací ilustrujících rychlost konvergence rozdělení statistiky U k normálnímu rozdělení vyplývá, že asymptotických vlastností lze využít již při n = 50, podrobněji viz [3]. V tom případě hypotézu H2 , resp. H1 , zamítáme na hladině α ve prospěch alternativy A2 , resp. A1 , jestliže U > u(1 − α), kde u(α) je 100α%-ní kvantil normovaného normálního rozdělení. Připomeňme, že testujeme proti jednostranné alternativě „zvětšeníÿ variability o δ 2 > 0, a proto uvažujeme pouze „horníÿ kvantil u(1 − α). Při analýze variability radiace máme k dispozici výběry s rozsahy n ≈ 550, můžeme tedy bez obav užít asymptotického rozhodovacího pravidla, přičemž dle očekávání hypotézu H1 na hladině α = 0, 05 jednoznačně zamítáme pro všechna pozorování beta i gama částic.

Detekce lineárního trendu v rozptylu normálního rozdělení

3

319

Odhad modelu

V předcházejících odstavcích jsme ukázali jak otestovat „adekvátnostÿ uvažovaného regresního modelu ̺ = Dβ + ε a zbývá nám tedy odhadnout jeho neznámé parametry; body změny s a t a složky rozptylu σ 2 a δ 2 . Odhad parametrů provedeme ve dvou krocích. Nejprve pro pevné hodnoty s a t, 1 ≤ s < t ≤ n, odhadneme parametry β(s, t) jako b t) = arg min β(s, β∈R2

n X i=1



Ψ ̺2i − d′i β = arg min RS β(s, t) , β∈R2

kde Ψ je vhodně zvolená funkce. Ve druhém kroku odhadneme
body změny s a t tak, abychom minimalizovali ztrátovou funkci RS β(s, t) , tedy 
b t) . sb, b t = arg min RS β(s, s=1,...,n−1 t=s+1,...,n


b sb, b Tím také dostaneme výsledný odhad β pomocí β t . Parametry lineárního modelu většinou odhadujeme metodou nejmenších čtverců (dále jen L2 regrese). V našem případě však nemáme splněn jeden ze základních předpokladů klasického lineárního modelu, a sice homoskedasticitu náhodné složky ε. „Neblahý vlivÿ heteroskedasticity náhodné složky na odhad parametrů β lze omezit užitím metody vážených nejmenších čtverců (WLS) s diagonální maticí vah Wn×n tvořenou prvky wii = 1/b τi2 , kde τbi2 je odhad rozptylu var εi = τi2 . Jako vhodný, „nezávislýÿ na metodě nejmenších čtverců se nabízí odhad pomocí již zmíněné jádrové regrese (1) ve tvaru
2 τbi2 = ̺2i − gbK (xi ) . S využitím informace o variabilitě náhodné složky ε dostáváme odhad neznámých parametrů β v podobě b β = (D ′ WD)−1 D ′ W̺. (2) WLS

Jako robustní alternativu k metodám nejmenších čtverců uveďme regresi v L1 normě (dále jen L1 regrese) odpovídající minimalizační úloze min2

β∈R

n X i=1

|̺2i − d′i β|.

(3)

Jedním z přístupů k řešení L1 regrese je numerická metoda iteračně vážených nejmenších čtverců (IWLS). Minimalizace úlohy (3) pomocí metody IWLS odpovídá řešení soustavy D′ W (β)̺ = D ′ W (β)Dβ, vzhledem k neznámým parametrům β ∈ R2 , přičemž matice vah W (β)n×n je diagonální s prvky wii (β) definovanými předpisem
sgn ̺2i − d′i β wii (β) = , ̺2i − d′i β 6= 0, (4) ̺2i − d′i β = 0,

̺2i − d′i β = 0.

320

Luboš Prchal

Jelikož váhy na rozdíl od metody WLS tentokrát závisejí na neznámých parametrech β, není možné pro odhad β užít přímo vztah (2), nýbrž je třeba přikročit k numerickému řešení. Metoda IWLS vychází z počátečního od(0) hadu βL1 , který v jednotlivých iteracích postupně „vylepšujeÿ předpisem (l+1)

β L1

 −1 (l) (l) = D ′ W (β L1 )D D′ W (βL1 )̺

(5) (l)

až do splnění vhodného zastavovacího pravidla. Dodejme, že β L1 značí L1 od(l)
had parametrů β po l iteracích dle vztahu (5) a W β L1 jsou známé váhy (l)

dány předpisem (4) pro již spočtenou hodnotu β L1 . V jednotlivých iteracích (l)

(l+1)

tedy známe matici vah W (β L1 ), a proto následující odhad parametrů βL1 získáme stejně jako u metody WLS vztahem (2). Druhým přístupem vedoucím k nalezení optimálního řešení problému L1 regrese, pokud takové řešení existuje, je přeformulovat regresní úlohu (3) jako standardní minimalizační úlohu lineárního programování ve tvaru n X
− min ǫ+ i + ǫi ǫ+, ǫ− i=1

za podmínek − 2 d′i β + ǫ+ i − ǫ i = ̺i ,

− ǫ+ i , ǫi ≥ 0,

i = 1, . . . , n,

σ 2 , δ 2 ≥ 0.

K jejímu vyřešení pak lze užít standardních nástrojů obsažených v matematických a statistických programech, např. funkci linprog implementovanou v Optimization Toolbox programového vybavení Matlab, Release13. Dodejme, že zde prezentované výsledky jsou získány metodu IWLS s tím, že výrazně rychlejší numerické řešení se jen nepatrně liší od přesného řešení úlohy lineárního programování. Podrobněji je volba metody diskutována v práci [3].

4

Porovnání metod

Podívejme se nyní na získané odhady prezentovanými metodami jednak z pohledu odhadu složek variability σ 2 a δ 2 , jednak z pohledu bodů, ve kterých dochází k jejich změnám. O odhadech bodů změny se obecně dá říci, že použití L2 regrese vede k odhadům pouze jednoho bodu změny, tedy na model se skokovou změnou v chování rozptylu. Tento typický rys L2 regrese je způsoben značnou citlivostí L2 odhadů na velké, „odlehléÿ hodnoty reziduí. Naproti tomu, užitím L1 ztrátové funkce získáme „očekávánýÿ odhad s postupným lineárním růstem variability σ 2 (x). Rozdíl v chování L2 a L1 odhadů ilustrujme na pozorování gama částic z 10. října 1995. Průběh čtverců reziduí ̺2i proložený L2 i L1 odhadem variability je znázorněn na obrázku 3 (levý graf). Na tomtéž obrázku 3 jsou kolečkem vyznačena dvě rezidua ve výšce asi 15 km, která způsobují skok

Detekce lineárního trendu v rozptylu normálního rozdělení

321

Obrázek 3: Čtverce reziduí proložené odhadnutým rozptylem metodou WLS (čárkovaně) a L1 regresí (plná čára). Levý graf znázorňuje odhady získané ze všech reziduí, pravý graf odhady po vynechání dvou zakroužkovaných „odlehlýchÿ reziduí. v L2 odhadu. Na pravém grafu obrázku 3, který odpovídá stejnému pozorování, vynecháme-li dvě vyznačená rezidua, vidíme, že nově odhadnuté body změny metodou WLS se „přiblížilyÿ nezměněnému L1 odhadu. Ačkoli nemáme a priori žádnou informaci o chování variability měření, zdá se rozumné přiklonit se k robustnějším L1 odhadům, vesměs podporujícím myšlenku lineární změny mezi dvěma konstantními hladinami rozptylu. Odhady bodů změny xs a xt pomocí L1 regrese více odpovídají také naší původní představě o průběhu rozptylu založené na neparametrické jádrové regresi (1). Připomeňme, že parametry σ 2 a δ 2 si můžeme představit jako rozptyl měření ve výškách do xs , resp. nad xt . Uvážíme-li dále, že jsme při testování předpokládali normální rozdělení dat, pak odhad neznámých parametrů σ 2 a δ 2 regresí v L1 normě se zdá být nevhodný. Pro odhad samotných složek rozptylu σ 2 a δ 2 bychom spíše měli volit metodu nejmenších čtverců, resp. její váženou variantu WLS. Ve světle předcházejících úvah se jako optimální metoda pro odhad chování variability měření radiace jeví kombinace L1 a L2 přístupu. Počítejme proto nejprve L1 odhad bodů změny obvyklým dvoukrokovým postupem, tj. v prvním kroku odhadněme pro pevné body změny s a t, 1 ≤ s < t ≤ n, složky rozptylu vztahem n X
b βL1 (s, t) = arg min RSL1 β(s, t) = arg min |̺2i − d′i β|, β∈R2

β∈R2

i=1


a na základě ztrátové funkce RSL1 β(s, t) pak v druhém kroku odhadněme body změny jako
b (s, t) . {s∗, t∗ } = arg min RSL1 β L1 s=1,...,n−1 t=s+1,...,n

322

Luboš Prchal

Složky rozptylu následně odhadněme metodou vážených nejmenších čtverců ′ ∗ ∗ −1 ′ b β∗ = β Ds∗t∗ W̺, WLS (s , t ) = (Ds∗t∗ WDs∗t∗ )

kde Ds∗t∗ je regresní matice odpovídající pevným bodů změny s = s∗ a t = t∗ . Matice vah W je diagonální s prvky wii = 1/b τi2 , i = 1, . . . , n, přičemž τbi2 je
2 2 odhad rozptylu ̺i získaný pomocí jádrové regrese τbi2 = ̺2i − gbK (xi ) . Představená kombinace L1 a L2 metod si zachovává výhody robustního odhadu bodů změny, přičemž oproti samotné L1 metodě „věrnějiÿ odhaduje složky rozptylu. Odhad průběhu rozptylu měření získaný touto kombinovanou metodou je znázorněn na obrázku 4.

Obrázek 4: Optimální odhad variability měření kombinovanou metodou.

Reference [1] Gupta A.K., Ramanayake A. (2001). Change points with linear trend for the exponential distribution. J. Statist. Plann. Inference 93, 181 – 195. [2] Hlubinka D. (1998). Metody pro prokládání křivek s použitím na reálných datech. In ROBUST’98 (Antoch J. a Dohnal G., eds.), JČMF, Praha, 55 – 75. [3] Prchal L. (2004). Neparametrické odhady pro analýzu funkcionálních dat. Diplomová práce, MFF UK, Praha. [4] Worsley K.J. (1986). Confidence regions and test for a change point in a sequence of exponential family random variables. Biometrika 73, 91 – 104. Poděkování: Autor děkuje prof. RNDr. Jaromíru Antochovi, CSc., za jeho nezištnou pomoc a neocenitelné rady v průběhu vzniku tohoto článku. Práce vznikla s podporou grantů GAČR 201/03/0945 a MSM 113200008. Adresa: L. Prchal, KPMS MFF UK, Sokolovská 83, 186 75 Praha 8 E-mail : [email protected]