一種基于ALD的貝葉斯分位數回歸估計

韋學永,金良瓊,沈婷,朱偉業(yè)成

(貴州民族大學 數據科學與信息工程學院,貴州 貴陽,550025)

自最小二乘回歸估計被提出后,其應用就遍及社會經濟的各個領域,但它只考慮了條件均值下因變量與自變量之間的關系,而在實際的應用中,有時需要分析數據的各個分位點下因變量與自變量之間的關系,因此1978年Koenker等第一次提出了“分位數回歸”的概念,補充了基于最小二乘估計對只關注條件均值估計的不足[1]。隨著現代計量經濟學研究的發(fā)展,分位數回歸估計引起了研究學者們的廣泛關注,在傳統(tǒng)的分位數回歸進行估計的方法中,基于貝葉斯分析的方法比頻率學派的方法更有優(yōu)勢,因此越來越多的學者更加關注基于貝葉斯的分位數回歸估計的研究。Koenker等研究了分位數與非對稱拉普拉斯分布(Asymmetric Laplace Distribution,ALD)之間的關系[2]。2001年,Yu等[3]基于文獻[1]提出的ALD對貝葉斯分位數回歸的參數進行估計,同時證明了即使先驗分布是不適當的,得出的后驗分布也是合適的,故可選無信息的先驗分布;2009年,王新宇等指出ALD的尺度參數不應該假定為常數1[4]。何鳳霞等對分位數回歸在R軟件中的應用領域方面做出了相關總結[5]。2015年肖北芳通過MCMC算法的M-H抽樣模擬得到參數的后驗分布,驗證了金融發(fā)展與城鄉(xiāng)差距的倒“U”關系在我國樣本縣域經濟中存在統(tǒng)計顯著[6]?；舸鋫ピ?018年基于貝葉斯分位數回歸對因變量為連續(xù)和離散的情況下作了研究[7]。2019年,方麗婷等提出了空間滯后分位數回歸模型的貝葉斯估計方法,該方法在小樣本條件下具有良好的估計效果和穩(wěn)健性[8]。同年,邸俊鵬等對基于貝葉斯估計方法的二元選擇分位數回歸模型進行研究,模擬實驗結果表明,小樣本下對二元選擇分位數回歸模型進行貝葉斯估計比頻率學派方法的估計效果更佳,統(tǒng)計推斷更準確[9]。Zhou等針對未知跳躍數的情況提出了一種使用局部多項式技術檢測跳躍的新程序,通過模擬實驗證明提出的方法在廣泛的實際環(huán)境中表現良好[10];Xu等提出了一種貝葉斯非參數方法來同時估計非交叉的非線性分位數曲線,當數據稀疏時,該模型可以更好地恢復響應分布的分位數[11];張發(fā)趕等利用非對稱拉普拉斯分布提出一種新的混合分位數回歸模型,并基于所提出模型及算法對城市房價數據進行了分析[12];潘瑩麗等在缺失數據的條件下對模型進行分位數回歸[13];Liu等基于貝葉斯方法對離散響應的回歸模型進行研究,同時證實了該方法得到的后驗分布不僅與響應的真實分布無關,而且對于未知模型參數的不正確先驗分布也是正確的[14]。

WinBUGS軟件在1989年由英國劍橋的生物研究所編譯,是一種以MCMC方法為基礎的貝葉斯分析方法,將所有的未知參數看作隨機變量,隨之被相關研究學者引用到分位數回歸的領域中,2018年,邸俊鵬等基于WinBUGS軟件用泊松分布指定ALD的似然函數,并且驗證了尺度參數應該被參數化,同時指出可用二項分布或負二項分布來指定ALD的似然函數[15]?；诖?本文采用二項分布指定ALD的似然函數,在WinBUGS軟件中實現貝葉斯分位數回歸估計的相關分析,當貝葉斯參數估計的后驗分布的表達式復雜或難以表示時,采用Gibbs抽樣算法可以得到后驗分布相關的統(tǒng)計性質。

1 Gibbs抽樣

在統(tǒng)計分析中,常用的計算方法有極大似然估計算法和貝葉斯算法。在現實應用中,由于后驗分布一般都是非標準形式、復雜、高維的分布,因此采用貝葉斯算法中的馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)方法。在MCMC方法中,常用的2種算法是:Gibbs抽樣和Metropolis-Hastings方法[17],Gibbs抽樣是MCMC方法里最簡單且應用最廣泛的抽樣方法,本文主要是基于Gibbs抽樣對貝葉斯分位數回歸的參數進行估計。

Gibbs抽樣的具體步驟如下:

當T只含有一個元素時,稱為單元素Gibbs抽樣(single-site Gibbs sampler),單元素Gibbs抽樣的具體步驟如下[16]:

第三步:令t=t+1,返回第二步,依此重復進行。

2 非對稱拉普拉斯分布的貝葉斯參數估計

在基于ALD的貝葉斯分位數回歸估計中,假定線性的隨機誤差項服從ALD,并且用二項分布指定ALD在WinBUGS軟件中實現,得到參數后驗分布的相關統(tǒng)計性質。

2.1 非對稱拉普拉斯分布

2.2 線性模型

假定分位數回歸線性模型為

令Z(x;β)=β0+β1x,則y=Z(x;β) +u,其中:因變量為y;自變量為x;待估參數為β0,β1;隨機誤差項u～ALD(μ,σ,q),即y～ALD(Z(x;β),σ,q)。

2.3 二項分布

記X為n重伯努利實驗中成功(記為事件A)的次數,則X的可能取值為0,1,…,n,記p為每次實驗中A發(fā)生的概率,即重伯努利實驗的基本結果記作ω=(ω1,ω2,…,ωn)。若某個樣本點則ω1,ω2,… ,ωn中有k個A,n-k個,由獨立性可得,P(ω)=pk(1-p)n-k,而事件{X=k}中共有個這樣的ω,故X的分布列為。1,…,n這個分布稱作二項分布,記為X～b(n,p)[4]

2.4 二項分布指定的非對稱拉普拉斯分布

在WinBUGS軟件中,只包含了23種經常用的分布來指定先驗分布的似然函數,然而對于泊松分布、二項分布和負二項分布等,則需要相關的一些編程技巧和數學的轉換才能在WinBUGS軟件中實現參數的估計。首先在WinBUGS軟件中編譯二項分布指定ALD的似然函數的程序;然后在Rstudio軟件中用“bugs()”函數調用WinBUGS編譯且已保存的程序,并基于WinBUGS軟件進行Gibbs抽樣,實現了貝葉斯分位數回歸的模擬。

3 數值模擬

生成隨機數據的模型為Y=1+ 2X+ε,其中X～U( 0,10),ε～ALD( 0 ,1,q)。利用生成的隨機數據對模型(1)進行模擬。

3.1 尺度參數σ設定為1

當尺度參數σ設定為1,待估參數為β0和β1時,假定待估參數β0和β1的先驗分布服從正態(tài)分布,由貝葉斯定理可知,參數的聯(lián)合后驗分布為其中,f()β為待估計系數β的先驗分布,則數據的對數密度函數可表示為

將式(4)代入式(2),可得參數pi的表達式

將式(5)代入式(3),故參數的聯(lián)合后驗分布可表示為

用二項分布指定的ALD,當迭代10 000次,預燒期為5 000次時對參數的聯(lián)合后驗分布進行Gibbs抽樣得到的軌跡圖和自相關圖更加平穩(wěn);在相同的迭代次數和預燒期不同的先驗分布、分位點及樣本容量下參數β0和β1的5 000個抽樣值的軌跡圖均在設定值上下波動,自相關系數隨著滯后期的增加而逐漸趨于零,故馬爾科夫鏈收斂。限于篇幅,僅展示當參數先驗分布設定為β～N(0,100),在0.75分位點下且樣本容量n為75時,參數β0和β1的5 000個抽樣值的軌跡圖、密度圖以及自相關圖,結果顯示參數β0和β1的5 000個抽樣值的軌跡圖均在設定值上下波動,密度圖的寬帶足夠小,且自相關系數隨著滯后期的增加而逐漸趨于零,如圖1所示。

圖1(a) β0抽樣值的軌跡圖、密度圖及自相關圖

圖1(b) β1抽樣值的軌跡圖、密度圖及自相關圖

如表1所示,對分位數回歸的線性模型參數進行Gibbs抽樣,迭代次數為10 000次,預燒期為5 000次,可得到參數后驗分布的均值、標準差和模擬誤差等。

表1 σ設定為1時參數β0和β1的均值、標準差及MC誤差

續(xù)表1

由表1可知:(1)當先驗分布和分位點相同時,隨著樣本容量n的增大,待估參數β0和β1的標準差和MC誤差逐漸變小;(2)當先驗信息、樣本容量n以及分位點相同時,β1的標準差和MC誤差總是小于β0的標準差和MC誤差,即在WinBUGS軟件中得到的系數項的估計精度往往高于常數項的估計精度;(3)在樣本容量和分位點相同條件下,隨著先驗分布的增強,參數β0和β1的標準差和MC誤差變小,即增強先驗信息時可提高參數的估計精度;(4)在這3種先驗設定中,當樣本容量和先驗信息相同時,除了先驗信息為β～N(0,100),樣本量為100,參數β0和β1在0.25分位點下的標準差和MC誤差比在0.5分位點、0.75分位點下的標準差和MC誤差更小,其他條件下,參數β0和β1在0.25分位點下的標準差和MC誤差均比在0.5分位點、0.75分位點下的標準差和MC誤差更大,即參數β0和β1在中位數和高分位數的估計精度較高。

3.2 尺度參數σ參數化

ALD中的尺度參數設定為1時,盡管參數的估計精度較高,但在實際應用中是不恰當的,因此,應當對ALD中的尺度參數進行參數化[4]。

當參數σ,β0和β1均待估計時,設定σ的初始值為1,假設待估參數β0和β1的先驗分布為正態(tài)分布和尺度參數σ的先驗分布為卡方分布,當尺度參數σ的卡方分布自由度為4時,3個待估參數的平穩(wěn)度更高,由貝葉斯定理可知,參數的聯(lián)合后驗分布為其中,f(β)為待估計系數β的先驗分布,g(σ)為尺度參數σ的先驗分布,且σ～χ2(4),則數據的對數密度函數為

將式(7)代入式(2),可得參數pi的表達式

將式(8)代入式(3),可得參數的聯(lián)合后驗分布的表達式

圖2為二項分布指定ALD的似然函數,當迭代10 000次,預燒期為5 000次時對參數的聯(lián)合后驗分布進行Gibbs抽樣得到的軌跡圖和自相關圖更加平穩(wěn);在相同的迭代次數和預燒期不同的先驗分布、分位點及樣本容量下參數β0、β1和σ的5 000個抽樣值的軌跡圖均在設定值上下波動,自相關系數隨著滯后期的增加而逐漸趨于零,因此抽樣構成的馬爾科夫鏈均收斂。限于篇幅,僅展示當參數先驗分布設定為β～N(0,100),在0.75分位點下且樣本容量n為75時,參數β0、β1和σ的5 000個抽樣值的軌跡圖、密度圖以及自相關圖,結果顯示參數β0、β1和σ的5 000個抽樣值的軌跡圖均在設定值上下波動,密度圖的寬帶足夠小,且自相關系數隨著滯后期的增加而逐漸趨于零。

圖2(a) β0抽樣值的軌跡圖、密度圖及自相關圖

圖2(b) β1抽樣值的軌跡圖、密度圖及自相關圖

圖2(c) σ抽樣值的軌跡圖、密度圖及自相關圖

如表2所示,對分位數回歸的線性模型參數進行Gibbs抽樣,迭代次數為10 000次,預燒期為5 000次,可得到參數后驗分布的均值、標準差和模擬誤差等,由表2可知:(1)當先驗分布和分位點相同時,隨著樣本容量n的增大,參數β0、β1和σ的標準差和MC誤差逐漸變小;(2)在先驗信息、樣本量以及分位點相同條件下,標準差的大小為β1＜σ＜β0,而MC誤差的大小為σ＜β1＜β0;(3)在樣本容量和分位點相同條件下,隨著先驗分布的增強,參數β0、β1和σ的標準差和MC誤差變小,即增強先驗信息時,可提高參數的估計精度;(4)當先驗信息相同且樣本量為100時,參數β0、β1和σ在0.75分位點下的標準差和MC誤差比在0.25分位點、0.5分位點下的標準差和MC誤差更小。

表2 參數σ,β0和β1均待估時的均值、標準差及MC誤差

續(xù)表2

對比表1和表2可知,當先驗信息、樣本量、參數與分位點相同時,表2中參數β0和β1的標準差與MC誤差均小于表1中參數β0和β1的標準差與MC誤差。因此,對尺度參數σ進行參數化可以提高參數β0和β1的估計精度,即尺度參數σ應參數化,而不是設定為1。

用二項分布指定ALD的實驗結果與邸俊鵬用泊松分布指定ALD的似然函數[15]的實驗結果相比,尺度參數進行參數化后,2種實驗的參數統(tǒng)計量性質效果一樣好。

綜上所述,當尺度參數未參數化時,利用二項分布和用泊松分布指定ALD的似然函數得到的參數統(tǒng)計量性質效果一樣;盡管先驗分布不適當,但是得到的后驗分布也是適當的。

4 結束語

本文利用二項分布指定ALD的似然函數,通過對分位數回歸線性模型進行Gibbs抽樣得到參數后驗分布的相關統(tǒng)計性質,通過模擬實驗證實了尺度參數被參數化得到參數后驗分布的估計量統(tǒng)計性質比尺度參數假定為1時的估計效果更好且更穩(wěn)健;同時也證實了適當的先驗分布可以提高待估參數的估計精度。接下來可考慮用負二項分布指定ALD的似然函數,分析所得結果是否與此方法一樣理想或比此方法更好。