劉亞新

（中南大學(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ù)估計的最小化問題等價于誤差項(xiàng)服從非對稱拉普拉斯分布的似然函數(shù)的最大化問題，進(jìn)而采用貝葉斯方法估計分位數(shù)回歸的系數(shù)，該方法在計量經(jīng)濟(jì)學(xué)領(lǐng)域受到了越來越多的關(guān)注。由于似然函數(shù)很復(fù)雜，參數(shù)的后驗(yàn)分布不是熟悉的函數(shù)形式，因此常常用MCMC方法對參數(shù)的后驗(yàn)分布進(jìn)行抽樣模擬，得到參數(shù)的貝葉斯估計。Yu和Moyeed采用的是隨機(jī)游走M(jìn)-H算法，建議分布為以當(dāng)前參數(shù)值為均值的正態(tài)分布。盡管該方法非常方便地產(chǎn)生候選值，但是接受率依賴于分位數(shù)，而且收斂速度很慢。Kozumi和Koboyashi[3]提出了基于非對稱拉普拉斯分布的位移-尺度模型的Gibbs抽樣算法。該算法根據(jù)參數(shù)的全條件后驗(yàn)分布進(jìn)行抽樣，而全條件后驗(yà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ù)的先驗(yàn)分布設(shè)成拉普拉斯先驗(yàn)，并用Gibbs算法進(jìn)行抽樣，實(shí)驗(yàn)證明貝葉斯Lasso分位數(shù)回歸比其他方法的Lasso分位數(shù)回歸更優(yōu)。Alhamzawi等[10]提出了貝葉斯自適應(yīng)Lasso分位數(shù)回歸，對不同的回歸系數(shù)賦予不同的懲罰參數(shù)，且懲罰參數(shù)的先驗(yàn)設(shè)為倒伽瑪分布，并把倒伽瑪分布的超參數(shù)設(shè)成未知量，采用Gibbs算法對后驗(yàn)分布抽樣來估計參數(shù)。

本文對帶有Elastic Net懲罰的貝葉斯分位數(shù)回歸提出了更簡單的估計方法，使所有參數(shù)的全條件后驗(yàn)分布都是熟知的分布形式，可以采用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是誤差項(xiàng)，且相互獨(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è)誤差項(xiàng)ε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ù)的全條件后驗(yàn)分布，將非對稱拉普拉斯分布用指數(shù)分布和正態(tài)分布混合表示。Kozumi和Kobayashi[11]已證明，若z服從指數(shù)分布exp(τ(1-τ)σ)，v服從標(biāo)準(zhǔn)正態(tài)分布N(0,1)，誤差項(xiàng)可以表示為:

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

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

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

其中Kc是第二類修正貝塞爾函數(shù)，d＞0，f＞0，c是實(shí)數(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)合后驗(yàn)密度為：

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

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

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

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

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

由各參數(shù)的全條件后驗(yàn)分布可以采用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) ，誤差項(xiàng)ε~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ù)測方面也具有較小的偏差。