基于逐步回歸的穩(wěn)健估計(jì)和異常值檢測

崔 樂, 吳 迪, 成麗波

(長春理工大學(xué) 理學(xué)院, 長春 130022)

0 引 言

近年來,異常值檢測一直是統(tǒng)計(jì)診斷中一個(gè)比較活躍的研究課題。在實(shí)際生活中,由于種種原因不能夠和其他數(shù)據(jù)一起用于多元統(tǒng)計(jì)的線性回歸模型中的值,稱為異常值。換句話說,參數(shù)估計(jì)、模型構(gòu)建及預(yù)測都會受到異常值的影響,因此,這篇文章主要研究的是基于逐步回歸的穩(wěn)健估計(jì)和異常值的檢測。

比較常見的統(tǒng)計(jì)量有殘差、學(xué)生化殘差、Cook距離以及W-K統(tǒng)計(jì)量等。殘差和學(xué)生化殘差[1-3]都可以用來檢測異常值,但是殘差沒有考慮到異方差性,而學(xué)生化殘差考慮到了這一點(diǎn),因此用學(xué)生化殘差檢測異常值比用普通殘差更有效。而Cook距離和W-K統(tǒng)計(jì)量[4-6]僅可以判斷出數(shù)據(jù)是否有強(qiáng)影響點(diǎn),但是異常點(diǎn)和強(qiáng)影響點(diǎn)之間沒有必然的關(guān)系,即異常點(diǎn)不一定是強(qiáng)影響點(diǎn),強(qiáng)影響點(diǎn)也不一定是異常點(diǎn)。

本文基于逐步回歸模型,從殘差平方和的角度出發(fā),研究出檢測異常值的方法,并估計(jì)異常值的大小。通過與傳統(tǒng)方法的比較,驗(yàn)證該方法的有效性。

1 逐步回歸

在多元線性回歸模型中,自變量的選擇實(shí)質(zhì)上就是模型的選擇?，F(xiàn)設(shè)一切可供選擇的變量是t個(gè),它們組成的回歸模型稱為全模型(記m=t+1),在獲得n組觀測數(shù)據(jù)后,有模型

其中:Y是n×1的觀測值,β是m×1未知參數(shù)向量,X是n×m結(jié)構(gòu)矩陣,并假定X的秩為m。

β=(βp,βq)′,X=(Xp?Xq)

下面從不同的角度給出自變量選擇的準(zhǔn)則[7]。

準(zhǔn)則1 平均殘差平方和達(dá)到最小

準(zhǔn)則2CP統(tǒng)計(jì)量達(dá)到最小

該準(zhǔn)則是由馬斯洛于1964年從預(yù)測的角度提出:

準(zhǔn)則3 AIC準(zhǔn)則

該準(zhǔn)則由日本統(tǒng)計(jì)學(xué)家赤池弘次于1974年根據(jù)極大似然估計(jì)原理提出:

AIC=nln(SSEP)+2P

其中:SSEm是m個(gè)自變量x1,x2,…,xm所對應(yīng)的殘差平方和;SSEP是p個(gè)自變量x1,x2,…,xp所對應(yīng)的殘差平方和。

快速選擇變量的算法有很多,例如:向前法、向后法、逐步回歸法[7]等,其中逐步回歸法是應(yīng)用最廣的一種方法。它的具體做法是先將變量一個(gè)一個(gè)的引入,當(dāng)引入到第3個(gè)變量之后的每一步,首先對已引入的變量進(jìn)行剔除。這樣,自變量將不斷的引入、剔除、再引入、再剔除……直到自變量不能被剔除,同時(shí)也無法引入自變量為止。

2 基于逐步回歸的殘差平方和

經(jīng)過逐步回歸后生成最優(yōu)多元線性回歸模型如下:

yi=β0+β1xi1+β2xi2+…+βpxip+εi,i=1,2,…,n

(1)

其中:p為解釋變量的數(shù)目;βj(j=0,1,…,p)為回歸系數(shù);εi為隨機(jī)誤差,且εi～N(0,σ2),i=1,2,…,n;β=(β0,β1,…βp)T;σ為未知參數(shù)。

均值漂移模型[8-9]是在第i個(gè)數(shù)據(jù)點(diǎn)上增加一個(gè)漂移項(xiàng)δ,即在這個(gè)數(shù)據(jù)點(diǎn)yi處的均值發(fā)生了非隨機(jī)漂移,若δ顯著不等于零,則yi處的均值發(fā)生了漂移,說明此點(diǎn)為異常點(diǎn)。

由于事先不知道在線性模型中出現(xiàn)異常值,因此,可以先假定模型中沒有異常值,其線性模型的矩陣形式為

Y=Xβ+ε

(2)

其中

SSE=YT(I-H(X))Y

殘差平方和(SSE)會受到異常值大小的影響,即隨著異常值的增加,殘差平方和(SSE)也增加。

3 異常值的檢測

基于均值漂移模型,假設(shè)模型中的第i個(gè)觀測值為異常值,即異常值的大小為δi,其余的皆為正常數(shù)據(jù),計(jì)算此時(shí)的殘差平方和(SSEi)為

(3)

(4)

(5)

4 穩(wěn)健估計(jì)和異常值的識別

4.1 標(biāo)準(zhǔn)差的穩(wěn)健估計(jì)

M估計(jì)是基于最小二乘估計(jì)發(fā)展起來的一種抗差估計(jì)(Robust Estimation)方法[10-12],是由huber于1964年最先提出來的,也稱為廣義最大似然估計(jì)。M估計(jì)已經(jīng)成為最經(jīng)典的一種穩(wěn)健估計(jì)方法。

M估計(jì)的估計(jì)方程寫成矩陣形式是這樣的:

XTWXβ=XTWY

(6)

迭代公式如下:

(7)

其中:W是以ωi,i=1,2,…,n為對角線的權(quán)矩陣;X是解釋變量矩陣,X=(x1,x2,…,xn);Y是因變量向量,Y=(y1,y2,…yn)T。

4.2 異常值的識別

5 模擬與實(shí)例分析

5.1 模擬實(shí)驗(yàn)

假設(shè)有一多元線性回歸模型y=3+5x1-4x2+4x3-3x4+x5-2x6+2x7-x8+3x9+6x10+ε,現(xiàn)利用計(jì)算機(jī)模擬產(chǎn)生100個(gè)數(shù)據(jù)樣本。下面考慮3種方案進(jìn)行實(shí)驗(yàn)。

方案1 將大小為5.5的異常值加入到第49個(gè)樣本觀測值中。

方案2 將大小為6,-5的異常值分別加入到第31和69個(gè)樣本觀測值中。

方案3 將大小為-6.5, 5.5, 4.5, -3的異常值分別加入到第48～52個(gè)樣本觀測值中。

1) 異常值檢測

對于3種方案的D統(tǒng)計(jì)量如圖1所示。圖1是基于M估計(jì)的D統(tǒng)計(jì)量,由圖1(a)可以看出,第50個(gè)觀測值對應(yīng)的D統(tǒng)計(jì)量遠(yuǎn)大于3,判定為異常值。同樣由圖1(b)可以看出,第32個(gè)和第70個(gè)觀測值為異常值;由圖1(c)或者看出,從第49～53個(gè)觀測值均被檢測出異常值。由此可以看出該檢測異常值的統(tǒng)計(jì)量可以很好地檢測出異常值,并且由于方案3是5個(gè)連續(xù)的異常值,因此說明此方法對異常值的遮蔽現(xiàn)象有一定的作用,能夠有效地檢測出連續(xù)幾個(gè)異常值。

圖1 方案1～3的基于M估計(jì)的D統(tǒng)計(jì)量Fig.1 D Statistics Based on M Estimation for Schemes 1～3

圖2是方案1～3中每個(gè)樣本觀測值的學(xué)生化殘差示意圖。按照傳統(tǒng)的方法,將觀測值的學(xué)生化殘差的絕對值大于3的認(rèn)為是異常值。從圖2(a)中可以看出,無法準(zhǔn)確地檢測出方案1中異常值,圖2(b)中,可以很好地檢測出方案2中的異常值,從圖2(c)中,學(xué)生化殘差只能檢測出第49, 50, 52, 53個(gè)觀測值為異常值。由此可以說明,用新構(gòu)建的基于穩(wěn)健估計(jì)的D統(tǒng)計(jì)量比傳統(tǒng)的學(xué)生化殘差檢測異常值更加有效。

圖2 方案1～3的學(xué)生化殘差Fig.2 Student Residual in Scheme 1～3

2) 異常值大小的估計(jì)

通過公式(4)來計(jì)算異常值的大小,表1是3種不同方案的異常值大小的評估值。從表中可以看出,對于方案1～2,本文構(gòu)建的統(tǒng)計(jì)量都可以很好地檢測出異常值的大小。對于方案3異常值大小估計(jì)精確度不如方案1～2,但仍然還是比較準(zhǔn)確的。

表1 基于方案1～3的異常值大小估計(jì)值Tab.1 Outlier size estimates based on Schemes 1～3

5.2 實(shí)例分析

在房地產(chǎn)行業(yè)中,影響房屋價(jià)值的因素有很多。本實(shí)例中,根據(jù)房屋價(jià)值的影響因素及人們的偏好,記錄的影響房屋價(jià)值的指標(biāo)包括面積、戶型、當(dāng)前樓層、總樓層、朝向、裝修情況。其中,戶型為幾室?guī)讖d,可以分開作為單獨(dú)變量;當(dāng)前樓層分為底層、中層和高層;朝向分為北、東北、東、東南、南、西南、南北、西、西北、東西;裝修情況分為毛坯、普通裝修、精裝修。

房價(jià)和房屋面積需要以10為底取對數(shù),目的是降低房價(jià)和面積的數(shù)量級,提高模型估計(jì)精度。戶型是一個(gè)連續(xù)的數(shù)值變量,可以不用修改。其他的字符型變量則需要進(jìn)行量化。樓層包括3種,根據(jù)計(jì)量經(jīng)濟(jì)學(xué)知識,本方案將采用2個(gè)0-1變量(構(gòu)成一個(gè)二維行向量)度量該3種類型。例如,(0,0)代表低層,(1,0)代表中層,(0,1)代表高層。其他分類指標(biāo)同樣需要進(jìn)行類似處理。本文選取河北省邯鄲市某小區(qū)的樣本數(shù)據(jù)。

表2 河北省邯鄲市某小區(qū)的樣本數(shù)據(jù)Tab.2 Sample data of a residential district in Handan city, Hebei province

該實(shí)例的線性回歸模型為

經(jīng)過逐步回歸之后所得到的“最優(yōu)”回歸模型為

圖3 基于M估計(jì)的D統(tǒng)計(jì)量Fig.3 D StatisticsbBased on M estimation

1) 異常值檢測

從圖3可以看出,第24,38,42,43個(gè)樣本的對應(yīng)的D統(tǒng)計(jì)量大于3,可判斷為異常值,其中雖然第42和43的樣本連續(xù)異常,但仍然能夠被很好地檢測出來,再一次證明此方法對異常值的遮蔽現(xiàn)象有一定的作用,能夠有效地將連續(xù)幾個(gè)異常值檢測出來。

2) 異常值大小的估計(jì)

同樣,由公式(4)估計(jì)出異常值的大小,估計(jì)結(jié)果如表3所示。

表3 異常值大小的估計(jì)值Tab.3 Estimates of the size of outliers

6 結(jié) 語

異常值檢測是當(dāng)前數(shù)據(jù)分析研究中的一個(gè)熱點(diǎn)問題。通過大量的模擬實(shí)驗(yàn)和實(shí)例分析,得到這樣的結(jié)論:基于M估計(jì)的D統(tǒng)計(jì)量可以很好地檢測出異常值,尤其是對異常值的遮蔽現(xiàn)象有一定的作用----能夠有效地檢測出連續(xù)的幾個(gè)異常值。