基于Pena距離的偏正態數據下位置回歸模型的統計診斷

聶興鋒,吳劉倉,邢伊琦

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

1.引言

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

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

2.偏正態分布下的位置回歸模型

Ⅰ偏正態分布

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

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

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

Ⅱ偏正態分布下的位置回歸模型

下面給出偏正態分布下的位置回歸模型為:

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

ⅢPena距離

給定一組觀測數據(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的設計矩陣,維數為p,β為p×1的參數向量,ε為n×1的向量.則

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

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

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

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

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

證

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

而當hjj ≥n1時,我們有

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

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

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

3.偏正態數據下位置回歸模型的統計診斷

Ⅰ數據刪除模型

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

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

Ⅱ極大似然估計

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

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

上式取自然對數,得其對數似然函數為:

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

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

Ⅲ基于數據刪除模型的診斷統計量

i) 似然距離及其計算

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

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

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

ii) Cook距離及其計算

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

iii) Pena距離及其計算

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

4.Monte Carlo模擬

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

圖1 樣本量為200時模擬數據的LD散點圖

圖2 樣本量為200時模擬數據的CD散點圖

圖3 樣本量為200時模擬數據的PD散點圖

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

5.實例分析

Ⅰ發動機性能數據[11]

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

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

表1 發動機性能數據

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

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

Ⅱ紅鱒鮭魚數據[12]

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

表2 紅鱒鮭魚數據

圖4 發動機性能數據的正態性檢驗QQ圖

圖5 發動機性能數據似然距離LD散點圖

圖6 發動機性能數據Cook距離CD散點圖

圖7 發動機性能數據Pena距離PD散點圖

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

圖8 正態分布下的Pena距離散點圖

圖9 偏正態分布下的Pena距離散點圖

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