不同先驗(yàn)下分位數(shù)回歸的貝葉斯估計(jì)：實(shí)現(xiàn)與比較

邸俊鵬，張曉峒(.上海社會(huì)科學(xué)院 數(shù)量經(jīng)濟(jì)研究中心，上海 0000；.南開(kāi)大學(xué) 數(shù)量經(jīng)濟(jì)研究所，天津 30007)

一、引 言

傳統(tǒng)最小二乘回歸模型是建立在隨機(jī)變量的條件均值函數(shù)(Conditional Mean Function)基礎(chǔ)之上的，即用一組變量來(lái)解釋因變量的期望。然而，在經(jīng)濟(jì)學(xué)、金融學(xué)、醫(yī)學(xué)、環(huán)境科學(xué)等研究領(lǐng)域，隨機(jī)變量在任意概率水平下與解釋變量的關(guān)系越來(lái)越受到重視，即考察分位數(shù)回歸模型。分位數(shù)回歸模型估計(jì)方法一般包括三類(lèi)：第一類(lèi)是直接優(yōu)化法，如單純形法和內(nèi)點(diǎn)法；第二類(lèi)是參數(shù)化方法，該方法借助貝葉斯統(tǒng)計(jì)理論，在得到參數(shù)估計(jì)量分布的基礎(chǔ)上，構(gòu)建置信區(qū)間并進(jìn)行統(tǒng)計(jì)推斷；第三類(lèi)是與參數(shù)化相對(duì)的非(半)參數(shù)方法。本文聚焦分位數(shù)回歸的貝葉斯估計(jì)方法。

貝葉斯推斷在一般的線(xiàn)性或擴(kuò)展模型中已經(jīng)得到廣泛的應(yīng)用，尤其是針對(duì)復(fù)雜的目標(biāo)函數(shù)，可以通過(guò)蒙特卡羅模擬得到參數(shù)的后驗(yàn)概率分布。然而，在分位數(shù)回歸領(lǐng)域，有關(guān)貝葉斯推斷的文獻(xiàn)并不多，在2000年之前只有為數(shù)不多的文章對(duì)貝葉斯分位數(shù)回歸方法及其應(yīng)用進(jìn)行了初步嘗試[1]。Yu和Moyeed發(fā)現(xiàn)，可以借助曾在Koenker和Bassett的一文中使用的非對(duì)稱(chēng)拉普拉斯分布(Asymmetric Laplace Distribution，ALD)，將分位數(shù)回歸納入貝葉斯推斷的框架，從而實(shí)現(xiàn)了采用貝葉斯方法分析分位數(shù)回歸模型[2-3]。

Yu和Moyeed借助ALD可謂開(kāi)啟了貝葉斯分位數(shù)回歸的參數(shù)估計(jì)之門(mén)，之后在計(jì)量經(jīng)濟(jì)學(xué)領(lǐng)域中該方法越來(lái)越受到關(guān)注。Yu、Philippe和Jin采用該方法對(duì)英國(guó)工資分布進(jìn)行了研究[4]，在貝葉斯推斷中，他們視待估參數(shù)為隨機(jī)變量，使用MCMC方法得到參數(shù)的后驗(yàn)分布。通過(guò)研究發(fā)現(xiàn)，與傳統(tǒng)的方法相比貝葉斯分位數(shù)回歸方法有以下優(yōu)勢(shì)：一是盡管兩種方法得到的點(diǎn)估計(jì)值相近，但由于傳統(tǒng)的方法基于“漸近性質(zhì)”，所以貝葉斯方法相對(duì)于頻率方法得到的標(biāo)準(zhǔn)誤差更??；二是貝葉斯分位數(shù)回歸方法是關(guān)于因變量的全后驗(yàn)分布(Full Posterior Distributions)，而不是單個(gè)值，因而對(duì)估計(jì)參數(shù)的了解更全面；三是頻率學(xué)派的假設(shè)檢驗(yàn)是建立在對(duì)參數(shù)分布設(shè)定的基礎(chǔ)之上的，而貝葉斯分位數(shù)回歸是通過(guò)參數(shù)后驗(yàn)分布的HPD(Highest Posterior Density)進(jìn)行假設(shè)檢驗(yàn)，因而功效更高。

根據(jù)貝葉斯定理，后驗(yàn)密度正比于參數(shù)先驗(yàn)密度和樣本似然函數(shù)的乘積，因而在實(shí)施貝葉斯分位數(shù)回歸估計(jì)過(guò)程中，有兩項(xiàng)關(guān)鍵的工作：一是設(shè)定參數(shù)的先驗(yàn)分布；二是通過(guò)抽樣得到參數(shù)的后驗(yàn)分布。

參數(shù)的先驗(yàn)分布在貝葉斯分析中占有重要地位，它代表了研究者對(duì)參數(shù)的主觀信念，無(wú)論是在貝葉斯分位數(shù)回歸的參數(shù)估計(jì)還是半?yún)?shù)估計(jì)過(guò)程中，都涉及先驗(yàn)分布的設(shè)定。在先驗(yàn)分布的選擇上，Yu和Moyeed推導(dǎo)得出，如果設(shè)定參數(shù)的先驗(yàn)分布為均勻分布，即使這種選擇是不適當(dāng)?shù)?improper)，由此得到的參數(shù)的聯(lián)合后驗(yàn)分布也是適當(dāng)?shù)?proper)[2]。因此，在貝葉斯分位數(shù)回歸中，先驗(yàn)分布可以設(shè)定為無(wú)信息先驗(yàn)，如均勻分布或者方差很大的正態(tài)分布。在通常的情形下，給不同分位數(shù)下的估計(jì)參數(shù)設(shè)定統(tǒng)一的先驗(yàn)分布，然而Alhamzawi和Yu認(rèn)為實(shí)際研究中高分位和低分位下參數(shù)的分布可能是不同的，因而應(yīng)該采用先驗(yàn)勢(shì)(Power Prior)的方法針對(duì)不同的分位數(shù)設(shè)定不同的先驗(yàn)分布[5]。

當(dāng)模型參數(shù)的后驗(yàn)分布為熟悉的函數(shù)形式(如正態(tài)分布、Gamma分布、Beta 分布等)時(shí)，可以很容易地對(duì)相關(guān)參數(shù)的后驗(yàn)分布進(jìn)行抽樣模擬，但在一般情形下模型參數(shù)后驗(yàn)分布的形式是未知的，此時(shí)可以運(yùn)用以下幾類(lèi)算法對(duì)不熟悉(非標(biāo)準(zhǔn))的參數(shù)后驗(yàn)分布進(jìn)行抽樣模擬。常用的貝葉斯抽樣算法有：Hamiltonian Monte Carlo(HMC)、Gibbs算法、M-H算法、拒絕抽樣(Reject Sampling)和重要抽樣算法(Important Sampling)。

在貝葉斯分位數(shù)回歸方法中一些細(xì)節(jié)問(wèn)題也受到人們的關(guān)注，比如關(guān)于尺度參數(shù)σ。在一般的貝葉斯分位數(shù)回歸分析中都將σ設(shè)定為1，如Yu和Moyeed針對(duì)連續(xù)因變量[2]、Benoit和Poel針對(duì)離散因變量的模型[6]，而在一些研究中則認(rèn)為這種設(shè)定是存在問(wèn)題的，不符合金融時(shí)間序列的特征，建議將σ參數(shù)化，與其他參數(shù)一起估計(jì)[7-8]。

本文在已有文獻(xiàn)的基礎(chǔ)上，首先借助泊松分布實(shí)現(xiàn)了在STAN中對(duì)分位數(shù)回歸進(jìn)行貝葉斯分析，從而不僅可以靈活地選擇先驗(yàn)分布，還可以參數(shù)化規(guī)模參數(shù)，提高貝葉斯分位數(shù)回歸方法的靈活性和可操作性；然后通過(guò)模擬實(shí)驗(yàn)研究先驗(yàn)分布的設(shè)定對(duì)特定抽樣算法下參數(shù)估計(jì)量統(tǒng)計(jì)性質(zhì)的影響。

二、基于ALD的貝葉斯分位數(shù)回歸估計(jì)方法

設(shè)有如下線(xiàn)性回歸模型：

yt=Qτ(yt|xt；β)+ut

與普通最小二乘法的殘差平方和最小不同，分位數(shù)回歸是最小化損失函數(shù)：

(1)

其中，Qτ(yt|xt；β)為yt在xt；β條件下的τ分位數(shù)；β為參數(shù)向量；ρτ(u)=u(τ-I(u<0))為損失函數(shù)，其中I(·)為指示函數(shù)，當(dāng)滿(mǎn)足()內(nèi)條件時(shí)取值為1，否則為0。采用式(1)進(jìn)行分位數(shù)回歸的方法有多種，較常用的線(xiàn)性規(guī)劃解法是直接優(yōu)化分位數(shù)回歸的目標(biāo)函數(shù)，但直接優(yōu)化一般存在困難，而采用貝葉斯方法相對(duì)容易實(shí)現(xiàn)。用貝葉斯方法對(duì)分位數(shù)回歸進(jìn)行估計(jì)，其基本步驟是：首先，通過(guò)假設(shè)分位數(shù)回歸模型(1)的誤差項(xiàng)為非對(duì)稱(chēng)拉普拉斯分布，寫(xiě)出似然函數(shù)；其次，采用基于HMC的貝葉斯估計(jì)方法可以得到分位數(shù)回歸參數(shù)的抽樣分布，并對(duì)抽樣結(jié)果進(jìn)行收斂性檢驗(yàn)。由于對(duì)應(yīng)的似然函數(shù)式連續(xù)但不可導(dǎo)，所以對(duì)參數(shù)求導(dǎo)沒(méi)有解析解，在這種情況下采用HMC方法，對(duì)該后驗(yàn)密度進(jìn)行抽樣得到各個(gè)分位數(shù)上的參數(shù)β和尺度參數(shù)σ的抽樣分布。

三、貝葉斯分位數(shù)回歸方法的STAN實(shí)現(xiàn)

STAN是一套程序化的貝葉斯統(tǒng)計(jì)分析軟件，其操作方便靈活且對(duì)用戶(hù)免費(fèi)。在使用過(guò)程中用戶(hù)只需指定數(shù)據(jù)的似然函數(shù)和參數(shù)的先驗(yàn)分布，就可通過(guò)Hamiltonian Monte Carlo(HMC)抽樣得到參數(shù)的后驗(yàn)分布，這對(duì)沒(méi)有顯性表達(dá)式的貝葉斯后驗(yàn)分布尤其重要。STAN軟件中提供了正態(tài)分布、t分布、指數(shù)分布、泊松分布等多種常見(jiàn)的分布，用來(lái)指定似然函數(shù)和先驗(yàn)分布[9]，這些分布能夠滿(mǎn)足大部分的貝葉斯統(tǒng)計(jì)分析。然而，如果貝葉斯分析中的分布為非標(biāo)準(zhǔn)的分布，如本文中的非對(duì)稱(chēng)拉普拉斯分布在軟件中并未提供，則需要應(yīng)用一定的編程技巧和數(shù)學(xué)轉(zhuǎn)換以適應(yīng)于STAN程序的要求。下面介紹通過(guò)泊松分布給出非對(duì)稱(chēng)拉普拉斯分布的似然函數(shù)，從而實(shí)現(xiàn)用STAN軟件對(duì)貝葉斯分位數(shù)回歸的參數(shù)進(jìn)行估計(jì)。

(一)用泊松分布指定非對(duì)稱(chēng)拉普拉斯分布

設(shè)數(shù)據(jù)的對(duì)數(shù)密度函數(shù)為li=logf(xi|θ)，密度函數(shù)為exp(li)，所以似然函數(shù)為：

一般情況下，為了保證λ為正數(shù)，可以在負(fù)對(duì)數(shù)似然函數(shù)中加一個(gè)較大的數(shù)c=10 000。一般的STAN程序如下：

model{

c<-10000

for i in 1:n{

l[i]<-…… (指定對(duì)數(shù)似然)

lambda[i]<- -l[i]+c

h[i]<-0 (所有的h取0)

}

對(duì)于表達(dá)式y(tǒng)i=β0+β1x+u，假設(shè)誤差項(xiàng)u服從非對(duì)稱(chēng)拉普拉斯分布，ALD(μ，σ，τ)；一般情況下σ=1，待估系數(shù)為β0、β1，假定它們的先驗(yàn)分布是正態(tài)分布，進(jìn)而可以利用泊松分布指定非對(duì)稱(chēng)拉普拉斯分布的似然函數(shù)，程序如下：

model {

int h[n]；

real u[n]；

real lamda[n]；

for(i in 1:n) {

h[i]=0；

u[i]=y[i]-beta0-beta1*x[i]；

lambda[i]=log(p*(1-p))+(fabs(u[i])+(2*p-1)*u[i])/2+10000；

}

如果非對(duì)稱(chēng)拉普拉斯分布的尺度參數(shù)σ也是未知的，待估參數(shù)是β0、β1、σ，并假設(shè)σ的先驗(yàn)分布是自由度為3的卡方分布，具體的STAN程序如下：

model{

int h[n]；

real u[n]；

real lamda[n]；

sigma～chi_square(3)；

for(i in 1:n) {

h[i]=0

u[i]=y[i]-beta0-beta1*x[i]

lambda[i]=log(p*(1-p))+log(sigma)+(fabs(u[i]/sigma)+(2*p-1)*u[i]

/sigma)/2+10000；

}

在具體的估計(jì)過(guò)程中，首先將上述寫(xiě)好的程序保存為“***.stan”文檔，然后在R軟件中用“rstan”調(diào)用STAN程序，從而得到待估參數(shù)的后驗(yàn)分布以及相關(guān)的統(tǒng)計(jì)量。除了用泊松分布指定非對(duì)稱(chēng)拉普拉斯分布的似然函數(shù)，還可以用二項(xiàng)分布和負(fù)二項(xiàng)分布指定。

(二)尺度參數(shù)是否參數(shù)化對(duì)貝葉斯分位數(shù)回歸估計(jì)結(jié)果的影響

在已有的文獻(xiàn)中，非對(duì)稱(chēng)拉普拉斯分布的尺度參數(shù)在估計(jì)過(guò)程中通常被指定為常數(shù)1。王新宇、宋學(xué)鋒指出這種做法在實(shí)際應(yīng)用中是不妥當(dāng)?shù)?，并證明如果隨機(jī)變量x來(lái)自非對(duì)稱(chēng)拉普拉斯分布ALD(μ,σ,τ)，密度函數(shù)為f(μ,σ,τ)，其方差記為φ(σ,τ)，則當(dāng)σ=1時(shí)，x的方差φ(1,τ)≥8，即如果尺度參數(shù)固定為1，將導(dǎo)致非對(duì)稱(chēng)拉普拉斯分布隨機(jī)變量的方差不小于8[7]，這對(duì)于實(shí)際的金融時(shí)間序列數(shù)據(jù)顯然是不符合常理的，應(yīng)將其視作待估的參數(shù)。

上文已經(jīng)在STAN中實(shí)現(xiàn)了貝葉斯分位數(shù)回歸估計(jì)中尺度參數(shù)的參數(shù)化。下面將通過(guò)一個(gè)模擬實(shí)驗(yàn)，比較尺度參數(shù)是否參數(shù)化對(duì)估計(jì)結(jié)果產(chǎn)生影響。

1.實(shí)驗(yàn)設(shè)計(jì)與抽樣結(jié)果。數(shù)據(jù)生成過(guò)程為Yi=β0+β1x+εi，i=1,2,…，n，其中β0=1，β1=2，x～unif(n=n，min=0，max=10)。設(shè)定β0和β1的先驗(yàn)分布均為正態(tài)分布N(0，10)，尺度參數(shù)σ的先驗(yàn)分布為自由度等于3的卡方分布。采用HMC抽樣算法共模擬10 000次，為消除初值對(duì)抽樣分布的影響，去掉前5 000個(gè)抽樣值。

值得說(shuō)明的是，在統(tǒng)計(jì)推斷之前需要對(duì)抽樣分布進(jìn)行檢驗(yàn)，常用的檢驗(yàn)方法有兩種：一種是判斷馬爾科夫鏈?zhǔn)欠袷諗浚@可以通過(guò)診斷統(tǒng)計(jì)量判斷，也可以通過(guò)抽樣的軌跡直觀地判斷；另一種是考察抽樣值的自相關(guān)圖。如果連續(xù)兩次的迭代值之間相關(guān)性太強(qiáng)將導(dǎo)致抽樣不能探測(cè)到參數(shù)的整個(gè)空間，測(cè)量?jī)蓚€(gè)抽樣值相依性的標(biāo)準(zhǔn)算法是計(jì)算抽樣值相隔一定滯后期L的相關(guān)系數(shù)。

圖1分別給出樣本容量為100時(shí)，β0、β1、σ在0.5分位點(diǎn)下的HMC抽樣值軌跡和核密度圖，結(jié)果顯示三個(gè)參數(shù)的5 000個(gè)抽樣值序列在某一個(gè)值的上下波動(dòng)，因此抽樣構(gòu)成的馬爾科夫鏈均收斂。

圖2分別給出了樣本容量為100時(shí)，β0、β1在0.5分位點(diǎn)下抽樣值的自相關(guān)圖*篇幅所限，σ在0.5分位點(diǎn)下抽樣的自相關(guān)圖省略，其結(jié)果類(lèi)似。。從圖2中可以看出，隨著滯后期的增加，自相關(guān)系數(shù)趨向于零。

圖1 β0、β1、σ抽樣值軌跡和核密度圖

綜上，馬爾科夫鏈?zhǔn)諗壳易韵嚓P(guān)系數(shù)趨于零，所以在樣本容量為100，各參數(shù)的貝葉斯中位數(shù)回歸估計(jì)結(jié)果是可靠的。其他樣本容量和分位點(diǎn)下的馬爾科夫鏈均收斂，自相關(guān)系數(shù)也均趨于零。

圖2 β0、β1抽樣值的自相關(guān)圖

2.結(jié)論。比較不同樣本容量下STAN軟件得到β0、β1在σ參數(shù)化與非參數(shù)化下的估計(jì)結(jié)果*因篇幅所限，圖表省略。如有讀者需要可向作者索要。，可以發(fā)現(xiàn)：1)在一定樣本容量、特定分位點(diǎn)下，當(dāng)尺度變量σ被參數(shù)化時(shí)，β0和β1的抽樣標(biāo)準(zhǔn)差均小于σ非參數(shù)化時(shí)各自的標(biāo)準(zhǔn)差，這表明尺度變量的參數(shù)化可以提高被估系數(shù)的估計(jì)精度，從而提高待估系數(shù)假設(shè)檢驗(yàn)的準(zhǔn)確性。因此，在采用HMC抽樣方法進(jìn)行貝葉斯分位數(shù)回歸估計(jì)時(shí)應(yīng)該將尺度變量σ參數(shù)化，與待估的系數(shù)一起估計(jì)，而不應(yīng)該設(shè)定為常數(shù)1。2)給定樣本容量，各分位點(diǎn)下，尺度變量是否參數(shù)化對(duì)β0和β1抽樣均值與真值間的偏誤沒(méi)有明顯的影響。因此，如果只考慮被估系數(shù)對(duì)真值的偏離程度就不必將尺度變量參數(shù)化。

此外，還可以發(fā)現(xiàn)：1)給定樣本容量，無(wú)論在尺度變量σ被參數(shù)化還是未被參數(shù)化的情形下，β1的偏誤總是小于同一分位點(diǎn)下β0的偏誤，β1的標(biāo)準(zhǔn)差總是小于同一分位點(diǎn)下β0的標(biāo)準(zhǔn)差，這表明在本模擬實(shí)驗(yàn)給定的參數(shù)先驗(yàn)分布條件下，通過(guò)STAN軟件得到解釋變量系數(shù)估計(jì)值的統(tǒng)計(jì)性質(zhì)要優(yōu)于常數(shù)項(xiàng)估計(jì)值的統(tǒng)計(jì)性質(zhì)。2)對(duì)于β1，無(wú)論是在尺度變量σ被參數(shù)化還是未被參數(shù)化的情形下，它的精度都隨著樣本容量的增加而提高，而β0的精度并未呈現(xiàn)出隨著樣本容量的增加而提高的特征。此外，尺度參數(shù)的精度也隨著樣本容量的增加而提高。

本實(shí)驗(yàn)的抽樣估計(jì)結(jié)果是在給定參數(shù)先驗(yàn)分布和抽樣算法的前提下得出的，而對(duì)于貝葉斯分析，先驗(yàn)分布的設(shè)定對(duì)抽樣結(jié)果會(huì)產(chǎn)生影響。下面將考察不同的參數(shù)先驗(yàn)分布對(duì)貝葉斯分位數(shù)回歸結(jié)果的影響。

四、不同先驗(yàn)設(shè)定下分位數(shù)回歸估計(jì)量性質(zhì)的比較

對(duì)于分位數(shù)回歸模型，不存在標(biāo)準(zhǔn)的共軛先驗(yàn)分布(Standard Conjugate Prior Distribution)，因此不容易獲得參數(shù)后驗(yàn)分布的解析表達(dá)式，但仍然可以通過(guò)模擬得到參數(shù)的估計(jì)，前提是在貝葉斯推斷中參數(shù)的先驗(yàn)分布應(yīng)該保證參數(shù)的后驗(yàn)分布是合適分布(Proper Distribution)。Yu和Moyeed指出，在貝葉斯分位數(shù)估計(jì)中，只要待估系數(shù)β的先驗(yàn)分布滿(mǎn)足p(β)∝1，那么β的后驗(yàn)分布會(huì)有合適的分布[2]。王新宇、宋學(xué)鋒進(jìn)一步指出，當(dāng)β服從不合適的先驗(yàn)分布時(shí)，即使尺度參數(shù)σ>0也未知，只要其密度函數(shù)f(σ)滿(mǎn)足：

那么，參數(shù)的聯(lián)合后驗(yàn)分布也是合適分布[7]。也就是說(shuō)，在理論上先驗(yàn)分布可以是任意的(如果沒(méi)有任何的先驗(yàn)信息，可以采用均勻分布)，盡管如此得到的參數(shù)聯(lián)合后驗(yàn)分布也是合適分布。然而，已有文獻(xiàn)中并未考察隨著先驗(yàn)分布信息的不同，參數(shù)后驗(yàn)分布的統(tǒng)計(jì)特征有什么不同，這對(duì)于選擇合適的先驗(yàn)信息來(lái)進(jìn)行貝葉斯分位數(shù)回歸估計(jì)具有指導(dǎo)意義。

(一)實(shí)驗(yàn)設(shè)計(jì)與抽樣結(jié)果

下面針對(duì)一般化的形式，考察不同先驗(yàn)設(shè)定對(duì)估計(jì)結(jié)果的影響。數(shù)據(jù)生成過(guò)程為，Yi=β0+β1x+εi，i=1,2,…，n，除了包括截距項(xiàng)還包括一個(gè)解釋變量，其中β0=1，β1=2，x～unif(n=n，min=0，max=10)?？紤]到誤差項(xiàng)的分布影響因變量的分布，由此可能影響不同分位點(diǎn)的回歸系數(shù)，于是討論誤差項(xiàng)服從正態(tài)分布和存在異方差，即β～N(0,1)、ε～(1+x2)N(0,1)，兩種不同的數(shù)據(jù)生成過(guò)程、不同的先驗(yàn)分布對(duì)估計(jì)結(jié)果的影響。

值得注意的是，由于回歸式中不是只有截距項(xiàng)，因此被估參數(shù)在τ分位點(diǎn)下的真值受誤差項(xiàng)的形式影響。在上面的回歸式中，qτ(y)=β0+β1x+qτ(ε)，當(dāng)誤差服從正態(tài)分布N(0,1)時(shí)，qτ(y)=β0+φ-1(τ)+β1x，其中φ-1(τ)是標(biāo)準(zhǔn)正態(tài)分布的第τ分位數(shù)，因此被估系數(shù)的真值分別是：β0=1+φ-1(τ)，β1=2。當(dāng)誤差項(xiàng)為異方差形式時(shí)，qτ(y)=β0+φ-1(τ)+(β1+φ-1(τ))x1，即被估系數(shù)的真值分別是：β0=1+φ-1(τ)，β1=2+φ-1(τ)，其中φ-1(τ)是標(biāo)準(zhǔn)正態(tài)分布的第τ分位數(shù)。

為比較不同先驗(yàn)分布條件下估計(jì)量小樣本性質(zhì)的差異，設(shè)定參數(shù)先驗(yàn)分布為正態(tài)分布N(μ，σ)。下面給出三種不同的參數(shù)先驗(yàn)分布設(shè)定，隨著σ減小，先驗(yàn)信息逐漸增強(qiáng)，具體形式如下：

先驗(yàn)設(shè)定1：β～N(0,100)

先驗(yàn)設(shè)定2：β～N(0,10)

先驗(yàn)設(shè)定3：β～N(0,1)

根據(jù)上述不同的數(shù)據(jù)生成過(guò)程，結(jié)合不同的先驗(yàn)信息，運(yùn)用HMC抽樣算法對(duì)模型參數(shù)進(jìn)行估計(jì)。為了考察特定參數(shù)的模擬值是否為合適后驗(yàn)分布的合理近似，對(duì)某一先驗(yàn)分布假設(shè)下各個(gè)分位點(diǎn)下的參數(shù)抽樣值軌跡和自相關(guān)圖進(jìn)行分析，結(jié)果表明所有情形下的馬爾科夫鏈均收斂、自相關(guān)圖均隨著滯后期的增加而趨于零。圖3給出了先驗(yàn)設(shè)定2下，樣本容量為75時(shí)，參數(shù)β0和β1的抽樣值軌跡和自相關(guān)圖，其他情形下也類(lèi)似。

(二)結(jié)論

限于篇幅，具體的輸出結(jié)果省略，主要的結(jié)論如下：其一，在一定樣本容量一定分位點(diǎn)下，隨著先驗(yàn)信息的增強(qiáng)，β0和β1抽樣分布的標(biāo)準(zhǔn)差越來(lái)越小，而對(duì)估計(jì)均值的偏離誤差沒(méi)有明顯影響，即增加先驗(yàn)信息可以提高參數(shù)估計(jì)量的精度。其二，0.5分位點(diǎn)下，β0和β1抽樣分布的標(biāo)準(zhǔn)差均小于0.25和0.75分位點(diǎn)下的標(biāo)準(zhǔn)差，因此在采用Gibbs抽樣算法進(jìn)行貝葉斯分位數(shù)回歸分析，尤其是預(yù)測(cè)時(shí)，在中位數(shù)下系數(shù)的預(yù)測(cè)精度更高。其三，在誤差項(xiàng)服從正態(tài)分布的情形，一定樣本容量一定分位點(diǎn)下，參數(shù)后驗(yàn)分布的標(biāo)準(zhǔn)差小于服從非對(duì)稱(chēng)拉普拉斯分布情形和誤差項(xiàng)存在異方差下對(duì)應(yīng)的標(biāo)準(zhǔn)差，這表明在貝葉斯分位數(shù)回歸分析中，要求誤差服從正態(tài)分布可以得到系數(shù)優(yōu)良的統(tǒng)計(jì)性質(zhì)。其四，在三種不同的先驗(yàn)分布下，β0和β1抽樣分布的標(biāo)準(zhǔn)差均在0～2之間，尤其是在先驗(yàn)設(shè)定1下，先驗(yàn)的標(biāo)準(zhǔn)差為100，幾乎是無(wú)信息先驗(yàn)的情形下，得到參數(shù)后驗(yàn)分布的標(biāo)準(zhǔn)差也不超過(guò)2，這也驗(yàn)證了已有的結(jié)論：貝葉斯分位數(shù)回歸分析中即使先驗(yàn)分布是非合適的，得到的后驗(yàn)分布卻是合適的[2，7]。

此外，還可以發(fā)現(xiàn)：無(wú)論設(shè)定哪種先驗(yàn)分布，β1的標(biāo)準(zhǔn)差總是小于同一分位點(diǎn)下β0的標(biāo)準(zhǔn)差，表明通過(guò)STAN軟件得到解釋變量系數(shù)估計(jì)值的統(tǒng)計(jì)性質(zhì)要優(yōu)于常數(shù)項(xiàng)估計(jì)值的統(tǒng)計(jì)性質(zhì)，而且在各種先驗(yàn)分布設(shè)定下，隨著樣本容量增加，β0和β1的精度也隨之提高。

圖3 β0和β1的抽樣值軌跡和自相關(guān)圖

綜上，使用STAN軟件采用HMC抽樣算法進(jìn)行貝葉斯分位數(shù)回歸分析時(shí)，針對(duì)先驗(yàn)分布的選擇問(wèn)題，可以使用無(wú)信息先驗(yàn)，如果已經(jīng)知道分位數(shù)回歸系數(shù)為正態(tài)的情況下，就可以使用正態(tài)先驗(yàn)分布，而且可以考慮使用標(biāo)準(zhǔn)正態(tài)分布以提高參數(shù)估計(jì)量的精度；而先驗(yàn)分布的選擇對(duì)于參數(shù)估計(jì)量的偏誤沒(méi)有顯著影響。

五、結(jié)束語(yǔ)

本文在STAN中編譯非對(duì)稱(chēng)拉普拉斯分布的似然函數(shù)，并在此基礎(chǔ)上實(shí)現(xiàn)了貝葉斯分位數(shù)回歸參數(shù)在STAN軟件中的模擬抽樣。同時(shí)，將尺度參數(shù)進(jìn)行參數(shù)化，通過(guò)HMC抽樣算法得到其后驗(yàn)分布。實(shí)驗(yàn)結(jié)果表明，與未參數(shù)化時(shí)相比，參數(shù)化后得到的估計(jì)量統(tǒng)計(jì)性質(zhì)更好。在貝葉斯分析中，設(shè)定不同的先驗(yàn)分布意味著不同的主觀信念，這種不同的信念最終會(huì)體現(xiàn)在參數(shù)的后驗(yàn)密度中，從而對(duì)估計(jì)結(jié)果產(chǎn)生影響。因此，本文針對(duì)連續(xù)因變量的貝葉斯分位數(shù)回歸估計(jì)方法，模擬不同先驗(yàn)分布設(shè)定對(duì)參數(shù)估計(jì)量統(tǒng)計(jì)性質(zhì)的影響。實(shí)驗(yàn)結(jié)果表明，合適的先驗(yàn)分布可以提高HMC抽樣估計(jì)量的統(tǒng)計(jì)性質(zhì)。

分位數(shù)回歸和貝葉斯分析方法作為計(jì)量經(jīng)濟(jì)學(xué)的兩大前沿分支，通過(guò)設(shè)定模型誤差項(xiàng)服從非對(duì)稱(chēng)拉普拉斯分布而巧妙地結(jié)合起來(lái)，實(shí)現(xiàn)了對(duì)分位數(shù)回歸模型的貝葉斯參數(shù)估計(jì)。本文僅探究了貝葉斯方法在分位數(shù)回歸中的某些優(yōu)越性，隨著研究的深入，未來(lái)對(duì)貝葉斯分位數(shù)估計(jì)方法的探索將越來(lái)越廣泛。

感謝波士頓學(xué)院肖志杰教授對(duì)本文的有益指導(dǎo)；感謝上海社會(huì)科學(xué)院研究生陳柯和王浩宇在軟件方面的協(xié)助。文責(zé)自負(fù)。

[1] Fatti L P，Senaoana E M.Bayesian Updating in Reference Centile Charts[J].Journal of the Royal Statistical Society，1998，161(1).

[2] Yu K，Moyeed R A.Bayesian Quantile Regression[J].Statistics and Probability Letters，2001，54(4).

[3] Koenker R，Bassett G S.Regression Quantiles[J].Econometrica，1978，46(1).

[4] Yu Keming，Philippe van Kerm，Jin Zhang.Bayesian Quantile Regression：An Application to the Wage Distribution in 1990s Britain[J].The Indian Journal of Statistics，2005(5).

[5] Alhamzawi R，Yu K.Power Prior Elicitation in Bayesian Quantile Regression[J].Journal of Probability and Statistics，2011(2).

[6] Benoit D F，Van den Poel D.Binary Quantile Regression：A Bayesian Approach Based on the Asymmetric Laplace Distribution[J].Journal of Applied Econometrics，2012，27(7).

[7] 王新宇，宋學(xué)鋒．基于貝葉斯分位數(shù)回歸的市場(chǎng)風(fēng)險(xiǎn)測(cè)度模型與應(yīng)用[J].系統(tǒng)管理學(xué)報(bào)，2009，18(1).

[8] 白仲林，隋雯霞，劉傳文.混合貝塔分布隨機(jī)波動(dòng)模型及其貝葉斯分析[J].統(tǒng)計(jì)與信息論壇，2013，32(4).

[9] Linzer D.Dynamic Bayesian Forecasting of Presidential Elections in the States[J].Journal of the American Statistical Association，2013，108(501).