聶興鋒,吳劉倉,邢伊琦

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

1.引言

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

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

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

Ⅰ偏正態(tài)分布

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

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

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

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

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

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

ⅢPena距離

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

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

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

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

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

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

模型(2.2)對應的Pena距離如下:

定理2.2當樣本中不含有異常點時,有

證

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

而當hjj ≥n1時,我們有

定理2.3當樣本中含有高杠桿異常點時,統(tǒng)計量Si的期望,有

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

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

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

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

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

為檢測第i個點是否為異常點或強影響點,可通過比較刪除第i個點前后統(tǒng)計推斷結(jié)果的變化,其中統(tǒng)計診斷量的變化可通過一些統(tǒng)計診斷量來刻畫.

Ⅱ極大似然估計

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

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

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

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

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

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

i) 似然距離及其計算

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

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

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

ii) Cook距離及其計算

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

iii) Pena距離及其計算

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

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 號,80 號,180號樣本點的被解釋變量的值做改變,即從樣本點中制造3個異常點,然后應用本文研究的方法如似然距離,Cook距離和Pena距離進行診斷.根據(jù)這3個異常點的診斷情況來判斷本文提出的方法是否行之有效.模擬結(jié)果如圖1,圖2和圖3所示,其中圖1為樣本量為200時模擬數(shù)據(jù)的似然距離LD散點圖,圖2樣本量為200時模擬數(shù)據(jù)的Cook距離散點圖,圖3樣本量為200時模擬數(shù)據(jù)的Pena距離散點圖.

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

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

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

從圖中我們可以很清晰的看出,第20,80,180號異常點均被診斷出來了,表明我們的方法是行之有效的,下面用實例進一步說明具體的應用.

5.實例分析

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

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

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

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

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

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

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

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

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

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

圖5 發(fā)動機性能數(shù)據(jù)似然距離LD散點圖

圖6 發(fā)動機性能數(shù)據(jù)Cook距離CD散點圖

圖7 發(fā)動機性能數(shù)據(jù)Pena距離PD散點圖

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

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

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

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