,,

(上海理工大學 理學院,上海 200093)

隨著數(shù)據(jù)獲取技術(shù)的迅猛發(fā)展,人們獲取的數(shù)據(jù)結(jié)構(gòu)越來越復雜,維數(shù)越來越高.統(tǒng)計學的主要任務就是對觀測數(shù)據(jù)的因變量和自變量建立模型,進而對數(shù)據(jù)進行分析、預測以及一些統(tǒng)計推斷.在現(xiàn)實問題中,因變量往往同時受多個自變量影響,但這些影響并不都很顯著.人們通常希望在模型中只引進對因變量有重要影響的自變量,所以,變量選擇就成了建模前的必要工作.但是,現(xiàn)實問題中,數(shù)據(jù)經(jīng)常被污染,往往存在異常值,這時用普通的變量選擇方法就會對模擬結(jié)果產(chǎn)生很大的偏差.

針對變量選擇的問題,統(tǒng)計學家已經(jīng)作出了大量研究.1996年,統(tǒng)計學家Tibshirani[1]提出了一種變量選擇方法LASSO,基本思想是在最小二乘法的基礎上施加L1懲罰.2001年,Fan等[2]提出了變量選擇的SCAD方法,并研究了該方法的Oracle性質(zhì).在某些LASSO不相合的情況下,Zou[3]又提出Adaptive LASSO,該方法是對LASSO的一種改進,能夠滿足Oracle性質(zhì).為了克服LASSO的一些缺點,Zou等[4]提出了Elastic Net變量選擇方法.針對高維數(shù)據(jù),Candès等[5]提出了Dantzig Selector方法.

針對數(shù)據(jù)中可能存在異常值這一情況,有許多文獻已經(jīng)研究了穩(wěn)健估計與穩(wěn)健變量選擇方法.文獻[6-7]率先提出當正態(tài)分布被污染時,估計位置參數(shù)的漸進理論.文獻[8-9]將最小一乘法用到穩(wěn)健估計中,之后文獻[10]進一步分析最小一乘法的優(yōu)良性質(zhì).文獻[11]提出了基于t函數(shù)的穩(wěn)健估計方法,考察了基于t函數(shù)估計量的優(yōu)良性.同時研究穩(wěn)健估計方法的還有文獻[12-14].文獻[15]提出了基于Huber函數(shù)的針對縱向數(shù)據(jù)的穩(wěn)健變量選擇方法.針對穩(wěn)健估計中常用的t函數(shù)和Huber函數(shù),文獻[16]提出了基于M估計的穩(wěn)健向前變量選擇方法,并進一步考察了t函數(shù)和Huber函數(shù)在穩(wěn)健向前變量選擇中的性質(zhì).

本文在前人研究的基礎上,提出一種新的基于t函數(shù)的穩(wěn)健變量選擇方法,并與文獻[15]中基于Huber函數(shù)的穩(wěn)健變量選擇方法進行比較.模擬結(jié)果顯示,t函數(shù)方法對數(shù)據(jù)中的異常值有更好的限制作用,可以達到更好的變量選擇效果.文章主要分為5個部分,第1部分介紹了穩(wěn)健的懲罰估計方程.第2部分將t函數(shù)和Huber函數(shù)的性質(zhì)進行比較分析,突出t函數(shù)在穩(wěn)健變量選擇方法中的優(yōu)勢.第3部分介紹本文中所用的算法.第4部分是數(shù)值模擬,通過3種污染方式來污染數(shù)據(jù),比較本文方法與文獻[15]中方法的模擬效果.第5部分為結(jié)論.

1 穩(wěn)健的懲罰估計方程

考慮如下線性模型:

y=xβ+ε

(1)

式中:y=(y1,y2,…,yn)T;x=(x1,x2…,xn)T,xi=(xi1,xi2,…,xip);β=(β1,β2,…,βp)T;ε=(ε1,ε2,…,εn)T,且i=1,2,…,n,εi的期望值為0,方差為σ2,ε的各分量相互獨立.

與文獻[12]類似,考慮如下穩(wěn)健估計方程

(2)

式中,Wi是權(quán)重矩陣W的第i個分量,權(quán)重Wi通過文獻[13]得來,用來降低自變量中異常值的影響,定義如下:

(3)

式中:r為大于1的常數(shù);p0為自由度與xi維數(shù)相同的卡方分布的0.95分位數(shù);取mx為xi的中位數(shù),則mx的第k個分量取為x第k列的中位數(shù);Sx的第k個對角元取為1.483(median|x(k)-mx(k)?In|),x(k)表示x的第k列,mx(k)表示mx的第k個分量,?表示kronecker乘積,In表示n維元素全為1的列向量.

式(2)中,函數(shù)φ(·)是一個有界得分函數(shù),用來限制因變量中異常值的影響,本文將此函數(shù)定義為自由度為2的t函數(shù),記作t2函數(shù).當φ(x)=x且Wi=1時,原穩(wěn)健估計方程就退化成一般的估計方程,不再具有穩(wěn)健性,即為非穩(wěn)健的估計方程,此時估計方程(2)會對異常值有較大的敏感性.

通過求解式(2),可以得到穩(wěn)健參數(shù)估計.為了同時達到變量選擇的效果,采用壓縮估計方法,即在估計方程中再添加一個懲罰項.因此考慮懲罰穩(wěn)健估計方程

nqλ(|β|)sgn(β)=0

(4)

本文所考慮的懲罰函數(shù)主要是SCAD懲罰函數(shù)[3].取懲罰函數(shù)為

2 t函數(shù)與Huber函數(shù)

t分布的密度函數(shù)為

Huber分布的密度函數(shù)為

其中

自由度不同,t分布密度函數(shù)的尾部厚度不同,從而t函數(shù)對異常值的抑制效果不同.圖1是自由度分別為2,6,10的t分布密度函數(shù),由圖1可見,自由度越小,密度函數(shù)的尾部越厚.文獻[16]已經(jīng)證明厚尾性對異常值有更好的抑制作用.圖2是自由度為2的t分布密度函數(shù)和Huber分布密度函數(shù)的圖像比較,由圖2顯然可見,t分布密度函數(shù)的尾部更厚.由此可以初步推斷,基于自由度為2的t函數(shù)的穩(wěn)健壓縮估計可以對異常值有更好的限制作用.下面,進一步分析t函數(shù)和Huber函數(shù)的圖像區(qū)別,以及通過圖像顯現(xiàn)出來的對異常值的作用效果,如圖3和圖4所示.

圖1 不同自由度的t分布密度函數(shù)Fig.1 t distribution density function with different kinds of degree of freedom

圖2 t分布和Huber分布密度函數(shù)比較Fig.2 Comparison density function between t distribution and Huber distribution

圖3 t函數(shù)Fig.3 t function

圖4 Huber函數(shù)Fig.4 Huber function

圖3是自由度分別為2,6,10的t函數(shù)的圖像,由圖3可見,t函數(shù)并不是單調(diào)的,而是一個回降函數(shù),而且隨著自由度的增大,在一定自變量范圍內(nèi),函數(shù)的變化范圍在變小.可以看出,當變量t的絕對值變大時,t函數(shù)將會對這些絕對值較大的變量產(chǎn)生作用,使其函數(shù)值接近于0,因此t函數(shù)可以很好地抑制數(shù)值模擬中異常值的影響.盡管自由度不同的t函數(shù)對較大異常值的抑制作用不相上下,但圖3表明,當變量t處在正常值范圍內(nèi)時,自由度越大的t函數(shù)會對變量施加越大的抑制作用,使原本正常的數(shù)據(jù)也受到更大的影響,從而破壞了數(shù)據(jù)原有的真實性.所以,綜合而言,t2函數(shù)是最優(yōu)的.本文將選取自由度v=2.

文獻[11]也模擬分析了自由度分別為1和4的t函數(shù)在M估計中的穩(wěn)健性,其模擬結(jié)果表明,自由度為1的t函數(shù)比自由度為4的t函數(shù)具有更好的穩(wěn)健性.本文也模擬分析了自由度分別為1和2的t函數(shù)在變量選擇中的穩(wěn)健性.結(jié)果表明,自由度為1的t函數(shù)和自由度為2的t函數(shù)的模擬結(jié)果比較接近.在模擬設置下,自由度為2的t函數(shù)比自由度為1的t函數(shù)在變量選擇和參數(shù)估計方面稍好一些,因此本文只報告自由度為2的t函數(shù)的結(jié)果.文獻[16]也得出了與本文相類似的結(jié)論,它們的研究表明,在M估計中,取自由度較小的t函數(shù)對異常值有較好的限制作用.

通過圖3和圖4的比較可見,當自變量趨于正無窮時,Huber函數(shù)值為+2,當自變量趨于負無窮時,Huber函數(shù)值為-2.而無論自變量趨于正無窮還是負無窮,t2函數(shù)值始終趨近于0.因此,t2函數(shù)的穩(wěn)健方法能減小異常值在模型估計中的作用,更好地削弱異常值的影響[14].所以,t2函數(shù)在變量選擇中比Huber函數(shù)具有更好的穩(wěn)健性.

3 算法

本文算法與文獻[15]類似,采用牛頓迭代法,具體算法如下:

a.對給定的一個λ值,首先計算β的初始值β(0).本文取最小二乘估計作為初始值.

b.用β(k)表示第k次迭代所得的估計值,k≥0,則

β(k+1)=β(k)-[Dβ(k)-Δλ(β(k))]-1Rp(β(k))

這里,D為p階導數(shù)陣?(Rp(β))/?β.

Δλ(β)=diag(qλ(|β1|)/(δ+|β1|),…,

qλ(|βp|)/(δ+|βp|))

式中,δ是一個非常小的正數(shù),文中取值為10-4.

4 數(shù)值模擬

通過式(1)產(chǎn)生真實數(shù)據(jù),這里令β=(3,2.5,0,0,3,0,0,0)T,xi取自p維多元正態(tài)分布,即xi～N(0,Ip),這里p=8.隨機誤差ε1,ε2,…,εn相互獨立,且都服從標準正態(tài)分布.

為了全面考察本文所提方法的穩(wěn)健性,這里將對數(shù)據(jù)進行2種污染力度、3種污染方式來產(chǎn)生異常值.通過2種污染力度在數(shù)據(jù)中分別產(chǎn)生較大異常值和較小異常值.污染方式1是隨機選取5個樣本只對x進行污染,將原本xi變?yōu)閤i+5(或xi+2.5);污染方式2是隨機選取5個樣本只對y污染,將原本yi變?yōu)閥i+10(或yi+3);污染方式3是隨機選取3個樣本對x污染,再隨機選取2個樣本對y污染,且所選樣本不重復,將選取的xi變?yōu)閤i+5(或xi+2.5),yi變?yōu)閥i+10(或yi+3).

為了排除隨機性對該方法的影響,本文將重復模擬200次,采用8個指標來衡量各種方法的效果.首先用估計的偏差、方差和均方誤差來衡量估計的優(yōu)劣;另一方面,用取對率、擬合不足率、擬合過度率、取對個數(shù)和取錯個數(shù)來衡量變量選擇的優(yōu)劣.其中,取對率就是將β中原本為0的分量估計為0,且將β中原本非0的分量估計為非0的總次數(shù)在200次模擬中的比率.擬合過度率為將β中原本的非0估計為非0,原本為0的也估計為非0的總次數(shù)在200次中的比率.擬合不足率即為200次模擬中,估計結(jié)果將β中原本不為0的估計為0的總次數(shù)在200次中的比率.取對個數(shù)是200次模擬中正確估出0的個數(shù)的平均值,取錯個數(shù)是200次模擬中將原本非0估計為0的個數(shù)的平均值,取對個數(shù)最好為5,取錯個數(shù)最差為3.

SCAD方法是變量選擇中比較常見的懲罰方法,由文獻[2]可知,SCAD方法在穩(wěn)健變量選擇中具有較好的優(yōu)良性,因此本文主要模擬不同SCAD方法的穩(wěn)健效果.表1～3分別表示污染1,2,3情況下各種方法變量選擇的模擬結(jié)果,表中字體傾斜的數(shù)據(jù)表示污染力度較小的結(jié)果.SCAD-NR表示方程(4)中φ(x)=x且wi=1時的非穩(wěn)健方法;SCAD-R (Huber)表示文獻[15]中所提的基于Huber函數(shù)的SCAD穩(wěn)健方法;SCAD-R (t2)表示本文所提的基于t2函數(shù)的SCAD穩(wěn)健方法.通過與非穩(wěn)健方法的比較,凸顯出本文所提方法與文獻[15]中方法是穩(wěn)健的.在數(shù)據(jù)中存在異常值時,穩(wěn)健估計結(jié)果明顯好于非穩(wěn)健的.而通過本文方法與文獻[15]中方法比較,模擬結(jié)果證明本文所提的方法在變量選擇上的效果明顯優(yōu)于基于Huber函數(shù)的方法.

表1 污染1情況下變量選擇的模擬結(jié)果Tab.1 Simulation results of variable selection under pollution 1

表2 污染2情況下變量選擇的模擬結(jié)果Tab.2 Simulation results of variable selection under pollution 2

由表1～3可見,在不同污染力度和不同污染方式下,非穩(wěn)健方法選擇的效果總體來說都沒有穩(wěn)健方法好.非穩(wěn)健方法擬合過度率比兩種穩(wěn)健的方法都高,也就是說,非穩(wěn)健方法總會把β中原本為0的估計為非0.由表中模擬結(jié)果可見,基于t函數(shù)的穩(wěn)健方法取對率總是最高的.無論是哪種污染方式,本文所提的方法與非穩(wěn)健方法和文獻[15]中的方法相比,取對率都遠大于其他兩種方法,而且取對個數(shù)也明顯大于其他兩種方法.雖然在擬合不足率上看不出本文所提方法的優(yōu)勢,但是,這里幾種方法的擬合不足率都非常接近,而在正確率、擬合過度率和取對率方面,卻明顯可以看出本文所提方法的優(yōu)勢.可見,本文方法對異常值的抵抗力比文獻[15]中基于Huber函數(shù)的穩(wěn)健方法更強,大大減少異常值在模型估計中引起的偏差.

表3 污染3情況下變量選擇的模擬結(jié)果Tab.3 Simulation results of variable selection under pollution 3

表4是在2種污染力度下第1種污染方式對β中非0分量的估計結(jié)果,統(tǒng)計了各非0分量的偏差、方差和均方誤差.從該結(jié)果中可見,兩種穩(wěn)健方法在估計方面也明顯優(yōu)于非穩(wěn)健方法,本文方法此時較文獻[15]并沒有很明顯的優(yōu)勢.在偏差的第1個分量和方差的第2,3個分量甚至稍微變大,但均方誤差卻稍有優(yōu)勢.相比較來說,基于t2函數(shù)的均方誤差比Huber函數(shù)的稍小.

表4 污染1情況下估計的模擬結(jié)果Tab.4 Simulation results of estimation under pollution 1

綜合表1～4,在非零參數(shù)估計上本文方法與文獻[15]中方法相比具有較小的均方誤差,同時本文方法在變量選擇方面的穩(wěn)健性和優(yōu)勢較為突出,其對異常值的抵抗力比文獻[15]中方法更強.

5 結(jié) 論

本文在前人研究的基礎上提出基于t2函數(shù)的穩(wěn)健變量選擇方法,并與文獻[15]中穩(wěn)健方法進行比較,通過數(shù)值模擬來驗證本文方法的有效穩(wěn)健性.文中首先詳細敘述了關于該方法的模型理論,在穩(wěn)健估計方程中施加一個懲罰函數(shù),來達到期望的穩(wěn)健變量選擇效果.然后在第2部分進一步考察穩(wěn)健方程中的有界得分函數(shù),通過比較不同自由度的t函數(shù)與Huber函數(shù)的性質(zhì),初步判斷方程(4)中的有界得分函數(shù)的穩(wěn)健性.第3部分介紹了本文所用的算法,采用牛頓迭代法.最后通過第4部分的數(shù)值模擬來驗證前文中預計的穩(wěn)健性.模擬結(jié)果體現(xiàn)了t2函數(shù)在變量選擇方面的明顯優(yōu)勢,雖然參數(shù)估計的結(jié)果并不是明顯好于Huber函數(shù)的結(jié)果,但是通過取對率、擬合不足率、擬合過度率、取對個數(shù)和取錯個數(shù)體現(xiàn)的選擇結(jié)果,說明t2函數(shù)在變量選擇上的穩(wěn)健性優(yōu)于Huber函數(shù).

本文主要通過模擬來考察各種方法的優(yōu)劣.基于t函數(shù)的變量選擇方法的大樣本性質(zhì),以及將該方法應用到更復雜的縱向數(shù)據(jù)中,或應用到超高維的橫截面數(shù)據(jù)中,這些問題還有待進一步研究.

[1] TIBSHIRANI R.Regression shrinkage and selection via the lasso:a retrospective[J].Journal of the Royal Statistical Society,2011,73(3):273-282.

[2] FAN J Q,LI R Z.Variable selection via nonconcave penalized likelihood and its oracle properties[J].Journal of the American Statistical Association,2001,96(456):1348-1360.

[3] ZOU H.The adaptive lasso and its oracle properties[J].Journal of the American Statistical Association,2006,101(476):1418-1429.

[4] ZOU H,HASTIE H.Regularization and variable selection via the Elastic Net[J].Journal of the Royal Statistical Society,2005,67(2):301-320.

[6] HUBER P J.Robust estimation of a location parameter[J].The Annals of Mathematical Statistics,1964,35(1):73-101.

[7] HUBER P J.Robust regression:asymptotics,conjectures and Monte Carlo[J].The Annals of Statistics,1973,1(5):799-821.

[8] PORTNOY S,KOENKER R.The Gaussian hare and the Laplacian tortoise:computability of squared-error versus absolute-error estimators[J].Statistical Science,1997,12(4):279-300.

[9] GILONI A,PADBERG M.Alternative methods of linear regression[J].Mathematical and Computer Modelling,2002,35(3/4):361-374.

[10] GILONI A,PADBERG M.The finite sample breakdown point ofL1-regression[J].Journal on Optimization,2004,14(4):1028-1042.

[11] HE X M,SIMPSON D G,WANG G Y.Breakdown points of t-type regression estimators[J].Biometrika,2000,87(3):675-687.

[12] HE X M,FUNG W K,ZHU Z Y.Robust estimation in generalized partial linear models for clustered data[J].Journal of the American Statistical Association,2005,100(472):1176-1184.

[13] SINHA S K.Robust inference in generalized linear models for longitudinal data[J].The Canadian Journal of Statistics,2006,34(2):261-278.

[14] MARONNA R A,MARTIN R D,YOHAI V J.Robust statistics:theory and methods[M].Chichester,England:John Wiley and Sons,2006.

[15] 樊亞莉,徐群芳.穩(wěn)健的變量選擇方法及其應用[J].上海理工大學學報,2013,35(3):256-260.

[16] SCHUMANN D H.Robust variable selection[D].Carolina:North Carolina State University,2009.