偏正態(tài)數(shù)據(jù)下聯(lián)合位置、尺度、偏度模型的統(tǒng)計診斷

吳劉倉, 聶興鋒, 鄭桂芬

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

1 引言

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

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

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

2 偏正態(tài)數(shù)據(jù)下聯(lián)合位置、尺度、偏度模型的極大似然估計

2.1 偏正態(tài)分布

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

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

2.2 基于偏正態(tài)數(shù)據(jù)的聯(lián)合位置、尺度、偏度模型

本文研究以下偏正態(tài)數(shù)據(jù)的聯(lián)合位置、尺度、偏度模型

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

2.3 極大似然估計

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

其中

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

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

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

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

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

2.4.1 似然距離及其計算

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

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

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

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

2.4.2 Cook 距離及其計算

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

2.4.3 Pena 距離及其計算

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

根據(jù)文獻[9]，Pena 距離定義如下

定理1偏正態(tài)數(shù)據(jù)下的Pena 距離

其中為第j個學(xué)生化殘差(標(biāo)準(zhǔn)化殘差).

證明 根據(jù)文獻[9]

定理2當(dāng)樣本中不含有異常點時，有

證明

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

記

我們有

綜上所述，當(dāng)樣本中含有高杠桿異常點時，統(tǒng)計量Si的期望，有：

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

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

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

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

3 局部影響分析

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

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

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

上式的二階近似表達式為

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

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

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

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

3.1 位置漂移擾動模型

位置漂移擾動模型如下

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

其中

L(θ|ω)的前二階導(dǎo)數(shù)可表示為

3.2 尺度加權(quán)擾動模型

尺度加權(quán)擾動模型如下所示

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

其中

L(θ|ω)的前二階導(dǎo)數(shù)可表示為

3.3 偏度加權(quán)擾動模型

偏度加權(quán)擾動模型如下

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

其中

L(θ|ω)的前二階導(dǎo)數(shù)可表示為

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

4 模擬研究及實例分析

4.1 模擬研究

為了檢驗本文提出方法的有效性，根據(jù)模型(2)我們產(chǎn)生隨機數(shù)據(jù)，其中yi ～SN(μi,,λi)，xi, zi, hi均產(chǎn)生于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 號點設(shè)為異常點.然后根據(jù)上述的診斷方法得出模擬結(jié)果，模擬結(jié)果如圖1 至圖3 所示.

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

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

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

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

4.2 虹鱒鮭魚數(shù)據(jù)[25]分析

魚卵數(shù)量x和當(dāng)年可捕撈成魚數(shù)量y之間的關(guān)系，是養(yǎng)殖者非常關(guān)心的問題.下表1 所示是1940 年至1967 年在Skeener 河中紅鱒鮭魚的產(chǎn)卵量x和可捕撈的成魚量y的測量數(shù)據(jù).

表1 虹鱒鮭魚數(shù)據(jù)

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

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

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

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

圖4 虹鱒鮭魚數(shù)據(jù)的正態(tài)性檢驗QQ 圖

圖5 虹鱒鮭魚數(shù)據(jù)的LD 散點圖

圖6 虹鱒鮭魚數(shù)據(jù)的CD 散點圖

圖7 虹鱒鮭魚數(shù)據(jù)的PD 散點圖

4.3 虹鱒鮭魚數(shù)據(jù)的局部影響分析

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

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

5 結(jié)論

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

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

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

圖10 數(shù)據(jù)尺度加權(quán)擾動下?i(i)散點圖

圖11 數(shù)據(jù)尺度加權(quán)擾動下|(dmax)i|散點圖

圖12 數(shù)據(jù)偏度加權(quán)擾動下?(i)散點圖

圖13 數(shù)據(jù)偏度加權(quán)擾動下|(dmax)i|散點圖