基于M－H抽樣算法的貝葉斯Probit分位回歸模型研究＊

朱慧明，李 榮，曾昭法，虞克明

（1.湖南大學(xué)工商管理學(xué)院，湖南長沙 410082； 2.湖南大學(xué)金融與統(tǒng)計學(xué)院，湖南長沙 410079；3.Brunel大學(xué)數(shù)學(xué)系，倫敦 UB8 3PH）

Probit模型廣泛應(yīng)用于經(jīng)濟布局、企業(yè)選址、交通問題、就業(yè)問題和購買決策等經(jīng)濟決策領(lǐng)域.在經(jīng)典計量經(jīng)濟學(xué)模型中，因變量通常為連續(xù)變量，但在經(jīng)濟分析，如銷售或購買某種商品、犯罪、遷移、生育和患病與否等，因變量只取兩個值，可以用Probit模型（0或1）來度量經(jīng)常面臨選擇問題，從而在可供選擇的方案中做出決策.離散Probit選擇模型在經(jīng)濟研究中大量應(yīng)用，許多學(xué)者用Probit模型回歸方法對實際問題進行分析.林相森，艾春榮（2008）［1］對中國居民醫(yī)療需求影響因素用有序的Probit回歸模型進行分析；陳磊等（2009）［2］運用序次Probit回歸模型對離散股票價格進行建模分析，結(jié)果表明該模型可反映價格離散特征；張欣［3］（2010）用Probit離散選擇模型對船舶溢油事故的選擇概率進行預(yù)測.Koenker和Bassett［4］于1978年提出了分位回歸理論，不少學(xué)者把Probit均值回歸模型推廣到Probit分位回歸模型，Delgado（2001）［5］從理論上證明了Probit分位回歸的重抽樣方法對參數(shù)最大得分估計方法的可行性，但此方法計算復(fù)雜，并且只能適用于低維小樣本問題.Korads（2006）［6］把核密度方法推廣到一般的Probit分位數(shù)回歸，但核密度估計方法對擾動項分布的光滑性要求很強，同時表明當(dāng)仿真樣本數(shù)等于1 000時，很難去估計參數(shù)的標(biāo)準(zhǔn)差.Koenker（2005）［7］指出，選擇光滑參數(shù)對最大得分函數(shù)的參數(shù)估計會產(chǎn)生較大的影響.Florios和Skouras（2008）［8］指出針對Probit分位回歸的最大得分函數(shù)的優(yōu)化問題，提出了混合整數(shù)規(guī)劃來解決此問題，但仍然不能得到待估參數(shù)的統(tǒng)計性質(zhì).總之，對Probit模型進行建模分析，存在參數(shù)隨機化條件下建模困難和參數(shù)求解的問題.

本文利用貝葉斯統(tǒng)計推斷理論和分位回歸技術(shù)對Probit模型進行分析，解決參數(shù)隨機化條件下的Probit模型構(gòu)建問題，對離散Probit模型分別用貝葉斯Probit分位回歸模型與分位回歸模型和光滑分位回歸模型對參數(shù)估計比較分析，得到不同分位點模型參數(shù)的后驗分布.

1 Probit分位回歸模型

Probit離散選擇回歸模型在社會生活、經(jīng)濟、統(tǒng)計決策和數(shù)據(jù)挖掘中有著廣泛的應(yīng)用，Probit離散選擇模型其具體形式如下：

其中此處yi為第i個樣本的指示變量，取值0或1，為潛變量，＝（x1，x2，…，xk），＝（β1，β2，…，βk），εi為隨機擾動項，i＝1，2，…，n.

下面考慮選擇yi＝1的概率.設(shè)F（ε｜x＝xi）表示給定x＝xi下隨機擾動項的條件分布函數(shù)，從式（1）可知

從普通的Probit均值回歸模型推廣到一般的Probit分位回歸模型過程中，待估參數(shù)隨著分位水平的變化而改變，從而Probit分位回歸模型可以更加詳細(xì)地描述解釋變量x對潛變量Y＊的條件分布的位置和形狀的影響.

假設(shè)隨機變量Y＊與解釋變量向量x之間具有線性函數(shù)關(guān)系

其中x′＝（x1，x2，…，xk），β＝（β1，β2，…，βk）′，ε為概率分布函數(shù)F的隨機擾動項.顯然，對于給定的x′，隨機變量Y＊的p分位數(shù)等于x′β＋F－1（p），即

不難看出，其等價形式為：

此處β（p）＝（β1（p），β2（p），…，βk（p））′，也就是說，β（p）與分位數(shù)p的取值大小有關(guān).顯然，對于不同的p值，可以得到相應(yīng)的多個分位點的回歸模型.給定一組觀察值β（τ）的估計為：

此處ρp（u）＝u（p－I（u＜0））為損失函數(shù)，I（u＜0）為示性函數(shù)；當(dāng)u＜0時，I（ui＜0）＝1；否則，I（u＜0）＝0.由于β（p））等價于

在式（1）的條件下，可求得Probit分位回歸參數(shù)估計形式如下：此處sgn（·）為符號函數(shù).

2 Probit分位回歸模型的貝葉斯分析

非對稱Laplace分布是應(yīng)用貝葉斯方法分析Probit離散選擇回歸模型的重要工具.如隨機變量y服從偏度參數(shù)p，尺度參數(shù)σ和位置參數(shù)μ的非對稱Laplace分布［9］，其密度為：

記作Y～ALD（μ，σ，p）；特別地，當(dāng)μ＝0，σ＝1時，記作Y～ALD（p）.

如ε～ALD（μ，σ，p），可得如下似然函數(shù)

根據(jù)貝葉斯定理，可知待估參數(shù)β的后驗分布

此處p（β）為β的先驗分布.

本文只考慮ALD（μ，σ，p）分布σ＝1時的情形，給定一組數(shù)據(jù)y＝（y1，y2，…，yn）和分位數(shù)p值，可得Probit分位回歸中待估參數(shù)β和潛變量Y＊的聯(lián)合后驗密度為：

此處p（β）為β的先驗分布，I（.）為指示變量.

由于式（10）聯(lián)合后驗密度的分布的形式未知，所以直接抽樣不可能完成，可以通過數(shù)據(jù)擴充來實現(xiàn)對β的抽樣，可以把式（10）分解兩個全條件后驗分布可以完成.從式（10），可得Y＊的全條件后驗密度為

當(dāng)yi＝1時，y＊i左截斷0；yi＝0時，y＊i右截斷0.

由于直接對非對稱Laplace分布抽樣困難，利用非對稱Laplace分布的性質(zhì)，可以分解成兩個常用的分布來完成上述抽樣.文獻［10］在對回歸模型進行貝葉斯分位分析時，研究證明選擇參數(shù)的先驗信息為無信息先驗分布，可以得到待估參數(shù)合適的后驗分布.從式（10），可得β的后驗密度為：

針對貝葉斯分位Probit回歸模型參數(shù)后驗分布的解析形式不確定的問題，采MCMC模擬解決貝葉斯估計過程中遇到的高維積分問題.利用M－H抽樣模擬Markov鏈，此鏈的平穩(wěn)分布可以近似看作模型參數(shù)的后驗條件分布，用MCMC模擬得到的樣本來估計模型的參數(shù)［11］.

為了使Markov鏈的平穩(wěn)分布就是目標(biāo)分布，接受概率通常選擇為：

選擇對稱建議分布（proposal distribution）為q（β′，β）＝q（β，β′），此時α（β，β′）＝min｛1，π（β′｜y）／π（β｜y）｝.如果Markov在時刻t處于狀態(tài)β，即β（t）＝β，則由q（.｜β）產(chǎn)生一個轉(zhuǎn)移β→β＇；然后根據(jù)α（β，β′）決定是否轉(zhuǎn)移.

假設(shè)給定初始值Θ（0）＝（y＊（0），β（0）（p）），并設(shè)（y＊（j－1），β（j－1）（p））為（y＊，β（p））的第j－1次抽樣結(jié)果.則第j次抽樣的過程包括以下兩個步驟：

由上完成了一次M－H迭代過程，即完成了由（y＊（j－1）（p），β（j－1）（p））到（y＊（j）（p），β（j）（p））的轉(zhuǎn)移；不斷重復(fù)以上步驟直到鏈條達到穩(wěn)定狀態(tài)，經(jīng)過t次迭代，可以得到（y＊（t）（p），β（t）（p））.

如果總共迭代N次，則獲得如下迭代結(jié)果：

一般說來，在抽樣過程的初始階段，過程處于非平穩(wěn)狀態(tài)；因此，在估計模型參數(shù)時，去掉最初抽樣的M個點｛（y＊（j），β（j）（p）），j＝1，2，…，M｝，利用隨后N－M次迭代所獲得的隨機數(shù)，形成一個容量大小為N－M的樣本｛（y＊（j），β（j）（p）），j＝M＋1，M＋2，…，N｝，據(jù)此求出β（p）的估計及其方差，即：

3 MCMC仿真分析

利用隨機模擬生成服從Probit回歸的200個數(shù)據(jù)，首先用這些數(shù)據(jù)來驗證貝葉斯方法應(yīng)用于Probit分位回歸模型的有效性，然后把Probit分位回歸（PQR）、光滑Probit分位回歸（SPQR）［6］和貝葉斯Probit分位回歸（BPQR）在0.5分位點進行參數(shù)估計結(jié)果的比較.實驗數(shù)據(jù)的生成過程如下.

對于異方差模型：

此處

選擇先驗分布為均勻分布，M－H抽樣算法中的建議分布為指數(shù)分布，由M－H抽樣可以模擬得到各個參數(shù)的邊緣后驗分布，把模型（13）用M－H抽樣10 000次估計待估參數(shù)，得到β1和β2在0.05，0.10，0.25，0.75，0.90，0.95等6個分位點參數(shù)β1的動態(tài)迭代軌跡圖（圖1）、參數(shù)β2的動態(tài)迭代軌跡圖（圖2），從參數(shù)的迭代軌跡圖可以發(fā)現(xiàn)馬爾可夫鏈?zhǔn)諗浚f明MCMC仿真過程是平穩(wěn)的.

圖1 不同分位數(shù)下參數(shù)β1動態(tài)迭代軌跡Fig.1 Dynamic iterative Chain track of theβ1 parameter in different quantiles

圖3和圖4分別給出不同分位數(shù)下參數(shù)后驗分布的核密度估計曲線，從回歸的后驗分布密度圖可以看出，其表現(xiàn)為比較平滑且呈鐘型，說明M－H算法有效地模擬了模型中各參數(shù)的邊緣后驗分布.在低分位點，X對y的影響為負(fù)，在高分位點，X對y的影響為正，說明雖然y的取值只是兩元，誤差項不是標(biāo)準(zhǔn)的Laplace分布，但通過貝葉斯probit分位回歸，可以得到在每個分位點的影響情況.

圖2 不同分位數(shù)下參數(shù)β2動態(tài)軌跡Fig.2 Dynamic iterative Chain track of theβ2 parameter in different quantiles

圖3 不同分位數(shù)下參數(shù)β1的后驗密度圖（樣本量10 000）Fig.3 Posterior density estimation of theβ1 parameter in different quantiles

圖4 不同分位數(shù)下參數(shù)β2的后驗密度圖（樣本量10 000）Fig.4 Posterior dersity estimation of theβ2 parameter in different quantiles

在樣本量分別為200，500和800的情形下，分別用PQR，SPQR和BPQR對模型（13）的參數(shù)進行估計，在0.5分位點，得到參數(shù)估計的偏差、單位均方誤差（RMSE）和95%的置信區(qū)間如表1所示.

表1 參數(shù)估計的偏差、單位均方誤差和95%的置信區(qū)間Tab.1 Biased，RMSE and the 95%confidance level of parameters estimation

4 結(jié) 論

本文針對Probit分位回歸在參數(shù)隨機化條件下的建模問題，利用非對稱Laplace分布實現(xiàn)對Probit分位回歸模型的貝葉斯推斷.假定模型中的待估參數(shù)先驗分布為無信息先驗，利用M－H抽樣算法得到待估參數(shù)的后驗邊緣分布，模擬在不同分位水平下參數(shù)的MCMC迭代動態(tài)軌跡圖，參數(shù)的貝葉斯估計和邊緣后驗分布，各分位水平下參數(shù)的MCMC迭代軌跡是收斂的，說明M－H抽樣很好地模擬了參數(shù)的后驗邊緣分布，各分位水平下參數(shù)的標(biāo)準(zhǔn)差比較小，且參數(shù)的后驗密度呈鐘型，貝葉斯Probit分位回歸模型解決了因變量為離散變量時參數(shù)估計不確定問題.與傳統(tǒng)的Probit分位回歸方法相比，貝葉斯分位回歸方法可以合理地解釋Probit模型中參數(shù)隨分位數(shù)變化的特點，得到更加準(zhǔn)確有效的參數(shù)估計.

［1］ 林相森，艾春榮.我國居民醫(yī)療需求影響因素的實證分析－有序probit模型的半?yún)?shù)估計［J］.統(tǒng)計研究，2008，25（11）：40－45.LIN Xiang－sen，AI Chun－rong.Determinants of Chinese resident’s demand for medical care－an application of semiparametric estimation of ordered probit model［J］.Statistical Research，2008，25（11）：40－45.（In Chinese）

［2］ 陳磊，李平，曾勇.基于序次probit模型的離散股價研究［J］.系統(tǒng)工程學(xué)報，2009，24（5）：561－567.CHEN Lei，LI Ping，ZENG Yong.Research on discrete stock price based on ordered probit model［J］.Journal of Systems Engineering，2009，24（5）：561－567.（In Chinese）

［3］ 張欣.船舶溢油事故報告選擇概率預(yù)測研究［J］.數(shù)理統(tǒng)計與管理，2010，29（2）：197－204.ZHANG Xin.Forecast of accident report choice probability in ship oil－spill emergency［J］.Journal of Applied Statistics and Management，2010，29（2）：197－204.（In Chinese）

［4］ KOENKER R，BASSETT G.Regression quantiles［J］.Econometrica，1978，46（1）：33－50.

［5］ DELGADO M A，RODRIGUEZ－POO J M，WOLF M.Subsampling inference in cube root asymptotics with an application to Manski’s maximum score estimator［J］.Economic Letters，2001，73（2）：241－250.

［6］ KORDAS G.Smoothed binary regression quantiles［J］.Journal of Applied Econometrics，2006，21（3）：387－407.

［7］ KOENKER R.Quantile regression［M］.New York：Cambridge University Press，2005.

［8］ FLORIOS K，SKOURAS S.Exact computation of max weighted score estimators［J］.Journal of Econometrics，2008，146（1）：86－91.

［9］ YU K，MOYEED R A.Bayesian quantile regression［J］.Statistics and Probability Letters，2001，54（4）：437－447.

［10］YU K，STANDER J.Bayesian analysis of a tobit quantile regression model［J］.Journal of Econometrics，2007，137（1）：260－276.

［11］曾慧芳，朱慧明，李素芳，等.基于MH算法的貝葉斯分位自回歸模型［J］.湖南大學(xué)學(xué)報：自然科學(xué)版，2010，37（2）：88－92.ZENG Hui－fang，ZHU Hui－ming，LI Su－fang，et al.Bayesian inference on the quantile autoregressive models with metropolis－h(huán)astings algorithm［J］.Journal of Hunan University：Nature Sciences，2010，37（2）：88－92.（In Chinese）