穩健的懲罰經驗似然方法及壓縮估計

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

1 引言

隨著科學技術的不斷發展，數據收集與處理的要求也日益提升.由于數據的來源之廣和數量的不斷積累，數據的結構越來越復雜，維數也越來越高.如何從海量的數據中快速高效的篩選并提取有價值的信息是我們值得研究和解決的問題.在統計學中，我們習慣的稱研究目標為響應變量，收集到的相關信息為解釋變量.變量之間的聯系具有普遍性，一般是借助統計模型推斷.線性回歸模型是研究變量之間相互關系的一種有效工具，但當統計模型中有較多的變量時，如果選擇的變量較少，就會發生變量的缺失導致模型會產生有偏差的估計和預測.為了避免這種困境，我們通常會選擇相對較多的自變量.但是自變量的過多又會產生模型的過度擬合而影響參數估計和預測的效率，因此變量選擇是統計模型和知識發現領域非常重要的課題之一.隨著高維數據分析需求的日益增加，這一問題的研究顯得更加突出而具有現實意義.

自20世紀六十年代以來，變量選擇一直是統計模型推斷的研究課題之一，許多文獻對其進行了詳細的研究.傳統的變量選擇有Akaike[1]在熵的概念上提供權衡估計模型復雜度和擬合數據優良性的標準AIC.從一組可供選擇的模型中選擇最佳模型時，通常選擇AIC值最小的變量子集.Schwarz[2]提出了BIC的變量選擇方法. BIC準則與AIC準則相似，都用于模型選擇，但在懲罰力度上，BIC的懲罰項比AIC的大.在樣本數量過多時，BIC準則可有效防止模型精度過高而造成的模型復雜度過高.Mallows[3]基于預測精度的角度下提出了CP準則，選擇使CP的值達到最小的模型作為最優模型.Hoerl和Kennard[4]提出了嶺回歸(ridge regression)的變量選擇方法.除了上述方法，常用的傳統變量選擇方法還有逐步回歸、最優子集回歸、嶺回歸等.許多研究表明，傳統的變量選擇方法總是由于某些客觀條件的限制，諸如協變量的維數很大等而表現的不盡如人意.近年來基于懲罰函數類型方法的出現在很大程度上使得變量選擇方法有較大的發展，開創了變量選擇方法的新景象.Tibshirani[5]提出了LASSO的變量選擇方法.Fan 和Li[6]提出了SCAD的變量選擇方法.Fran 和 Friedman[7]提出了橋回歸(bridge regression)的變量選擇方法.Zhang[8]發展了MCP的變量選擇方法.Yohai[9]研究了穩健的FPE準則的變量選擇.樊亞莉等[10]作出了縱向數據的穩健變量選擇方法.Welsh[11]開創了基于bootstrap分層的穩健變量選擇.Sophie 等[12]提出了結合懲罰項的穩健變量選擇.

經驗似然是構造置信區間和檢驗非參數的一門統計技術，類似于bootstrap的抽樣特性，這一方法與其它的統計方法相比有很多突出的優點.最早，Owen[13]在完全樣本下提出了一種進行非參數類型的統計推斷方法，并且將經驗似然方法的應用推廣到線性回歸模型的統計推斷中.Kolaczyk[14]把經驗似然方法應用到了廣義的線性模型.Chen[15-16]研究了非參數的回歸的經驗似然統計推斷.近些年來，關于經驗似然的研究也在不斷的發展中.Qin 和Lawless[17]研究了約束情況下參數的估計方程在經驗似然中的統計推斷.還有在復雜數據下所用到的一些經驗似然方法，如Wang 等[18]在刪失數據的一類生存函數中用到了經驗似然；Rao等[19]研究了帶缺失數據的線性模型的經驗似然統計推斷.

本文對普通的懲罰經驗似然方法提出了改進，詳細討論基于經驗似然方法的穩健變量選擇問題.具體來說，我們把普通經驗似然約束中的估計方程穩健化，提出了穩健經驗似然估計，再把懲罰經驗似然方法和穩健估計方程結合起來，提出了一種懲罰穩健經驗似然的壓縮估計方法，從而達到穩健估計和變量選擇同時進行.通過數值模擬分析，該新方法在變量的準確度和非零參數估計的均方誤差上優于普通的懲罰經驗似然估計.因此，本文的方法在變量選擇等數據分析中具有一定的優勢.

以下的內容安排如下.第2節介紹本文考慮的模型，以及基于懲罰穩健經驗似然的壓縮估計方法，第3節介紹本文的迭代算法和調節參數的選取方法，第4節進行數值模擬，第5節是實證分析，第6節為本文的總結.

2 模型和方法

2.1 線性回歸模型

考慮簡單的線性模型：

Y=XTβ+ε,

(2.1)

其中Xi=(Xi1,Xi2,…,Xip)T為解釋變量，i=1,2,…,n，Y=(Y1,Y2,…,Yn)T為響應變量，β=(β1,β2,…,βp)T為回歸系數，ε=(ε1,ε2,…,εn)T為隨機誤差，ε1,ε2,…,εn獨立同分布，且E(εi)=0，Var(εi)=σ2，i=1,2,…,n.待估計的參數為β，σ2.

2.2 普通的懲罰經驗似然估計

普通的最小二乘估計是基于以下估計方程得來的：

(2.2)

根據Owen[13，20]，模型(2.1)的普通經驗似然(EL)函數可表示為

滿足上式的β與pi，可以通過拉格朗日乘數法求解.令

(2.3)

普通的經驗似然函數為

對數經驗似然函數為

(2.4)

為了達到變量選擇的目的，我們在對數經驗似然函數的基礎上加懲罰項，得到如下的懲罰經驗似然函數

(2.5)

懲罰函數的選取有多種形式，如文獻[5]提出的LASSO，文獻[8]提出的MCP.本文選取文獻[6]提出的SCAD.SCAD懲罰函數具有無偏性、稀疏性、連續性.無偏性保證了系數較大的變量其估計系數是漸近無偏的，這樣避免了不必要的模型偏差.稀疏性保證了將不太重要的系數壓縮為零，從而起到選擇變量、降低模型復雜度的效果.連續性是減少了模型預測時的不穩定性，保證了數據本身的連續性.因此本文選取SCAD懲罰函數.

SCAD懲罰函數具有如下形式：

其中λ>0，a>2.根據文獻[6]，其一階導數為

2.3 穩健的懲罰經驗似然估計

(2.2)式所定義的估計方法對數據中的異常值是異常敏感的，因此，基于(2.2)式的普通懲罰經驗似然估計也會由于異常值的存在而有較大偏差.為了降低異常值的影響，把(2.2)式穩健化，采用有界的得分函數來限制估計方程中異常值的影響.記ψc(x)為對應于Huber函數的得分函數.函數ψc(x)的定義為

ψc(x)=min{c,max{-c,x}},

其中調節參數c用來調節估計的效率和穩健性.對于學生化殘差Xi,c的取值一般介于1到2之間.

與文獻[21]類似，我們采用如下的權重函數來限制杠桿點的影響：

記

其中σr為σ的絕對離差估計.

穩健經驗似然函數可表示為：

Qr(β)的最大值點就是β的穩健經驗似然估計.

加上懲罰函數的穩健經驗似然函數為

(2.6)

3 算法

3.1 基于SCAD的穩健懲罰經驗似然估計的牛頓迭代算法

由于SCAD懲罰函數是非凸且不可導，直接最小化目標函數并不容易.利用局部二次逼近方法，對SCAD懲罰函數在初始值β0展開，有

記

我們使用修改的牛頓迭代算法來執行懲罰經驗似然函數的優化.具體而言，對于k=0,1,2,…，我們生成迭代序列

其中

重復上述迭代過程，當|β(k+1)-β(k)|<ξ時就停止迭代.這里我們指定ξ=10-6.

3.2 調節參數的選取

為了選擇合適的調節參數，我們對普通的懲罰經驗似然估計和穩健的懲罰經驗似然估計都運用BIC準則來選取調節參數.為了簡單起見，在模擬時令所有的λ相等，即假設λ(j)=λ,j=1,2,…,p.我們在區間Ω=(0,λmax)中選擇最佳的λ.根據文獻[23]，

BICλ=log(RSS)+dfλlog(n)/n,

4 數值模擬

為了比較穩健的懲罰經驗似然方法及壓縮估計在有限樣本下的性質,我們通過數據模擬試驗來進行說明.

4.1 數據的產生

我們從模型(2.1)產生模擬數據，這里參數設定為β=(3,3,0,0,2.5,0,0,0)T，p=8，Xi=(Xi1,Xi2,…,Xip)T的每個分量服從N(0,1)，樣本大小n=50，對每組數據重復模擬N=200次.對于誤差項εi的分布，我們考慮三種情況：

(1)εi服從標準正態分布N(0,1)，且數據無污染；

(2)εi服從標準正態分布N(0,1)，但數據被污染；

(3)εi服從厚尾的t(2)分布.

4.2 數據的3種污染方式

污染1：僅對解釋變量X的污染, 隨機污染3%的Xi，并將其替換為Xi-1.2；

污染2：僅對響應變量Y的污染, 隨機污染3%的Yi，并將其替換為Yi-1.2；

污染3：解釋變量X和響應變量Y同時污染, 隨機污染3%的Xi，并將其替換為Xi-1.2，隨機污染3%的Yi，并將其替換為Yi-1.2.

數值模擬的結果如下面各表所示.

表1 正態分布時數據沒有污染時的模擬結果

由表1可知，在模擬數據無污染的情況下，本文所提出的穩健懲罰經驗似然方法及壓縮估計在非零參數估計的變量選擇準確度上相同于最小二乘估計和普通的懲罰經驗似然估計，但在參數為零的估計上明顯優于最小二乘估計，在非零參數的均方誤差下優于普通的懲罰經驗似然估計.

表2 正態分布時數據在污染1下的模擬結果

表3 正態分布時數據在污染2下的模擬結果

由表2和表3可知，無論是在解釋變量污染，還是在響應變量污染的情況下，本文所提出的方法在非零參數估計的變量選擇準確度上相同于最小二乘估計和普通的懲罰經驗似然估計，但在參數為零的估計上明顯優于最小二乘估計，在非零參數的均方誤差下優于普通的懲罰經驗似然估計.

表4 正態分布時數據在污染3下的模擬結果

由表4可知，即使在解釋變量和相應變量都污染的情況下，本文所提出的方法在非零參數估計的變量選擇準確度上相同于最小二乘估計和普通的懲罰經驗似然估計，但在參數為零的估計上明顯優于最小二乘估計，在非零參數的均方誤差下優于普通的懲罰經驗似然估計.

表5 t(2)分布時數據的模擬結果

由表5可知，不管是模擬數據在正態分布的情況下，還是在自由度為2的厚尾t分布情況下，本文所提出的方法在非零參數估計的變量選擇準確度上仍然相同于最小二乘估計和普通的懲罰經驗似然估計，但在參數為零的估計上還是明顯優于最小二乘估計，在非零參數的均方誤差下優于普通的懲罰經驗似然估計.

在模擬研究中的變量選擇效果上，本文提出的方法相較于普通的懲罰經驗似然估計的優勢似乎不太明顯，這可能和參數β=(3,3,0,0,2.5,0,0,0)T的設置有關，因為我們取了非零的值比較大，信號比較明顯.

5 實證分析

在本小節中，我們將穩健的懲罰經驗似然方法及壓縮估計應用到一組實際數據： Boston房屋數據.Harrison等[24]用這組數據在各種方法下研究了人們對于干凈空氣的需求，這組數據來源于網址:http://t.cn/RfHTAgY.該數據的每個類的觀察值是均等的，共有506個觀察值，13個解釋變量和1個響應變量，其中解釋變量包括：城鎮人犯罪率(CRIM)，住宅用地超過25000 sq.ft.的比例(ZN)，城鎮非零售商用土地的比例(INDUS)，查理斯河空變量(如果邊界是河流取1，否則取0)(CHAS)，一氧化氮濃度(NOX)，住宅平均房間數(RM)，1940年之前建成的自用房屋比例(AGE)，到波士頓五個中心區域的加權距離(DIS)，輻射性公路的接近指數(RAD)，每10000美元的全值財產稅率(TAX)，城鎮師生比例(PTRATIO),1000(BK-0.63)^2(其中BK指代城鎮中黑人的比例)(B),人口中地位地下者的比例(LSTAT).響應變量：自住房的平均房價，以千美元計(MEDV).考慮到實際的數據各項指標單位不同的數值差異較大，因此先將各項指標標準化，使其均值為0，方差為1.

為了比較，我們將這組數據做了交叉驗證來比較穩健的方法和普通的非穩健方法的表現.每次剔除1個觀察值，然后用剩下的505個觀察值的回歸系數，并用模型在交叉驗證過程中的均方誤差(MSEcv)來衡量各種方法的優劣：

表6 穩健與非穩健方法對自住房平均房價的MSEcv和Varcv 的估計比較

MSEcv值越小，表示該變量選擇的效果越好.在表6中可以看出最小二乘估計的MSEcv最小，但最小二乘估計不做參數壓縮和變量選擇，它把所有的參數值估計為非零，而穩健的方法和非穩健的方法都做變量選擇，并且做壓縮估計.綜合來看，本文提出的穩健懲罰方法最好.

6 結論

本文提出穩健的懲罰經驗似然方法及壓縮估計，給出了相應的算法，并作了大量的模擬研究.根據數值模擬研究表明，在數據無污染時，我們提出的方法在衡量變量選擇的優劣中具有較高的準確度，在非零分量估計的均方誤差下相對較??；在數據污染的情況下，本文的方法在選擇準確性和估計的均方誤差方面的優勢尤為顯著.所以，本文的穩健壓縮估計方法在變量選擇中具有一定可行性的優勢.

本論文穩健的懲罰經驗似然方法及壓縮估計可以推廣到縱向數據，也可以結合其他的懲罰函數，如LASSO[5]的變量選擇方法，MCP[8]的變量選擇方法等與SCAD的變量選擇方法來做比較.當然，證明本文方法的理論性質,如變量選擇的相合性，以及壓縮估計的漸近正態性，還需要做進一步研究.