偏正態數據下聯合位置、尺度、偏度模型的統計診斷

吳劉倉, 聶興鋒, 鄭桂芬

(昆明理工大學理學院，昆明 650093)

1 引言

目前，對均值參數建模的理論和方法都研究得比較透徹[1-6]，但在現實生活中我們會發現，同方差數據只占少數部分，大多數還是異方差數據，這表明對方差參數建模同樣很重要，對方差參數建模能夠很好的了解方差的來源，以此來達到有效的控制方差[7].另一方面，在這大多數的異方差數據中大部分并不具有嚴格的對稱性，而是具有一定偏斜的，這時我們再用正態分布去描述它們的性質就不太適合了.近年來，偏正態分布[8]成為非對稱分布研究的重要分支.因此，偏態數據的統計推斷也隨之成為統計學探索的熱點問題.

我們知道，統計診斷是數據分析的第一步，主要目的就是對樣本數據中異常點或強影響點的診斷與識別.傳統的判斷異常點的方法有Cook 距離、似然距離、W-K 統計量、AP 統計量等.而文獻[9]中提出了一種新方法Pena 距離，文獻[9]所提方法與傳統方法有較大差別.傳統的方法是刪除一個樣本點，對估計值的影響，或者是某個樣本點的微小擾動對估計值的影響，而Pena 距離則是研究樣本中某一個樣本點受其余各個樣本點的影響，簡單來說，就是樣本中各點刪除后，對某一特定的點的預測值的影響.

Pena 距離的研究方面：孟麗麗和盧志義[10]基于Pena 距離關于加權最小二乘估計的影響分析做了研究；胡江等[11-13]基于Pena 距離研究了非線性回歸模型以及廣義回歸模型的影響分析.異方差的研究方面：Aitkin[14]基于異方差模型研究了正態分布下的極大似然估計；戴琳等[15]基于聯合均值與方差模型研究了統計診斷；馬婷等[16]基于SN 分布聯合位置、尺度、偏度模型研究了極大似然估計；吳劉倉等[17,18]基于StN 分布下聯合位置、尺度、偏度模型研究了極大似然估計以及基于偏正態數據下聯合位置與尺度混合專家回歸模型研究了參數估計；Lachos 等[19]基于SN 分布的混合尺度異方差非線性回歸模型研究了參數估計與統計診斷.偏態數據的統計診斷方面：基于Cook 距離、似然距離等，Xie 等[20]研究了SN 分布下非線性均值回歸模型的統計診斷；萬文等[21]基于偏正態數據下聯合位置與尺度模型研究了統計診斷.李玲雪等[22]缺失偏態數據下聯合位置與尺度模型研究了統計推斷；李世凱等[23]偏態數據下混合非線性回歸模型研究了統計推斷.但是偏正態數據下聯合位置、尺度、偏度模型的異常點診斷和識別還沒有人研究，而統計診斷又是處理數據必不可少的一部分.因此，基于SN 分布下聯合位置、尺度、偏度模型的統計診斷進行研究，得出了比較有價值的相關結果.

2 偏正態數據下聯合位置、尺度、偏度模型的極大似然估計

2.1 偏正態分布

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

其中φ(·),Φ(·)分別表示標準的正態分布的密度函數和分布函數.當偏度參數λ= 0 時，密度函數(1)退化為正態分布的密度函數，即此時偏正態分布退化為正態分布.在偏正態分布中有

2.2 基于偏正態數據的聯合位置、尺度、偏度模型

本文研究以下偏正態數據的聯合位置、尺度、偏度模型

其中yi是被解釋變量，服從SN 分布，xi=(xi1,xi2,··· ,xip)T, zi=(zi1,zi2,··· ,ziq)T,hi= (hi1,hi2,··· ,hir)T是與yi有關的解釋變量.β= (β1,β2,··· ,βp)T是與位置模型有關的p× 1 維向量，γ= (γ1,γ2,··· ,γq)T是與尺度模型有關的q× 1 維向量，α=(α1,α2,··· ,αr)T是與偏度模型有關的r×1 維向量，xi, zi, hi三個解釋變量的觀測值可以相同，可以不同或部分相同.本篇文章主要研究模型(2)的統計診斷方法.

2.3 極大似然估計

假設(yi,xi,zi,hi)為樣本中的第i個樣本點，由密度函數(1)及模型(2)可知其密度函數為

其中

由(3)式可得其似然函數表示為

上式兩邊取自然對數，就得到對數似然函數為

令θ=(βT,γT,αT)T，則L(β,γ,α)=L(θ)，因此

由Gauss-Newton 迭代方法可獲得相應參數的估計值.設我們未刪除模型參數估計值用,,表示，.刪除模型的參數估計值則可以用表示，，即刪除第i個點后的參數估計值.，則表示刪除第j個數據點后第i個數據點的參數估計值.

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

2.4.1 似然距離及其計算

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

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

由于˙L(?θ)=0，從而得到似然距離(LD)的近似表達式如下

其中I()為Fisher 信息陣，為方便計算，本文使用Fisher 陣計算似然距離L.

2.4.2 Cook 距離及其計算

Cook 距離是統計診斷中非常重要的診斷統計量之一，是Cook 于1977 年基于參數置信域的統計意義提出來的.針對本文中的刪除模型，第i個點的Cook 距離定義如下

2.4.3 Pena 距離及其計算

Pena 距離是Pena 于2005 年提出的一種診斷統計量.常見的統計診斷方法有數據刪除診斷和局部影響分析，數據刪除是針對完全數據的，主要考察刪除某一個或某一組樣本點對既定模型參數估計的影響，即對模型回歸分析和預測分析的影響.局部影響分析是對模型施加一個微小擾動，然后研究樣本點對模型參數估計和預測的影響.而Pena 距離則是研究樣本中某一點受其余各點的影響，即刪除樣本中的各點對某一既定樣本點的預測值的影響，也是基于數據刪除模型的診斷分析，是對診斷統計量的重要補充.

根據文獻[9]，Pena 距離定義如下

定理1偏正態數據下的Pena 距離

其中為第j個學生化殘差(標準化殘差).

證明 根據文獻[9]

定理2當樣本中不含有異常點時，有

證明

由文獻[24]可知：E()=1，故

記

我們有

綜上所述，當樣本中含有高杠桿異常點時，統計量Si的期望，有：

1)E(Si)→0，高杠桿異常點；

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

Pena 距離其向量形式定義如下

其中H=X(XTX)?1XT為帽子矩陣，p為對應解釋變量的維數，exp(zTi?γi)為刪除第i個點后模型方差的估計值.?θi(j)則表示刪除第j個數據點后第i個數據點的參數估計值.具體分析時，同樣是先算出刪除各點后某一點的Si，畫出散點圖，Si較大的則可能是異常點.

3 局部影響分析

局部影響分析是1986 年Cook 首次提出來的一種很有用的統計診斷方法，其主要思想是引入一個微小擾動(干擾)的概念，而把異常點或強影響點歸結為“比其他點受到更大干擾的點”.若數據集中的一個或幾個點比其他數據點受到的擾動大，則這個或這幾個數據點就是異常點.

由文獻[24]可知，將未受到干擾的模型(2)記為D，其擾動模型記為D(ω), ω=(ω1,ω2,··· ,ωn)T為刻畫各樣本點擾動大小的向量.對于受到擾動的模型，記其對數似然函數為L(θ/ω)，參數的極大似然估計為?θω.另外，存在ω0=ω，使得D(ω0)=D，表示模型未受到擾動.

對于擾動模型D(ω)其似然距離定義為

上式的二階近似表達式為

其中?¨F稱為影響矩陣，表達式如下

其中?是(p+q+r)×n階矩陣，是n×n階矩陣.

LD(ω)?反映的是第i個樣本點的擾動對極大似然估計的影響，其數值越大，就表示這個樣本點對估計值的影響越大，如果存在某點j的擾動特別大，則這個點就是異常點.

最大特征向量法：LD(ω)?關于方向d=ω?ω0的最大值，并設d=dmax時，LD(ω)?達到最大值.dmax= (d1,d1,··· ,dn)，并假設其中有一分量|dj|的值比其他分量大得多，則說明dj對于是的似然距離達到最大值做出了最大貢獻，因而對應的數據點(yi,xi,zi,hi)即為異常點.因此，可作(i,|(dmax)i|)的散點圖來找出異常點.

3.1 位置漂移擾動模型

位置漂移擾動模型如下

其中ω=(ω1,ω2,··· ,ωn)T, ω0=(0,0,··· ,0)T，表示模型沒有擾動.θ=(βT,γT,αT)T，對應的對數似然函數表達示為

其中

L(θ|ω)的前二階導數可表示為

3.2 尺度加權擾動模型

尺度加權擾動模型如下所示

其中ω=(ω1,ω2,··· ,ωn)T, ω0=(0,0,··· ,0)T，表示模型沒有擾動.θ=(βT,γT,αT)T，對應的對數似然函數可表達示為

其中

L(θ|ω)的前二階導數可表示為

3.3 偏度加權擾動模型

偏度加權擾動模型如下

其中ω=(ω1,ω2,··· ,ωn)T, ω0=(0,0,··· ,0)T，表示模型沒有擾動.θ=(βT,γT,αT)T，對應的對數似然函數可表達示為

其中

L(θ|ω)的前二階導數可表示為

其中?(θ)是(p+q+r)×n階矩陣，是n×n階矩陣.

4 模擬研究及實例分析

4.1 模擬研究

為了檢驗本文提出方法的有效性，根據模型(2)我們產生隨機數據，其中yi ～SN(μi,,λi)，xi, zi, hi均產生于U(?1,1), β, γ, α的真值分別取為(0.6,0.5,0.8)T,(1,2,?1)T,(2,0.8,?0.5)T，樣本量n為200，并把38、132、180 號點設為異常點.然后根據上述的診斷方法得出模擬結果，模擬結果如圖1 至圖3 所示.

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

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

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

從圖中我們可以很清晰地看出38、132、180 號異常點均被診斷出來了，說明我們的方法是行之有效的.下面用實例進一步說明.

4.2 虹鱒鮭魚數據[25]分析

魚卵數量x和當年可捕撈成魚數量y之間的關系，是養殖者非常關心的問題.下表1 所示是1940 年至1967 年在Skeener 河中紅鱒鮭魚的產卵量x和可捕撈的成魚量y的測量數據.

表1 虹鱒鮭魚數據

用QQ 圖對虹鱒鮭魚數據進行正態性檢驗，得到如下圖4，利用Matlab 中的偏度函數skewness( )、峰度函數kurtosis( )得到實例數據的偏度為0.7063、峰度為3.0568，而正態分布的偏度值為0，峰度值為3.綜合分析可知，虹鱒鮭魚數據近似服從偏正態分布，可用本文研究的方法進行統計診斷.

考慮魚卵數量x與可捕撈成魚量y之間的聯合位置、尺度、偏度模型，在該模型中xi, zi, hi完全相同，通過Matlab 計算得到完全數據下模型(2)的參數估計結果如下

為了判斷該數據集中哪些點是異常點，我們通過似然距離、Cook 距離、Pena 距離三種診斷統計量來診斷，診斷結果如圖5 至圖7 所示.

從圖5 我們可以看出5、12、18、22、25 號點可能為異常點，從圖6 可以看出12、18號點可能為異常點，而從圖7 可以看出5、12、16、18 號點可能為異常點.由文獻[24]統計診斷例6.4 可知5、12 號點為異常點，這是合理的，因為在原始數據中，第5、12 號點分別是被解釋變量的最大值點和最小值點.而16、18 號點從表1 中我們也可以看出魚卵量和可捕撈成魚量與其他點明顯異常.相較而言，Pena 距離診斷效果要比似然距離和Cook 距離要更精確一點.

圖4 虹鱒鮭魚數據的正態性檢驗QQ 圖

圖5 虹鱒鮭魚數據的LD 散點圖

圖6 虹鱒鮭魚數據的CD 散點圖

圖7 虹鱒鮭魚數據的PD 散點圖

4.3 虹鱒鮭魚數據的局部影響分析

采用位置漂移擾動、尺度加權擾動、偏度加權擾動三種方法對虹鱒鮭魚數據進行診斷，通過計算得到如下圖8 至圖13 所示結果.

從圖8 和圖9 可以看出12 號點為強影響點或異常點；從圖10 可以看出5、12、18 號點為異常點或強影響點，從圖11 可以看出4、5、12、18 號點為強影響點或異常點；從圖12、圖13 可以看出12 號點為強影響點或異常點.而根據4.2 實例可知5、12、18 號點為異常點或強影響點，所以我們的局部影響分析中，尺度加權擾動模型的診斷效果比較好，位置漂移擾動模型和偏度加權擾動模型的效果略差一點.

5 結論

本文研究了偏正態數據的聯合位置、尺度、偏度模型，將Pena 距離從正態推廣到了偏正態，適用范圍更廣.利用Pena 距離、Cook 距離、似然距離以及局部影響分析進行診斷，得到了在一定條件下Pena 距離優于Cook 距離和似然距離.并做了局部影響分析，結果表明尺度加權模型的診斷效果較好，位置漂移擾動模型和偏度加權模型的效果略差一點.

圖8 數據位置漂移擾動下?i(i)散點圖

圖9 數據位置漂移擾動下|(dmax)i|散點圖

圖10 數據尺度加權擾動下?i(i)散點圖

圖11 數據尺度加權擾動下|(dmax)i|散點圖

圖12 數據偏度加權擾動下?(i)散點圖

圖13 數據偏度加權擾動下|(dmax)i|散點圖