XXVI. ASR '2001 Seminar, Instruments and Control, Ostrava, April 26 - 27, 2001
Paper 44
Odhad optimálního stupně regresního polynomu MORÁVKA, Jan Ing., Ph.D., Město, 739 61,
Třinecký inženýring, a.s., Středisko projekce, Frýdecká 126, Třinec–Staré
[email protected],
http://www.inzenyr.trz.cz
Abstrakt: V praxi je často třeba najít regresní závislost ve tvaru polynomu. V příspěvku je uveden přehled a otestování vhodných kritérií odhadu optimálního stupně regresního polynomu Byly testovaná informační kritéria autorů: Akaike, Schwarz-Rissanen, HannanQuinn, Geweke-Meese a další numericko-statistická kritéria: kRSC, D_dif, F-test, DW (Durbin-Watson), MSE, MEP, Rp2, W (Wald), znaménkový test, standardizovaná asymetrie a exces, JB (Jarque-Berra), BP (Breusch-Pagan), CW (Cook-Weisberg), RSC = RSS (Residual Sum of Squares) a sk2 (reziduální rozptyl). Klíčová slova: regrese, regrese polynomická, optimální stupeň regresního polynomu
1 Úvod V praxi je třeba poměrně často najít regresní závislost, která má na základě fyzikálněchemické analýzy, či na základě expertního odhadu, tvar regresního polynomu (ei ~ N(0, σe)) s
y i = b0 + b1 x i + ... + bs x i + ei = y mi + ei ,
i = 1,..., n ,
(1)
avšak stupeň s regresního polynomu není znám. Prvním úkolem je proto zjistit odhad parametru s. Odhady regresních koeficientů bj (j = 0,...,s) lze pak získat standardními postupy linearizované regrese pomocí metody nejmenších čtverců (MNČ).
2 Metody řešení a kritéria hodnocení A. Zdánlivě logický postup přidávání statisticky významných členů s vyšší mocninou [ECKSCHLAGER, K. AJ. 1980], [LEPŠ, J. 1996] (nebo opačně jejich ubírání při startování od maximálně možného stupně) až k statisticky nevýznamnému koeficientu určité mocniny, nevede k uspokojivým výsledkům [ANDĚL, J. 1993], zvláště z důvodu přirozené multikolinearity polynomických modelů, která způsobuje nevýznamnost regresních koeficientů na základě t-testů [MELOUN, M. & MILITKÝ, J. 1994]. B. Hledání optimálního stupně regresního polynomu podle reziduálního součtu čtverců (RSC) také nevede k cíli, protože RSC monotónně klesá (index determinace R2 monotónně stoupá) se zvyšujícím se stupněm polynomu [MELOUN, M. & MILITKÝ, J. 1994]. C. Z literatury – např. [RALSTON, A. 1973], [ANDĚL, J. 1993], [CIPRA, T. 1986] je známo, že při dosažení skutečného (správného) stupně regresního polynomu dochází k ustálení hodnot reziduálního rozptylu. Matematicky lze tuto formulaci vyjádřit následovně: model (1) obsahuje k = s + 1 regresních parametrů b0,...,bs. Označme s2k výběrový reziduální rozptyl. Předpokládejme, že skutečný počet regresních parametrů je k0 (takže skutečný stupeň regresního polynomu je s0 = k0 – 1). Dá se ukázat, že platí: -1-
k < k 0 pro k ≥ k 0
2 2 s k > σ e . 2 2 s k ≈ σ e
(2)
Otázkou však dlouho zůstávalo, jak z grafu hodnot s2k vyčíst právě tu hodnotu k0, od níž počínaje již graf dostává přibližně vodorovný (konstantní) charakter [RALSTON, A. 1973], [ECKSCHLAGER, K. AJ. 1980]. Nakonec se ukázalo, že je třeba zavést třídu vhodných funkcí penalizujících počet členů polynomu. Tyto funkce nabývají svého globálního extrému (minima) v případě správné volby stupně regresního polynomu. Mezi nalezená optimalizační kritéria hodnocení správného stupně polynomu patří dle [ANDĚL, J. 1993], [ARLT, J. 1999], [CIPRA, T. 1986], [MELOUN, M. & MILITKÝ, J. 1994] následující: 2
AIC k = ln(s k ) +
2k , n
(3)
což je známé tzv. Akaikeovo informační kritérium (Akaike’s Information Criterion, používané např. v programech ADSTAT, QC Expert, TSP, JMP IN a Matlab), kde symbol ln ve vztahu označuje přirozený (naturalis) logaritmus. Toto kritérium však občas nadhodnocuje odhad stupně polynomu, tj. odhad může s určitou kladnou pravděpodobností konvergovat i k nějaké vyšší hodnotě, než je k0. Modifikované AIC dle Ozakiho je uvedeno např. v [ARLT, J. 1999], [CIPRA, T. 1986]. Dalším oblíbeným kritériem, které navrhli Schwarz a Rissanen je 2
SRk = ln(s k ) +
k ln(n) . n
(4)
Používá se také kritérium, které odvodili Hannan a Quinn (s volbou c > 1, c = 2, 3) 2
HQk = ln(s k ) +
2kc ln(ln(n)) . n
(5)
V praxi se osvědčilo i kritérium dle autorů Geweke a Meese 2
GM k = s k (1 +
k ). 4 n
(6)
Obecně jsou výše uvedená (a další) kritéria používána hlavně v analýze časových řad pro stanovení optimálních ARMA, ARIMA modelů, či ve vícerozměrné regresní analýze. D. Vzhledem k monotónnímu klesání RSC při zvyšování stupně polynomu je možné zkusit použít také autorem příspěvku navržené jednoduché kritérium obsahující součin počtu členů polynomu k (které lineárně stoupají) a RSC s předpokladem, že u tohoto kritéria by mohlo docházet k ostrému lokálnímu minimu při volbě optimálního stupně regresního polynomu: n
kRSC k = k .RSC = k ∑ eik . 2
(7)
i =1
-2-
Obdobná kritéria obsahující součin k2 a RSC, či k, k2 a s2k byly také autorem testovány, ovšem jejich výsledky nevykazovaly tak dobré vlastnosti, jako má kritérium kRSC. E. Při analýze časových řad (zejména v ekonometrii) se k určení stupně trendového polynomu, nebo k určení řádu diferencování nestacionárního modelu ARIMA, používá metoda postupných diferencí [CIPRA, T. 1986]. Zde se vychází z faktu, že při postupném diferencování hodnoty odhadnutých rozptylů vysvětlované proměnné y klesají až k minimu (kdy je dosažena stacionarita), a pak opět začnou růst. Sledovaným kritériem je tedy: D _ dif d = σˆ 2 ∆d y ≈ s 2 ∆d y = cov(diff ( y, d )) .
(8)
kde d = k – 1 ∈ (0, ..., kmax – 1) je řád diferencování proměnné y, s2 označuje výběrový rozptyl a poslední tvar kritéria je zapsán v notaci programu Matlab. F. V moderních statistických postupech a programech regresní analýzy se doporučují a používají ještě další kritéria [MELOUN, M. & MILITKÝ, J. 1994]: • hodnota statistiky Fisherova F-testu (v základním i upraveném tvaru), • kritérium MEP (Mean Error Prediction – střední kvadratická chyba predikce), • predikovaný koeficient determinace Rp2. Pozn.: Scottovo kritérium přeurčenosti a tím i multikolinearity regresního modelu nemá význam použít z důvodu přirozené (nezbytné) multikolinearity polynomických modelů.
3 Analyzovaný model a data V literatuře [ANDĚL, J. 1993] je analyzován případ regresního modelu 3. řádu pro jednu hodnotu směrodatné odchylky aditivního šumu ei ~ N(0, σe), σe = 1/2 a jsou zde použita jenom kritéria AIC, SR, HQ a GM včetně uvedení hodnot s2k. V příspěvku jsou uvedeny výsledky pro více hodnot směrodatné odchylky aditivního šumu a pro další diskutovaná kritéria. V programu Matlab byl simulován polynomický regresní model 3.stupně tvaru 2
3
y i = 10 + 2 x i + 2 x i − x i + ei = y mi + ei ,
i = 1 ... 30,
ei ~ N (0, σ e ), σ e ∈ {0, 0.01, 0.1, 0.2, 0.5, 1, 2, 5},
,
(9)
x i = 0.1(i − 1), s uvažováním hodnoty násady generátoru pseudonáhodných posloupností (randn s normálním rozdělením) seed = 1357, maximálního stupně polynomu kmax = 9 a hodnotou c = 3 kritéria HQ. Na obr.1 je pro názornost uveden graf výstupu (připomínající tvarem momentovou charakteristiku asynchronního motoru) samotného modelu ym spolu s průběhem výstupu y zatíženého aditivním výstupním šumem e se směrodatnou odchylkou σe = 1/2.
-3-
Model polynomu 3.stupne 15 14
y __ ym _ _
13 12 11 10 smerodatná odchylka sumu = 0.5 9 8 7 0
0.5
1
1.5 x
2
2.5
3
Obr. 1. Model analyzovaného polynomu 3.stupně se σe = 1/2 Na obrázku je vidět, že tzv. oscilace funkce ym_osc samotného modelu (rozdíl ym_max – ym_min) v daném intervalu hodnot nezávisle proměnné x má hodnotu asi 6. Prakticky to znamená, že šum se zvětšující se směrodatnou odchylkou tuto funkci postupně „překryje”, bude tedy docházet ke snižování odhadnutého stupně polynomu a při hodnotě srovnatelné s hodnotou oscilace dojde k degeneraci polynomu. Zde budou zřejmě kritéria signalizovat konstantní průběh funkce, tj. polynom 0.stupně (s = 0, což znamená konstantu) s jedním absolutním členem (k = 1) – viz obr.2: Model polynomu 3.stupne 25
y __ ym _ _
20
15 smerodatná odchylka sumu = 5 pocet hodnot n = 30
10
5
0
0
0.5
1
1.5 x
2
2.5
3
Obr. 2. Model analyzovaného polynomu 3.stupně se σe = 5
4 Vyhodnocení polynomické regrese V tab.1 jsou uvedeny souhrnné přehledné výsledky odhadů optimálního stupně regresního polynomu podle výše uvedených i níže diskutovaných kritérií v závislosti na směrodatné odchylce výstupního šumu modelu.
-4-
Tab. 1. Odhad počtu členů k (k = s+1) regresního polynomu v závislosti na σe σe → 0 0.01 0.1 0.2 0.5 1 2 5 Kritérium ↓ AIC 4 4 4 4 4 3 3 1 4 4 4 4 4 3 3 1 SR 4 4 4 4 4 3 1 1 HQ 4 4 4 4 4 3 1 1 GM RSC 4 4 9 9 9 9 9 9 kRSC 4 4 4 4 4 3 1 1 2 sk 4 4 4 4 4 6/8 6/8 8 4 4 3 3 2 2 1 1 D_dif 4 4 4 4 4 3 3 3(6) F-statistika (>10) 4 4 4 4 3 3 1 DW Pozn.: Kurzívou jsou označeny hodnoty „ustálení” kritérií (neostré, ploché optimum), normálním písmem pak ostrá optima těchto kritérií. Tučným písmem jsou zvýrazněna spolehlivá a správná kritéria, kurzívou pak kritéria použitelná s určitými omezeními. DW ... hodnota Durbin-Watsonova koeficientu. Hodnota kritéria F-testu dosahuje při optimální hodnotě stupně polynomu (na rozdíl od jiných kritérií) svého maxima a hodnoty kritéria DW zde spadají dovnitř kritických mezí nezávislosti veličiny. Ve statistických programech Statgraphics, EasyReg a QC Expert byly pro model s σe = 0.5 otestovány následující kritéria hodnocení kvality modelu: • hodnota statistiky F-testu • střední kvadratická chyba odhadu (MSE ... Mean Squared Error) • střední kvadratická chyba predikce (MEP ... Mean Error of Prediction) • predikovaný koeficient determinace Rp2 • autokorelace reziduí: Durbin-Watsonův koeficient (DW), Waldův test (W) • trend reziduí: znaménkový test • normalita reziduí: standardizovaná asymetrie a exces, Jarque-Berraův test (JB) • homoskedasticita reziduí: Breusch-Paganův (BP) a Cook-Weisbergův (CW) test. Lze však konstatovat, že očekávané identifikační vlastnosti se projevily pouze u hodnot kritéria F-test (ostré maximum), JB (ostré minimum, avšak rezidua přitom pro všechny stupně polynomu vykazovaly normalitu) a DW (kde se projevily statisticky nevýznamné hodnoty autokorelace od optimálního stupně polynomu výše) – obr.3: F -krité rium
krité rium Du rb in a -W a tso na
2.5
1.E + 06 σ e = 0.01 1.E + 05
σe = 5
2
1.E + 04 1.E + 03
1.5
1.E + 02 1
1.E + 01 1.E + 00 1.E -01
1
2
3
4
5
6
7
8
9
10
0.5
k
1.E -02
σe = 0
0 1
2
3
4
5
6
7
8
9
10 k
1.E -03
Obr. 3. Průběhy kritérií F-test a DW -5-
Tento závěr je celkem logický: při polynomech nižšího stupně se dá očekávat autokorelace reziduí avšak narušení jejich normality a homoskedasticity spíše nemá opodstatnění. Hodnoty Rp2 a MEP vykazovaly nejednoznačný, nebo nevhodný průběh (více lokálních maxim pro k = 6, 8 a 9, či minimum pro k = 8). Proto jsou v tab.1 uvedeny pouze kritéria F-test a DW. Na obr.4 a obr.5 lze vidět průběhy uvedených základních numericko-statistických a informačních kritérií pro σe = 0.5. Reziduální soucet ctvercu
k * RSC
log10 RSC
2
200 150
1.5
100 1 0.5
50 0
5 Reziduální rozptyl
0
10
0 -0.5 -1
5 Rozptyl diferencí
10
0
5 k=rád+1
10
4 log10 D_dif
log10 sk^2
0.5
0
0
5 k=stupen+1
10
2 0 -2
Obr. 4. Průběhy základních kritérií (pro σe = 1/2) AIC : Akaike
GM : Geweke-Meese
2
6
1
4
0 2
-1 -2
0
5 10 SR : Schwarz-Rissanen
0
2
2
1
1
0
0
-1
0
5 k
10
-1
0
0
5 10 HQ : Hannan-Quinn
5 k
10
Obr. 5. Průběhy informačních kritérií (pro σe = 1/2) Na základě údajů v tabulce a v grafech lze souhrnně konstatovat, že: • kritéria RSC a s2k jsou spíše nepoužitelná pro stanovení stupně regresního polynomu • kritérium F-test dává nejednoznačný výsledek pro největší směrodatnou odchylku šumu z hlediska polohy globálního maxima. Pokud však budeme uvažovat zásadně pouze první lokální maximum, pak jsou jeho výsledky jednoznačné • použití kritéria DW pro nulový šum je nesmyslné (zde kritérium reaguje na numerické chyby výpočtu), v ostatních případech poskytuje správné hodnocení • ostatní kritéria dávají dobré a robustní odhady -6-
•
podle průběhů závislostí k na σe se uvedená kritéria sdružují do čtyř skupin: {F-test}, {AIC, SR, DW}, {HQ, GM, kRSC} a {D_dif} - viz obr. 6: O d h a d p o čtu čle n ů p o lyn o m u
4 F-tes t
k = s+ 1
3 D_dif
HQ
2
A IC
1 A IC
HQ
1
2
D_dif
F-tes t
0 0
3
4
5
σ e - sm ě ro d a tn á o d ch ylka e
Obr. 6. Průběhy kritérií F-test, AIC, HQ a D_dif pro σe ∈ <0, 5> Z průběhů zobrazených skupin kritérií je zřejmé, že jejich robustnost (odolnost) vzhledem k velikosti šumu klesá v už uvedeném pořadí, tj. {F-test}, {AIC, SR, DW}, {HQ, GM, kRSC} a {D_dif}, které je nejcitlivější.
5 Závěr Z uvedených rozborů je vidět, že: 1. Pro odhad optimálního stupně regresního polynomu jsou vhodnými kritérii: • obecná informační kritéria AIC, SR, HQ, GM, • kRSC, • D_dif, s určitými omezeními lze použít i kritéria: • F-test (při větších amplitudách šumu je třeba uvažovat první maximum), • DW (nelze použít pro nulovou amplitudu šumu) 2. Správný odhad stupně regresního polynomu je možné uskutečnit pomocí uvedených kritérií (kromě příliš subtilního D_dif) pouze do velikosti aditivního šumu asi o řád nižší, než je oscilace funkce na daném intervalu, čili pro:
σ
e
≈ se ≤
y m _ osc 5 ÷ 10
=
y m _ max − y m _ min 5 ÷ 10
(10)
3. Při větších velikostech šumu dochází k podhodnocování stupně polynomu z důvodu problematické rekonstrukce (identifikace) průběhu výstupu skutečného modelu.
-7-
6 Literatura ANDĚL, J. 1993. Statistické metody. 1.vyd. Praha : Matfyzpress MFF UK Praha, 1993. 246 s. ARLT, J. 1999. Moderní metody modelování ekonomických časových řad. 1.vyd. Praha: Grada Publishing, s.r.o., 1999. 312 s. ISBN 80-7169-539-4. CIPRA, T. 1986. Analýza časových řad s aplikacemi v ekonomii. 1.vyd. Praha : SNTL/ALFA, 1986. 248 s. ECKSCHLAGER, K. AJ. 1980. Vyhodnocování analytických výsledků a metod. 1.vyd. Praha : SNTL/ALFA, 1980. 224 s. LEPŠ, J. 1996. Biostatistika. 1.vyd.-dotisk. České Budějovice : skriptum BF JU České Budějovice, 1996. 166 s. MELOUN, M. & MILITKÝ, J. 1994. Statistické zpracování experimentálních dat. 1.vyd. Praha : PLUS, 1994. 839 s. ISBN 80-85297-56-6. RALSTON, A. 1973. Základy numerické matematiky. 1.vyd. Praha : Academia, 1973. 636 s.
-8-
