基于Pena距離的偏正態(tài)數(shù)據(jù)下位置回歸模型的統(tǒng)計(jì)診斷

聶興鋒,吳劉倉(cāng),邢伊琦

(昆明理工大學(xué)理學(xué)院,云南 昆明650093)

1.引言

在日常生活中,我們遇到的大多數(shù)數(shù)據(jù)并不具有嚴(yán)格的對(duì)稱性,而具有一定的偏斜,如果此時(shí)再用正態(tài)分布等對(duì)稱分布去描述它們的性質(zhì)就有點(diǎn)不恰當(dāng)了.目前,偏態(tài)數(shù)據(jù)的統(tǒng)計(jì)推斷成為統(tǒng)計(jì)學(xué)研究的一個(gè)熱點(diǎn)問(wèn)題之一.

我們知道,統(tǒng)計(jì)診斷是數(shù)據(jù)分析的第一步,主要目的就是對(duì)樣本數(shù)據(jù)中異常點(diǎn)或強(qiáng)影響點(diǎn)的識(shí)別和診斷.傳統(tǒng)的判斷異常點(diǎn)的常用統(tǒng)計(jì)量有Cook 距離、似然距離等.Pena[1]提出了一種度量線性回歸模型影響的新的方法,這種方法與之前的方法有較大區(qū)別,之前的方法是研究刪除一個(gè)(組)點(diǎn)對(duì)回歸分析的影響以及對(duì)模型預(yù)測(cè)值的影響,或者是某個(gè)(組)樣本點(diǎn)的微小擾動(dòng)對(duì)參數(shù)估計(jì)的影響或是對(duì)模型預(yù)測(cè)的影響；而Pena距離這一統(tǒng)計(jì)量是研究樣本中的某一點(diǎn)受其余各點(diǎn)的影響,也就是度量樣本中各點(diǎn)刪除后對(duì)某一特定點(diǎn)回歸值以及預(yù)測(cè)值的影響.孟麗麗等[2]研究了基于Pena距離的加權(quán)最小二乘估計(jì)的影響分析;胡江等[3]?[5]研究了基于Pena距離的非線性回歸模型和廣義線性回歸模型的影響分析.針對(duì)偏態(tài)數(shù)據(jù)的統(tǒng)計(jì)診斷方面,基于Cook距離、似然距離等,Xie等[6]研究了偏正態(tài)分布下非線性均值回歸模型的統(tǒng)計(jì)診斷;萬(wàn)文等[7]研究了偏正態(tài)數(shù)據(jù)下聯(lián)合位置與尺度模型的統(tǒng)計(jì)診斷.但是基于Pena距離的偏正態(tài)數(shù)據(jù)的統(tǒng)計(jì)診斷還沒(méi)有人研究,而統(tǒng)計(jì)診斷又是數(shù)據(jù)分析必不可少的一部分.本文對(duì)Pena距離在偏正態(tài)數(shù)據(jù)下位置回歸模型的影響分析進(jìn)行了討論,得出了比較有價(jià)值的相關(guān)結(jié)果.

2.偏正態(tài)分布下的位置回歸模型

Ⅰ偏正態(tài)分布

1985年Azzalini[8]首次研究提出偏正態(tài)分布,若隨機(jī)變量Y服從偏正態(tài)分布,即Y ～SN(μ,σ2,λ) 其中μ表示位置參數(shù),σ表示尺度參數(shù),λ表示偏度參數(shù).則其概率密度函數(shù)可表示為

其中φ(·),Φ(·)分別為標(biāo)準(zhǔn)正態(tài)分布的密度函數(shù)與分布函數(shù).當(dāng)偏度參數(shù)λ= 0 時(shí),密度函數(shù)(2.1)退化為正態(tài)分布的密度函數(shù),即此時(shí)偏正態(tài)分布退化為正態(tài)分布.

從E(Y)中我們可以看出μ只是均值的一部分.當(dāng)λ≠ 0時(shí),E(Y) =μ,此時(shí)分布不對(duì)稱;當(dāng)λ >0時(shí),E(Y)>μ,此時(shí)分布右偏;當(dāng)λ <0時(shí),E(Y)<μ,此時(shí)分布左偏.所以,偏正態(tài)分布是正態(tài)分布的進(jìn)一步推廣.

Ⅱ偏正態(tài)分布下的位置回歸模型

下面給出偏正態(tài)分布下的位置回歸模型為:

其中yi是被解釋變量,服從位置參數(shù)為μ,尺度參數(shù)為σ,偏度參數(shù)為λ的偏正態(tài)分布,xi=(xi1,xi2,...,xip)T是與yi有關(guān)的解釋變量.本文主要研究模型(2.2)的統(tǒng)計(jì)診斷方法.

ⅢPena距離

給定一組觀測(cè)數(shù)據(jù)(xi,yi)，i=1,...,n,yi為獨(dú)立服從SN分布的隨機(jī)變量,則位置回歸模型(2.2)可表示為:

其中xi=(xi1,xi2,...,xip)T.其向量形式為:

其中Xi= (1,xi1,xi2,xi3,...,xip)，X= (X1,X2,X3,...,Xn)為n×p的設(shè)計(jì)矩陣,維數(shù)為p,β為p×1的參數(shù)向量,ε為n×1的向量.則

其中(H=X(XTX)?1XT)是一個(gè)帽子矩陣,且有H2=H,HT=H.

定理2.1模型(2.2)的Pena距離為:

證根據(jù)文[1],我們定義Pena距離如下:其中由文[1]知:其中為第個(gè)i點(diǎn)的擬合值,是刪除第j個(gè)點(diǎn)后第i個(gè)點(diǎn)的擬合值為帽子矩陣H的對(duì)角元素(杠桿值),p為帽子矩陣H的維數(shù).所以有:

模型(2.2)對(duì)應(yīng)的Pena距離如下:

定理2.2當(dāng)樣本中不含有異常點(diǎn)時(shí),有

證

由韋博成等[9]可知: E(?r2j)=1,故

而當(dāng)hjj ≥n1時(shí),我們有

定理2.3當(dāng)樣本中含有高杠桿異常點(diǎn)時(shí),統(tǒng)計(jì)量Si的期望,有

1) E(Si)→0,高杠桿異常點(diǎn);

由定理2.3可知,當(dāng)數(shù)據(jù)中含有一簇相同的高杠異常點(diǎn)時(shí),可根據(jù)Si的值很容易找到它們但Cook 距離不能識(shí)別.特別,當(dāng)λ=0時(shí),g(0)=1,即可得到文[1-5]類似的結(jié)論.所以,本文進(jìn)一步拓展了文[1-5]在偏態(tài)數(shù)據(jù)的實(shí)際應(yīng)用.

3.偏正態(tài)數(shù)據(jù)下位置回歸模型的統(tǒng)計(jì)診斷

Ⅰ數(shù)據(jù)刪除模型

數(shù)據(jù)刪除是統(tǒng)計(jì)診斷中最常用的方法之一,比較第i個(gè)點(diǎn)刪除前后模型參數(shù)估計(jì)量之間的差異,能得出一些很好的結(jié)論.偏正態(tài)數(shù)據(jù)下位置回歸模型的刪除模型可表示為:

為檢測(cè)第i個(gè)點(diǎn)是否為異常點(diǎn)或強(qiáng)影響點(diǎn),可通過(guò)比較刪除第i個(gè)點(diǎn)前后統(tǒng)計(jì)推斷結(jié)果的變化,其中統(tǒng)計(jì)診斷量的變化可通過(guò)一些統(tǒng)計(jì)診斷量來(lái)刻畫.

Ⅱ極大似然估計(jì)

對(duì)于模型(2.2),假設(shè)(yi,xi)為數(shù)據(jù)集中的第i個(gè)數(shù)據(jù)點(diǎn),由模型(2.2)可知其密度函數(shù)為:

由(3.2)式可得似然函數(shù)為:

上式取自然對(duì)數(shù),得其對(duì)數(shù)似然函數(shù)為:

令θ=(βT,σ2,λ)T,則L(β,σ2,λ)=L(θ).因此

利用Gauss-Newton迭代法[10]可得到參數(shù)極大似然估計(jì)的估計(jì)值.設(shè)未刪除模型的參數(shù)估計(jì)值用表示刪除模型的參數(shù)估計(jì)值則可以用表示,即刪除第i個(gè)點(diǎn)后的參數(shù)估計(jì)值則表示刪除第j個(gè)數(shù)據(jù)點(diǎn)后第i個(gè)數(shù)據(jù)點(diǎn)的參數(shù)估計(jì)值.

Ⅲ基于數(shù)據(jù)刪除模型的診斷統(tǒng)計(jì)量

i) 似然距離及其計(jì)算

在數(shù)據(jù)刪除模型框架下,似然距離是與Cook距離同等重要的診斷統(tǒng)計(jì)量.由于似然距離的定義并不限于線性模型,故而可以用于相當(dāng)廣泛的統(tǒng)計(jì)模型,諸如非線性模型、廣義線性模型等.針對(duì)本文中的刪除模型(3.1),第i個(gè)點(diǎn)的似然距離定義為:

由于L()為全局最優(yōu)解,因此LDi ≥0恒成立.似然距離反應(yīng)了第i個(gè)數(shù)據(jù)點(diǎn)(xi,yi)對(duì)參數(shù)θ的極大似然估計(jì)的影響.對(duì)于遠(yuǎn)大于其似然距離的點(diǎn),則為異常點(diǎn)或強(qiáng)影響點(diǎn).由于似然距離沒(méi)有顯示解,因此需要用近似計(jì)算得出其數(shù)值解.對(duì)(3.5)式在處進(jìn)行泰勒展開可得:

其中I()為Fisher信息陣,為方便計(jì)算,本文使用Fisher觀測(cè)陣計(jì)算似然距離LD?i.

ii) Cook距離及其計(jì)算

Cook距離是當(dāng)今統(tǒng)計(jì)診斷中最重要的診斷統(tǒng)計(jì)量之一.針對(duì)本文中的刪除模型(3.1),第i個(gè)點(diǎn)的Cook距離定義如下:其中H=X(XTX)?1XT為帽子矩陣,p為對(duì)應(yīng)解釋變量的維數(shù),?σ2為未刪除模型方差的估計(jì)值.具體分析時(shí),先計(jì)算出各點(diǎn)的Cook距離,通過(guò)畫散點(diǎn)圖,找出其中特別大的,對(duì)應(yīng)的數(shù)據(jù)點(diǎn)可能就是異常點(diǎn)或強(qiáng)影響點(diǎn).

iii) Pena距離及其計(jì)算

Cook距離研究的是刪除一個(gè)(組)點(diǎn)后對(duì)估計(jì)值或預(yù)測(cè)值的影響,而Pena距離則研究的是樣本中的某一點(diǎn)受其余各點(diǎn)的影響,簡(jiǎn)單來(lái)說(shuō),就是樣本中各點(diǎn)刪除后,對(duì)某一特定的點(diǎn)的估計(jì)值或預(yù)測(cè)值的影響,Pena距離定義如下：其中H=X(XTX)?1XT為帽子矩陣,p為對(duì)應(yīng)解釋變量的維數(shù),為刪除第i個(gè)點(diǎn)后模型方差的估計(jì)值.則表示刪除第j個(gè)數(shù)據(jù)點(diǎn)后第i個(gè)數(shù)據(jù)點(diǎn)的參數(shù)估計(jì)值.具體分析時(shí),同樣是先算出刪除各點(diǎn)后某一點(diǎn)的Si,畫出散點(diǎn)圖,Si較大的則可能是異常點(diǎn)或強(qiáng)影響點(diǎn).

4.Monte Carlo模擬

為了比較Pena距離與Cook距離、似然距離的診斷效果,根據(jù)模型(2.2),產(chǎn)生偏正態(tài)數(shù)據(jù),其中xi ～U(?1,1),取β=(1,1,1),σ=2,λ=0.5.將第20 號(hào),80 號(hào),180號(hào)樣本點(diǎn)的被解釋變量的值做改變,即從樣本點(diǎn)中制造3個(gè)異常點(diǎn),然后應(yīng)用本文研究的方法如似然距離,Cook距離和Pena距離進(jìn)行診斷.根據(jù)這3個(gè)異常點(diǎn)的診斷情況來(lái)判斷本文提出的方法是否行之有效.模擬結(jié)果如圖1,圖2和圖3所示,其中圖1為樣本量為200時(shí)模擬數(shù)據(jù)的似然距離LD散點(diǎn)圖,圖2樣本量為200時(shí)模擬數(shù)據(jù)的Cook距離散點(diǎn)圖,圖3樣本量為200時(shí)模擬數(shù)據(jù)的Pena距離散點(diǎn)圖.

圖1 樣本量為200時(shí)模擬數(shù)據(jù)的LD散點(diǎn)圖

圖2 樣本量為200時(shí)模擬數(shù)據(jù)的CD散點(diǎn)圖

圖3 樣本量為200時(shí)模擬數(shù)據(jù)的PD散點(diǎn)圖

從圖中我們可以很清晰的看出,第20,80,180號(hào)異常點(diǎn)均被診斷出來(lái)了,表明我們的方法是行之有效的,下面用實(shí)例進(jìn)一步說(shuō)明具體的應(yīng)用.

5.實(shí)例分析

Ⅰ發(fā)動(dòng)機(jī)性能數(shù)據(jù)[11]

如下表1所示是一組檢驗(yàn)?zāi)撤N工業(yè)用發(fā)電機(jī)性能試驗(yàn)的數(shù)據(jù),該試驗(yàn)使用的原料是柴油和從有機(jī)原料中通過(guò)蒸餾產(chǎn)生的氣體的混合物,在各種不同的速度x(計(jì)量單位：百轉(zhuǎn)/分鐘)下,測(cè)量發(fā)動(dòng)機(jī)的馬力y.

用QQ圖對(duì)發(fā)動(dòng)機(jī)的馬力y數(shù)據(jù)進(jìn)行正態(tài)性檢驗(yàn),結(jié)果如圖4所示,表明數(shù)據(jù)具有一定的偏斜.利用MATLAB中的偏度函數(shù)skewness(),峰度函數(shù)kurtosis()得到發(fā)動(dòng)機(jī)的馬力y的偏度為-0.3332,峰度為1.9679,而正態(tài)分布的偏度值為0,峰度值為3.綜合分析可知,發(fā)動(dòng)機(jī)性能數(shù)據(jù)服從偏態(tài)分布,可用本文研究的方法進(jìn)行統(tǒng)計(jì)診斷.

表1 發(fā)動(dòng)機(jī)性能數(shù)據(jù)

本文考慮發(fā)動(dòng)機(jī)的馬力y與在各種不同的速度x(計(jì)量單位：百轉(zhuǎn)/分鐘)的位置回歸模型.經(jīng)過(guò)計(jì)算得到完全數(shù)據(jù)下模型(2.2)的參數(shù)估計(jì)結(jié)果如下:

由圖5可知第2,10,17,24號(hào)點(diǎn)可能為異常點(diǎn)或強(qiáng)影響點(diǎn),由圖6可知第2,10,24號(hào)點(diǎn)可能為異常點(diǎn)或強(qiáng)影響點(diǎn),由圖7可知第2,24號(hào)點(diǎn)可能為強(qiáng)影響點(diǎn)或異常點(diǎn).由韋博成等[9]的例5.4可知第2,24號(hào)點(diǎn)為異常點(diǎn)或強(qiáng)影響點(diǎn).比起似然距離和Cook距離,Pena距離很好的診斷出了這兩個(gè)點(diǎn).

Ⅱ紅鱒鮭魚數(shù)據(jù)[12]

魚卵數(shù)量x當(dāng)年可捕撈的成魚數(shù)量y之間的關(guān)系,是經(jīng)營(yíng)漁場(chǎng)者十分關(guān)心的問(wèn)題.下表2所示是1940年至1967年在Skeener河中紅鱒鮭魚的產(chǎn)卵量x和可捕撈的成魚量y的測(cè)量數(shù)據(jù).

表2 紅鱒鮭魚數(shù)據(jù)

圖4 發(fā)動(dòng)機(jī)性能數(shù)據(jù)的正態(tài)性檢驗(yàn)QQ圖

圖5 發(fā)動(dòng)機(jī)性能數(shù)據(jù)似然距離LD散點(diǎn)圖

圖6 發(fā)動(dòng)機(jī)性能數(shù)據(jù)Cook距離CD散點(diǎn)圖

圖7 發(fā)動(dòng)機(jī)性能數(shù)據(jù)Pena距離PD散點(diǎn)圖

利用MATLAB中的偏度函數(shù)skewness()、峰度函數(shù)kurtosis()得到紅鱒鮭魚當(dāng)年可捕撈的成魚數(shù)量y的偏度為0.7063,峰度為3.0568,而正態(tài)分布的偏度值為0,峰度值為3.綜合分析可知,紅鱒鮭魚當(dāng)年可捕撈的成魚數(shù)量y服從偏態(tài)分布.我們分別用正態(tài)分布下的Pena距離和偏正態(tài)分布下的Pena距離診斷做比較,比較結(jié)果如圖8,圖9所示.

圖8 正態(tài)分布下的Pena距離散點(diǎn)圖

圖9 偏正態(tài)分布下的Pena距離散點(diǎn)圖

從圖8我們可以看出第5號(hào)點(diǎn)為異常點(diǎn)或強(qiáng)影響點(diǎn),而從圖9可以看出第5,12號(hào)點(diǎn)為異常點(diǎn)或強(qiáng)影響點(diǎn).由文[9]中例6.4可知第5,12號(hào)點(diǎn)為異常點(diǎn)或強(qiáng)影響點(diǎn),這是合理的,因?yàn)樵谠紨?shù)據(jù)中,第5,12 號(hào)點(diǎn)分別是被解釋變量的最大值點(diǎn)和最小值點(diǎn).偏正態(tài)分布下的Pena距離很好的診斷出了這兩個(gè)點(diǎn),而正態(tài)分布下的Pena 距離則只診斷出了一個(gè)點(diǎn).相比較而言,偏正態(tài)分布下的Pena距離診斷效果比正態(tài)分布下的Pena 距離要好.