王新軍，朱永忠，李佳蓓

（河海大學(xué) 理學(xué)院，南京 211100）

0 引言

考慮普通最小二乘回歸問(wèn)題：

其中y為解釋變量，y∈Rn，X為回歸矩陣，X∈Rn×p，其中p為預(yù)測(cè)變量數(shù)，n為樣本觀測(cè)數(shù)，‖?‖2是Rn中的歐幾里得空間并且β∈Θ?Rn。為了更好地提出問(wèn)題，假設(shè)p＜n，XTX非奇異，對(duì)于上述問(wèn)題，普通最小二乘估計(jì)（OLS）是最常見(jiàn)的解決方法。然而，在高維問(wèn)題中，普通最小二乘估計(jì)很難處理，一些微小的變化可能導(dǎo)致估計(jì)不穩(wěn)定。比如，觀測(cè)數(shù)據(jù)的微小變化就可能導(dǎo)致參數(shù)估計(jì)的極端變化。

1996年，Tibshirani提出了Lasso估計(jì)[1]，其本質(zhì)是通過(guò)L1懲罰對(duì)模型系數(shù)進(jìn)行壓縮，將沒(méi)有影響或影響較小的自變量的回歸系數(shù)自動(dòng)壓縮到零，從而不僅在一定程度上能消除多重共線性的影響，而且在對(duì)參數(shù)進(jìn)行估計(jì)的同時(shí)也實(shí)現(xiàn)了對(duì)變量的選擇。

對(duì)于如下有約束的最優(yōu)化問(wèn)題：

Lasso估計(jì)以方差小幅度的降低來(lái)?yè)Q取偏差的少量增加，從而使得均方誤差整體上降低。Lasso估計(jì)不僅避免了普通最小二乘估計(jì)在高維問(wèn)題中的不穩(wěn)定性，而且能夠達(dá)到變量選擇的目的。然而Lasso引進(jìn)的L1懲罰函數(shù)不可微[2]，從而導(dǎo)致估計(jì)的協(xié)方差陣不存在緊集，對(duì)于此問(wèn)題，需要采用重采樣方法解決。

2008年，Park T和Casella G提出了Bayesian Lasso方法(BLasso)[3]，在貝葉斯理論框架下，對(duì)回歸系數(shù)引進(jìn)如下的獨(dú)立Laplace條件先驗(yàn)：

發(fā)現(xiàn)其邊際后驗(yàn)眾數(shù)與Lasso給出的估計(jì)一致，彌補(bǔ)了Lasso的不足。上述σ2的分布如果不合適會(huì)導(dǎo)致最優(yōu)化求解有多個(gè)滿足條件的估計(jì)值，而不是唯一的，因此σ2的后驗(yàn)分布必須是單峰的。BLasso法避免了采用重抽樣方法對(duì)協(xié)方差陣估計(jì)，然而，其未解決正態(tài)誤差項(xiàng)在極大似然假設(shè)下的穩(wěn)健性。

在很多情況下（基因組學(xué)、計(jì)量經(jīng)濟(jì)學(xué)、金融數(shù)據(jù)等）經(jīng)常會(huì)遇到數(shù)據(jù)中存在異常值或者當(dāng)數(shù)據(jù)不是服從正態(tài)分布假設(shè)，而是服從重尾分布。在這些情況下，回歸系數(shù)的估計(jì)很有可能受到影響，使得模型預(yù)測(cè)精度受到影響。

為此，本文改進(jìn)BLasso方法，在BLasso的基礎(chǔ)之上引入兩組潛變量，這樣做的優(yōu)勢(shì)在于：一是允許模型存在誤差或者異常值，能夠模擬出常見(jiàn)的偏態(tài)和重尾分布，從而能夠確保模型的穩(wěn)健性，二是保持模型的共軛性，能夠更好地應(yīng)用Gibbs抽樣。在引進(jìn)潛變量之后，利用邊緣極大似然估計(jì)方法，對(duì)潛變量邊際優(yōu)化，從而得出未知參數(shù)的極大似然函數(shù)。本文選用層次模型的誤差項(xiàng)是逆伽馬正態(tài)混合分布的尺度參數(shù)，這與BLasso層次模型采用的服從t分布的誤差項(xiàng)相類似[4]，但是選用了不同的層次模型以及懲罰函數(shù)結(jié)構(gòu)，從而能夠在已知參數(shù)的似然Laplace先驗(yàn)分布的情況下實(shí)現(xiàn)更加穩(wěn)健的Gibbs抽樣。通過(guò)數(shù)據(jù)仿真發(fā)現(xiàn)，當(dāng)數(shù)據(jù)中出現(xiàn)異常值，或者數(shù)據(jù)服從重尾分布的情況下能夠給出相比Lasso，BLasso更加精確的估計(jì)以及較好的變量選擇。

1 模型及方法

1.1 分層模型

為了使估計(jì)過(guò)程更穩(wěn)健，考慮如下回歸問(wèn)題，并假設(shè)貝葉斯方法是其解決方法：

上述最優(yōu)化問(wèn)題通過(guò)引入懲罰系數(shù)λ對(duì)損失函數(shù)L1以及回歸系數(shù)懲罰，其中λ∈[ ]0,1。之所以在損失函數(shù)L1中加入懲罰項(xiàng)，是考慮到保持模型預(yù)測(cè)和數(shù)據(jù)壓縮之間的平衡，從而獲得較好的穩(wěn)健估計(jì)。當(dāng)λ=0時(shí)，；當(dāng)λ=1時(shí)，

通過(guò)邊際似然方法，在假設(shè)兩組超參數(shù)以及懲罰參數(shù)的先驗(yàn)分布為L(zhǎng)aplace先驗(yàn)分布的情況下，得到服從邊際似然Laplace先驗(yàn)分布的重尾數(shù)據(jù)，利用貝葉斯理論得到各參數(shù)的后驗(yàn)分布，將這種方法稱之為平衡的貝葉斯Lasso（BBLasso）。

首先，為了給每一個(gè)參數(shù)構(gòu)建一個(gè)具有封閉形式的條件后驗(yàn)分布，所有參數(shù)采用如下的全條件分層結(jié)構(gòu)：

1.2 Gibbs抽樣

考慮到在所有參數(shù)的先驗(yàn)分布已知的情況下，各參數(shù)的全條件分布容易抽樣，于是可通過(guò)Gibbs抽樣來(lái)推導(dǎo)出各參數(shù)的全條件后驗(yàn)分布。通過(guò)邊緣優(yōu)化，去除含有潛變量τ和γ相關(guān)的表達(dá)式，得到如下的初始邊緣貝葉斯模型：

最后推導(dǎo)出β的全條件后驗(yàn)分布服從均值為y,方差為σ2A-1的多元正態(tài)分布，其中A=XT。潛變量和的倒數(shù)的全條件后驗(yàn)分布服從逆高斯分布，如下是的逆高斯參數(shù)：

利用共軛性推導(dǎo)出σ2的全條件后驗(yàn)分布服從參數(shù)分別為的逆伽

馬分布。

2 懲罰參數(shù)選擇

在頻率學(xué)派框架下，懲罰參數(shù)的選擇一般采用交叉驗(yàn)證法進(jìn)行選擇，然而在貝葉斯理論框架下，可以采用全新的辦法，比如利用無(wú)信息先驗(yàn)分布的假設(shè)，從分布中模擬出參數(shù)的合理估計(jì)值。本文采用兩種不同的懲罰函數(shù)參數(shù)的選擇方法。一種方式通過(guò)邊緣極大似然方法，在基于Gibbs連續(xù)抽樣的MCEM算法[5]上不斷地更新懲罰參數(shù)λ，另一種則是通過(guò)給定懲罰參數(shù)λ合適的先驗(yàn)信息，并且利用MWG（Metropolis within G）算法[6]從其后驗(yàn)分布中不斷抽樣，從而達(dá)到更新懲罰參數(shù)λ的目的。

2.1 經(jīng)驗(yàn)貝葉斯估計(jì)

將參數(shù)β，τ，γ和λ視為缺失數(shù)據(jù)，通過(guò)邊緣極大似然來(lái)估計(jì)懲罰參數(shù)λ，然后通過(guò)EM迭代算法來(lái)更新懲罰參數(shù)，由于邊緣極大似然方法涉及到積分運(yùn)算，當(dāng)模型維數(shù)很高時(shí)，很難直接通過(guò)積分運(yùn)算出λ的解析形式。為了解決此問(wèn)題，采用MCEM算法來(lái)估計(jì)λ的值。該算法的E-步均采用吉布Gibbs采樣，根據(jù)M-步可以導(dǎo)出λ(k)的估計(jì)值是如下三次多項(xiàng)式的解：

2.2 分層貝葉斯估計(jì)

為了得到所有參數(shù)的全貝葉斯模型，本文假設(shè)λ的先驗(yàn)分布為π(λ)=1，此方法能夠在缺少懲罰參數(shù)信息時(shí)保證λ∈(0 ,1)，并且使得方差部分的后驗(yàn)估計(jì)為通常的極大似然估計(jì)。使用上述先驗(yàn)時(shí)，很難計(jì)算出λ的全條件后驗(yàn)分布，但是通過(guò)計(jì)算推導(dǎo)出λ的全條件后驗(yàn)分布與如下形式的表達(dá)式成比例：

3 數(shù)據(jù)仿真

通過(guò)在不同的仿真情況下比較BBALASSO、LASSO、BLASSO的預(yù)測(cè)精度和變量選擇。LASSO選用Efron提出的5層交叉檢驗(yàn)的LARS算法[2]，BLASSO選用Gramacy的BLASSO算法來(lái)選擇壓縮系數(shù)[6]。利用R軟件來(lái)實(shí)現(xiàn)BBLASSO的Gibbs采樣。

選用如下線性模型：

在不同的預(yù)測(cè)值和誤差項(xiàng)分布情形下進(jìn)行仿真。每一個(gè)仿真數(shù)據(jù)集被分為一個(gè)驗(yàn)證集、一個(gè)訓(xùn)練集和一個(gè)測(cè)試集。從驗(yàn)證集中選擇LASSO的懲罰參數(shù)λj。在選擇完懲罰參數(shù)之后，混合驗(yàn)證集和訓(xùn)練集來(lái)估計(jì)LASSO的β。

仿真一：本例中選定系數(shù)β=(3,1.5,2,0,0,0,0,0)′，訓(xùn)練集中含有20組觀測(cè)值，驗(yàn)證集中含有20組觀測(cè)值，測(cè)試集中含有200組觀測(cè)值。設(shè)計(jì)矩陣X生成自均值為0，方差為1的多元正態(tài)分布，對(duì)于所有的i和j，向量xij和xjk的相關(guān)系數(shù)為0.5|j-k|。這里誤差項(xiàng)ε生成自均值為0，方差分別為1，3，5的多元正態(tài)分布[7]。

仿 真 二 ：本 例 中 選 定 系 數(shù)β=(5 .6,5.6,5.6,0)′，j＜k＜4時(shí)，cov (xij,xik)=-0.39，j＜4時(shí)，cov(xij,xi4)=0.23。其他參數(shù)設(shè)置與仿真一一致，包括數(shù)據(jù)集設(shè)置[8,9]。

仿真三：本例中選定系數(shù)：

三組相互獨(dú)立的預(yù)測(cè)因子P1，P2和P3生成自均值為0，方差為1的多元正態(tài)分布。令Xi=P1+εi，i=1,…,5；這里誤差項(xiàng)ε生成自均值為0，方差分別為1，3，5，10，15的多元正態(tài)分布，訓(xùn)練集中含有100組觀測(cè)值，驗(yàn)證集中含有100組觀測(cè)值，測(cè)試集中含有400組觀測(cè)值。

為了模擬出異常值，本文通過(guò)將生成的回歸矩陣X中的任意一個(gè)元素值修改成較大的值，從而可將模型視作存在異常值模型。將仿真一、仿真二及仿真三在不同異常值個(gè)數(shù)的詳細(xì)模擬結(jié)果分別列于表1、表2及表3中，模擬過(guò)程中均選定λ=0.5，迭代10000次，每組比較都將給出不同的預(yù)測(cè)精度和變量選擇效果。

表1 仿真一中四種方法重復(fù)50次的預(yù)測(cè)精度平均值

3.1 預(yù)測(cè)精度

對(duì)四種方法重復(fù)50次仿真模擬來(lái)比較各個(gè)方法的預(yù)測(cè)精度和變量選擇。利用各仿真模擬的預(yù)測(cè)均方誤差（PSEs）平均值來(lái)判別各個(gè)方法的預(yù)測(cè)精度。

表2 仿真二中四種方法重復(fù)50次的預(yù)測(cè)精度平均值

表3 仿真三中四種方法重復(fù)50次的預(yù)測(cè)精度平均值

從表1中可以看出在模型中不論是否存在異常值，Lasso預(yù)測(cè)精度相對(duì)其他三種方法較差。當(dāng)模型中存在異常值時(shí)，BBlasso預(yù)測(cè)精度明顯高于Blasso預(yù)測(cè)精度，并且可以明顯看出在存在異常值模型中，BBlasso的PSEs值幾乎一致以及隨著方差以及異常值個(gè)數(shù)的增加，均方誤差也隨著增加，從而可以印證BBlasso在存在異常值乃至極端異常值的模型中是比較穩(wěn)健的。

從表2中可以看出仿真二的結(jié)果與仿真一幾乎一致，結(jié)合仿真1與仿真2的模型可以看出BBlasso的兩種方法在小樣本乃至中等樣本中的表現(xiàn)相對(duì)其他方法有明顯優(yōu)勢(shì)。然而在大樣本模型中，就不會(huì)存在著明顯的優(yōu)勢(shì)，甚至?xí)鄬?duì)其他方法較差。

從表3中可以看出，當(dāng)模型中不存在異常值時(shí)，BBlasso-em的預(yù)測(cè)精度與Lasso以及Blasso相差不大，然而當(dāng)存在異常值時(shí)，BBlasso-em相比其他三種方法的預(yù)測(cè)精度有明顯的優(yōu)勢(shì)，BBlasso-prior相對(duì)其他三種方法，在預(yù)測(cè)精度上表現(xiàn)較差。

3.2 變量選擇

在多元線性回歸問(wèn)題中，模型變量選擇是其重要環(huán)節(jié)，本文為了比較模型變量選擇的好壞，采用ROC（Receiver Operating Characteristic）分析對(duì)四種方法的變量選擇進(jìn)行評(píng)判，并通過(guò)R軟件計(jì)算AUC值（Area Under Roc Curve）直觀地體現(xiàn)各個(gè)方法變量選擇的效果。AUC值就是處于ROC曲線下方的那部分面積的大小。通常，AUC的值介于0.5到1之間，較大的AUC代表了在變量選擇上表現(xiàn)得更好。

這里僅計(jì)算了仿真一以及仿真三各方法的AUC值，之所以不計(jì)算仿真二的AUC值在于仿真二僅有四個(gè)預(yù)測(cè)因子其中還包括了一個(gè)零，計(jì)算出來(lái)的AUC值不具代表性。對(duì)于仿真一以及仿真三的AUC值見(jiàn)表4和表5。

表4 仿真一中四種方法的AUC值

表5 仿真三中四種方法的AUC值

從表4可以看出，在仿真一中，BBlasso-em在模型不存在異常值的情況下，其AUC值與Blasso以及Lasso的AUC值幾乎一致，然而當(dāng)模型存在異常值時(shí),BBlasso-em的AUC值相比Blasso以及Lasso的AUC值有一定的優(yōu)勢(shì)。

從表5可以看出，在仿真三中，Blasso的AUC值均比其他方法的AUC值大，然而不得不考慮Blasso方法在未能很好地處理異常值的情況下，所得的AUC值與BBlasso-prior的AUC值相差不大。因此也可以說(shuō)明BBlasso在很好地處理異常值的情況下，變量選擇效果與其他方法相接近。

4 實(shí)例分析

4.1 前列腺數(shù)據(jù)

數(shù)據(jù)來(lái)源于Stamey 等[10]的研究，并在文獻(xiàn)[7,9，11]中多次被作為歷史數(shù)據(jù)使用。數(shù)據(jù)包含97例前列腺癌病例，感興趣結(jié)果變量為PSA水平，8個(gè)預(yù)測(cè)變量分別為：前列腺癌體積、前列腺重量系數(shù)、年齡、良性前列腺增生量、精囊入侵、包膜穿透、前列腺評(píng)分、前列腺評(píng)分4/5所占百分比。本文將67組前列腺數(shù)據(jù)作為訓(xùn)練集，30組前列腺數(shù)據(jù)作為測(cè)試集。為了判斷誤差項(xiàng)是否服從正態(tài)分布，本文通過(guò)K-S檢驗(yàn)，得出K-S檢驗(yàn)P值為0.7012，其遠(yuǎn)遠(yuǎn)大于0.05，可以認(rèn)為誤差項(xiàng)服從正態(tài)分布。此外，還通過(guò)繪制誤差項(xiàng)Q-Q圖，見(jiàn)圖1，可以看出誤差項(xiàng)與正態(tài)分布并沒(méi)有很大的偏差。

圖1 前列腺數(shù)據(jù)的誤差項(xiàng)Q-Q圖

通過(guò)前列腺數(shù)據(jù)，研究當(dāng)誤差項(xiàng)服從正態(tài)分布的情況下，Lasso、Blasso、Blasso-em、blasso-prior四種方法的預(yù)測(cè)精度，利用R軟件，最終得出四種方法的預(yù)測(cè)均方誤差值(PSE)，見(jiàn)表6。

表6 前列腺數(shù)據(jù)在四種方法下的PSE值

從表6中可以看出，當(dāng)數(shù)據(jù)誤差項(xiàng)服從正態(tài)分布時(shí)，四種方法的預(yù)測(cè)精度相差不大。

4.2 上證指數(shù)與宏觀經(jīng)濟(jì)指標(biāo)數(shù)據(jù)

隨著股市不斷發(fā)展，股市與宏觀經(jīng)濟(jì)之間的密切關(guān)系完全展現(xiàn)出來(lái)，為了研究股市與宏觀經(jīng)濟(jì)之間的關(guān)系[12]，本文選取2007年9月1日到2016年9月15日時(shí)間段上證指數(shù)共3300組數(shù)據(jù)，每日收盤價(jià)Pt作為基礎(chǔ)數(shù)據(jù)，將股指收益作為感興趣結(jié)果變量。一般地,收益率(也稱為對(duì)數(shù)收益率)序列定義為：

1,2,…,其中,P0,Pt分別為初始時(shí)刻及t時(shí)刻的資產(chǎn)價(jià)格。以8個(gè)宏觀經(jīng)濟(jì)指標(biāo)作為預(yù)測(cè)變量，宏觀經(jīng)濟(jì)指標(biāo)選用工業(yè)增加值、貨幣供應(yīng)量、年利率、人民幣儲(chǔ)蓄存款、國(guó)內(nèi)生產(chǎn)總值、股票成交金額、固定資產(chǎn)投資、進(jìn)出口差額[13]。本文將2300組數(shù)據(jù)作為訓(xùn)練集，1300組數(shù)據(jù)作為測(cè)試集。為了判斷誤差項(xiàng)是否服從正態(tài)分布，通過(guò)K-S檢驗(yàn)，得出K-S檢驗(yàn)P值遠(yuǎn)遠(yuǎn)小于0.05，可以認(rèn)為誤差項(xiàng)不服從正態(tài)分布，并通過(guò)繪制誤差項(xiàng)Q-Q圖，見(jiàn)圖2，可以看出誤差項(xiàng)Q-Q圖呈現(xiàn)出重尾現(xiàn)象。

圖2 上證指數(shù)與宏觀經(jīng)濟(jì)指標(biāo)數(shù)據(jù)的誤差項(xiàng)Q-Q圖

通過(guò)上述數(shù)據(jù)，研究當(dāng)誤差項(xiàng)服從重尾分布的情況下，Lasso、Blasso、Blasso-em、Blasso-prior四種方法的預(yù)測(cè)精度，利用R軟件，最終得出四種方法的預(yù)測(cè)均方誤差值（PSE），見(jiàn)下頁(yè)表7。

表7 上證指數(shù)與宏觀經(jīng)濟(jì)指標(biāo)數(shù)據(jù)在四種方法下的PSE值

從表7中可以看出，在數(shù)據(jù)誤差項(xiàng)服從重尾分布的情況下，LASSO與BLASSO預(yù)測(cè)精度相差不大，然而相比本文提出的BBLASSO兩種方法，還有較大差距。

5 總結(jié)

本文在Bayesian Lasso的基礎(chǔ)上，引入兩組潛變量，從而能夠允許模型存在誤差，通過(guò)選用未知參數(shù)以及潛變量合適的先驗(yàn)分布，提出分層模型，并通過(guò)Gibbs抽樣與邊際似然，獲得未知參數(shù)的后驗(yàn)分布，以及采用不同的懲罰參數(shù)結(jié)構(gòu)，來(lái)不斷地更新懲罰參數(shù)。從模擬研究和實(shí)例分析可以看出，本文提出的方法相比現(xiàn)有的方法，在模型中存在異常值以及重尾數(shù)據(jù)的情形下更加精確，并且也可以勝任變量選擇。

當(dāng)然本文也存在一些問(wèn)題，即假設(shè)回歸系數(shù)先驗(yàn)分布所依賴的條件參數(shù)對(duì)所有解釋變量都是相同的，也即對(duì)所有分量選用一樣的壓縮系數(shù)，為此可以進(jìn)一步研究對(duì)不同分量采用不同壓縮系數(shù)，從而研究出更加有效的估計(jì)方法。

參考文獻(xiàn)：

[1]Tibshirani R.Regression Shrinkage and Selection Via the Lasso[J].Journal of the Royal Statistical Society,1996,58(1).

[2]Efron B,Hastie T,Johnstone I,et al.Least angle Regression[J].Ann Statist,2004,32(2).

[3]Park T,Casella G.The Bayesian Lasso[J].J Amer Statist Assoc,2008,103(482).

[4]Gramacy R B,Pantaleio E.Shrinkage Regression for Multivariate Inference With Missing Data,and an Application to Portfolio Balancing[J].Bayesian Anal,2010，5(2).

[5]袁晶.貝葉斯方法在變量選擇問(wèn)題中的應(yīng)用[D].濟(jì)南：山東大學(xué)博士學(xué)位論文,2013.

[6]魏艷華,王丙參,孫永輝.分組數(shù)據(jù)場(chǎng)合Burr-XII分布參數(shù)貝葉斯估計(jì)的混合Gibbs算法[J].統(tǒng)計(jì)與決策,2015,(8).

[7]Li Q,Lin N.The Bayesian Elastic Net[J].Bayesian Anal,2010,5(1).

[8]Zou H.The Adaptive Lasso and Its Oracle Properties[J].J Amer Statist Assoc,2006,101(476).

[9]Leng C,Tran M N,Nott D.Bayesian Adaptive Lasso[J].Ann Inst Statist Math,2014,66(2).

[10]Stamey T A,Kabalin J N,McNeal J E,et al.Prostate Specific Antigen in the Diagnosis and Treatment of Adenocarcinoma of the Prostate.ii.Radical Prostatectomy Treated Patients[J].Journal of Urology，1989,141(5).

[11]Zou H,Hastie T.Regularization and Variable Selection via the Elastic net[J].J R Stat Soc Ser B Statist Methodol,2005,67(2).

[12]尚鵬岳,李勝宏.上證指數(shù)與宏觀經(jīng)濟(jì)指標(biāo)協(xié)整關(guān)系的實(shí)證分析[J].預(yù)測(cè),2002,(4).

[13]王曉燕,呂效國(guó),范玉,張媛.上證指數(shù)與宏觀經(jīng)濟(jì)指標(biāo)關(guān)系的實(shí)證分析[J].商業(yè)時(shí)代,2012,(10).