帶有Elastic Net懲罰的貝葉斯分位數(shù)回歸及其數(shù)值模擬

劉亞新

（中南大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，長沙 410083）

0 引言

自1978年Koenker和Bassett[1]提出分位數(shù)回歸的概念以來，分位數(shù)回歸因能提供更全面的信息以及優(yōu)良的性質(zhì)，在理論和應(yīng)用領(lǐng)域都得到了廣泛的研究和應(yīng)用。Yu和Moyeed[2]最早提出了貝葉斯分位數(shù)回歸，將分位數(shù)回歸問題與非對稱拉普拉斯分布（Asymmetric Laplace Distribution,ALD）聯(lián)系起來，分位數(shù)回歸系數(shù)估計的最小化問題等價于誤差項服從非對稱拉普拉斯分布的似然函數(shù)的最大化問題，進(jìn)而采用貝葉斯方法估計分位數(shù)回歸的系數(shù)，該方法在計量經(jīng)濟(jì)學(xué)領(lǐng)域受到了越來越多的關(guān)注。由于似然函數(shù)很復(fù)雜，參數(shù)的后驗分布不是熟悉的函數(shù)形式，因此常常用MCMC方法對參數(shù)的后驗分布進(jìn)行抽樣模擬，得到參數(shù)的貝葉斯估計。Yu和Moyeed采用的是隨機(jī)游走M(jìn)-H算法，建議分布為以當(dāng)前參數(shù)值為均值的正態(tài)分布。盡管該方法非常方便地產(chǎn)生候選值，但是接受率依賴于分位數(shù)，而且收斂速度很慢。Kozumi和Koboyashi[3]提出了基于非對稱拉普拉斯分布的位移-尺度模型的Gibbs抽樣算法。該算法根據(jù)參數(shù)的全條件后驗分布進(jìn)行抽樣，而全條件后驗分布是已知的函數(shù)形式，這大大提高了抽樣的收斂速度，因此得到了廣泛應(yīng)用。

對分位數(shù)回歸中的變量選擇問題，Koenker[4]首次將Lasso的思想應(yīng)用于分位數(shù)回歸中。Wang等[5]將最小絕對差估計與自適應(yīng)Lasso懲罰結(jié)合起來進(jìn)行變量選擇。Li和Zhu[6]將Lasso的思想應(yīng)用于分位數(shù)回歸中進(jìn)行變量選擇，將系數(shù)的絕對值和作為懲罰部分，提出了計算Lasso分位數(shù)回歸的全部正則化路徑的有效算法，同時又對擬合模型的有效維數(shù)進(jìn)行估計，可以用來選擇正則化參數(shù)。Wu和Liu[7]證明了SCAD方法和自適應(yīng)Lasso分位數(shù)回歸的Oracle性質(zhì)。Park和Casella[8]從貝葉斯角度研究Lasso分位數(shù)回歸問題，提出了分層模型，并用Gibbs抽樣進(jìn)行參數(shù)估計。Li等[9]從貝葉斯的角度研究Lasso分位數(shù)回歸，提出了一個分層模型的框架，將參數(shù)的先驗分布設(shè)成拉普拉斯先驗，并用Gibbs算法進(jìn)行抽樣，實驗證明貝葉斯Lasso分位數(shù)回歸比其他方法的Lasso分位數(shù)回歸更優(yōu)。Alhamzawi等[10]提出了貝葉斯自適應(yīng)Lasso分位數(shù)回歸，對不同的回歸系數(shù)賦予不同的懲罰參數(shù)，且懲罰參數(shù)的先驗設(shè)為倒伽瑪分布，并把倒伽瑪分布的超參數(shù)設(shè)成未知量，采用Gibbs算法對后驗分布抽樣來估計參數(shù)。

本文對帶有Elastic Net懲罰的貝葉斯分位數(shù)回歸提出了更簡單的估計方法，使所有參數(shù)的全條件后驗分布都是熟知的分布形式，可以采用Gibbs算法進(jìn)行抽樣，避免了Gibbs抽樣中又包含M-H抽樣的問題[9]，提高了迭代效率。

1 帶有Elastic Net懲罰的貝葉斯分位數(shù)回歸

1.1 分位數(shù)回歸

給定訓(xùn)練數(shù)據(jù){(xi,yi),i...,n},xi為預(yù)測變量，yi為響應(yīng)變量，預(yù)測變量xi與響應(yīng)變量yi滿足：

其中β=(β1，β2，...，βk)T∈Rk,εi是誤差項，且相互獨(dú)立，滿足εi的第τ分位數(shù)為0，即的τ分位數(shù)回歸方程可以寫為：

根據(jù)Koenker和Bassett[1]的理論，分位數(shù)回歸系數(shù)β的求解可轉(zhuǎn)化為下列優(yōu)化問題：

其中ρτ(u)=u(τ-I(u＜0))是損失函數(shù)，I(·)表示示性函數(shù)。

Yu和Moyeed[2]指出，在假設(shè)誤差服從非對稱拉普拉斯分布的前提下，最小化問題（3）可以轉(zhuǎn)化為最大化似然函數(shù)問題。即假設(shè)誤差項εi服從非對稱拉普拉斯分布ALD(0，σ，τ)，其中σ是尺度參數(shù)，εi的密度函數(shù)為：

可以證明，服從該分布的變量的τ分位數(shù)是0。進(jìn)而yi|xi，β，σ~ALD(，σ，τ)，密度函數(shù)為：

樣本集y=(y1，y2，...，yn)T的似然函數(shù)為：

需要說明的是，誤差ε并不是真的服從ALD分布，這樣的假設(shè)只是為了將最小化問題（3）轉(zhuǎn)化為最大化似然函數(shù)（4）。不管數(shù)據(jù)的原始分布是什么，使用非對稱拉普拉斯分布都是一種有效的擬合貝葉斯分位數(shù)回歸的方式，即使誤差不是真的服從非對稱拉普拉斯分布，參數(shù)估計仍能取得很好的效果。

1.2 基于Elastic Net懲罰的貝葉斯分層模型

帶有Elastic Net懲罰的分位數(shù)回歸的參數(shù)估計模型為：

由于：

式（5）可以轉(zhuǎn)化為:

其中X=(x1，x2，…，xn)T，Ik×k表示k階單位矩陣，則式（6）將變?yōu)?

令:

則:

可以看出，最小化式（7）相當(dāng)于最大化式（8）。令η=，為了得到各參數(shù)的全條件后驗分布，將非對稱拉普拉斯分布用指數(shù)分布和正態(tài)分布混合表示。Kozumi和Kobayashi[11]已證明，若z服從指數(shù)分布exp(τ(1-τ)σ)，v服從標(biāo)準(zhǔn)正態(tài)分布N(0,1)，誤差項可以表示為:

從而:yi|xi，zi，β，σ~N(+(1-2τ)zi，2σzi)

為了得到參數(shù)的貝葉斯估計，需要指定各參數(shù)的先驗分布。本文對尺度參數(shù)σ，懲罰參數(shù)(λ1，η)分別采用各自的共軛先驗分布，σ的先驗分布為倒伽瑪分布IG(a，b)，λ1的先驗分布為N(μ，δ2)I(λ1＞0)，η的先驗分布為廣義倒高斯分布GIG(c，d，f)，其概率密度函數(shù)為:

則式（1）將轉(zhuǎn)化為:

其中Kc是第二類修正貝塞爾函數(shù)，d＞0，f＞0，c是實數(shù)。

外墻保溫性能的高低，將會直接影響舊工業(yè)建筑改造后日常運(yùn)營中的能耗，利用原有厚磚墻的畜熱性能在此基礎(chǔ)上采取修復(fù)、增加保溫層等措施進(jìn)行外墻的改造設(shè)計。一般情況下，舊工業(yè)建筑的外窗氣密性和保溫性能差，應(yīng)進(jìn)行節(jié)能計算采取相應(yīng)的改造更換措施，將原單層平板玻璃更換位為中空低輻射節(jié)能玻璃。屋頂保溫性能設(shè)計。采取更新原來保溫層的辦法進(jìn)行改造設(shè)計。

綜上所述，得到如下的貝葉斯分層模型：

1.3 Gibbs抽樣

由貝葉斯分層模型可得，參數(shù) (β，σ，λ1，η，z)的聯(lián)合后驗密度為：

根據(jù)貝葉斯定理，回歸系數(shù)βj的全條件后驗密度為：

其中β-j表示去除參數(shù)βj之后的參數(shù)向量，xi，-j也類似，則βj的全條件后驗分布為正態(tài)分布N(ba′，1a′)。

參數(shù)σ的全條件后驗密度為：

故參數(shù)σ的全條件后驗分布為倒伽瑪分布，即：

同理可得其他參數(shù)的全條件后驗分布分別為：

由各參數(shù)的全條件后驗分布可以采用Gibbs算法進(jìn)行抽樣，其中廣義倒高斯分布的抽樣可以參考Dagpunar[12]和J?rgensen[13]，也可以借助 R軟件中GeneralizedHyperbolic程序包中的rgig函數(shù)。

2 數(shù)值模擬

2.1 獨(dú)立同分布情形

對于獨(dú)立同分布情形，本文采用在這一研究領(lǐng)域經(jīng)常使用的數(shù)值模擬模型：

其中回歸系數(shù)分為以下兩種形式：

Case1：β=(3，1.5，0，0，2，0，0，0)T

Case2：β=(0.8，0.8，0.8，0.8，0.8，0.8，0.8，0.8)T

Case1和Case2分別對應(yīng)稀疏和密集的情形。在Case1中，自變量的樣本數(shù)據(jù)產(chǎn)生于多元正態(tài)分布Nk( )0，Σx，誤差ε~N(0，1)。在Case2中，自變量的樣本數(shù)據(jù)產(chǎn)生于多元正態(tài)分布誤差ε~N(0，1)。模型生成20組數(shù)據(jù)集，每組數(shù)據(jù)集中有400條訓(xùn)練樣本，100條測試樣本。評價指標(biāo)為對測試樣本計算的均值絕對差的均值（MMAD）和標(biāo)準(zhǔn)差（SD），即：

分別估計分位數(shù)τ=(0.3，0.5，0.8)時的分位數(shù)回歸模型。

表1和表2（見下頁）分別是獨(dú)立同分布情形下Case1的參數(shù)估計結(jié)果和MMAD（SD）值。可以看出，在參數(shù)估計方面，帶Elastic Net懲罰的分位數(shù)回歸方法比QRLasso方法和QR-SCAD方法更準(zhǔn)確，但是其MMAD值和SD值比這兩種方法都大，而只比QR方法小，說明當(dāng)預(yù)測變量之間的相關(guān)性較低時，帶Elastic Net懲罰的分位數(shù)回歸方法的優(yōu)勢并不十分明顯。

表1 獨(dú)立同分布情形下Case1的參數(shù)估計

表2 Case1各方法的MMAD(SD)值

表3和表4分別是獨(dú)立同分布情形下Case2的參數(shù)估計結(jié)果和MMAD（SD）值。Case2對應(yīng)密集的情形，并且變量之間存在較高的相關(guān)性，可以看出帶Elastic Net懲罰的分位數(shù)回歸方法凸顯出了優(yōu)勢，其MMAD(SD)值比其他四種方法都要小，而且參數(shù)估計也相對更加準(zhǔn)確。這說明帶Elastic Net懲罰的分位數(shù)回歸方法在變量之間存在比較嚴(yán)重的多重共線性時具有明顯的效果，同時也具有較強(qiáng)的穩(wěn)健性。

表3 獨(dú)立同分布情形下Case2的參數(shù)估計

2.2 非獨(dú)立同分布情形

對于非獨(dú)立同分布下的數(shù)值模擬，本文采用如下回歸模型：

預(yù)測變量x1，x2，…，x7的樣本數(shù)據(jù)產(chǎn)生于均勻分布U(0 ，1) ，誤差項ε~N(0 ，1)。與獨(dú)立同分布情形類似，模型生成20個數(shù)據(jù)集，每個數(shù)據(jù)集中有400條訓(xùn)練樣本，100條測試樣本，估計分位數(shù)τ=0.2，0.3，0.5下的分位數(shù)回歸方程。

表5和表6（見下頁）分別是非獨(dú)立同分布下的參數(shù)估計結(jié)果和MMAD（SD）值。可以看出，在非獨(dú)立同分布下，五種方法仍能得到比較準(zhǔn)確的估計。其中，本文的方法在預(yù)測方面具有較小的標(biāo)準(zhǔn)差(SD)，總體上要優(yōu)于QR方法和QRLasso方法。Alasso方法雖然比帶Elastic Net懲罰的分位數(shù)回歸方法有更小的MMAD值和SD值，但是其不能解決預(yù)測變量個數(shù)大于樣本量這種情況，而且Alasso方法的懲罰參數(shù)多，抽樣復(fù)雜，沒有帶Elastic Net懲罰的分位數(shù)回歸方法便于理解，易于操作。

表4 Case2各方法的MMAD(SD)值

表5 非獨(dú)立同分布情形下的參數(shù)估計

表6 非獨(dú)立同分布情形下各方法的MMAD(SD)值

3 結(jié)論

帶有Elastic Net懲罰的參數(shù)估計方法是對Lasso方法的改進(jìn)，可以解決預(yù)測變量大于樣本量情況下的變量選擇問題，并且當(dāng)預(yù)測變量間存在著較高相關(guān)性時，Elastic Net懲罰仍然能得到比較準(zhǔn)確的參數(shù)估計。本文研究了帶有Elastic Net懲罰的貝葉斯分位數(shù)回歸問題，建立了相應(yīng)的貝葉斯分層模型，使Gibbs抽樣成為可能，提高了馬爾科夫鏈的收斂速度，同時在參數(shù)估計和預(yù)測方面也具有較小的偏差。