半?yún)?shù)順序回歸的貝葉斯推斷

趙煥麗, 何幼樺

(上海大學(xué)理學(xué)院, 上海 200444)

Cox[1]提出的logit 模型是二分類問題最常用的解決方案.而對于多分類或多等級問題,McCullagh[2]擴展二分類logit 模型為比例優(yōu)勢模型.在比例優(yōu)勢模型中, 設(shè)響應(yīng)變量z為具有有限個等級1<2< ··· < K的順序變量, 在自變量x ∈Rk×1處, 記pj(x) =P(z=j|x),j=1,2,··· ,K, 其累積概率

由此, 在x處,z所屬等級落入{1,2,··· ,j}或{j+1,j+2,··· ,K}的概率分別為γj(x)和1?γj(x).將{1,2,··· ,j}和{j+1,j+2,··· ,K}作為兩個分類, 由logit 模型[1]得到針對順序變量的比例優(yōu)勢模型

式中:αj關(guān)于等級j單調(diào)遞增;β ∈Rk×1為待估計參數(shù).式(2)是一個參數(shù)回歸模型.

當(dāng)響應(yīng)變量與解釋變量之間的關(guān)系不能用有限個參數(shù)描述時, 一般引入非參數(shù)模型.當(dāng)響應(yīng)變量與部分解釋變量滿足線性關(guān)系, 而與其他解釋變量的關(guān)系不能用有限個參數(shù)描述時, 則建立半?yún)?shù)回歸模型是合適的.熊笛等[3]用一個連續(xù)非線性函數(shù)替換式(2)中的線性部分, 建立了一種半?yún)?shù)順序回歸模型

并構(gòu)造了參數(shù)部分αj的最小二乘估計和非參數(shù)部分f(·)的局部線性估計量.

1986 年, Engle等[4]提出了一般形式的半?yún)?shù)模型

式中: (x,u)為解釋變量;β和f(·)為待估計參數(shù)和待估計函數(shù).

考慮到半?yún)?shù)模型適用范圍廣泛, 本工作在模型(4)的基礎(chǔ)上, 建立了更一般形式的半?yún)?shù)順序回歸模型

式中:γj(x,u) =P(z≤j|x,u)為給定解釋變量x,u時z所屬等級不超過j的概率.模型(5)是模型(2)與模型(3)的推廣.

對于半?yún)?shù)模型的估計, 研究人員給出了多種經(jīng)典估計方法[5-7].關(guān)于貝葉斯方法,Koop等[8]對局部線性正態(tài)回歸模型的半?yún)?shù)推斷發(fā)展了一種新的貝葉斯方法, 但該方法只討論了取共軛先驗下的正態(tài)線性模型, Koop等[9]使用了類似的方法.李琪琪[10]導(dǎo)出了半?yún)?shù)回歸模型中參數(shù)的貝葉斯最小風(fēng)險線性無偏估計, 并討論了相對于最小二乘加權(quán)估計(least squares weighted estimation, LSWE)的優(yōu)良性.Dimitrakopoulos[11]在誤差項服從Dirichlet 過程時發(fā)展了一種估計含隨機波動的時變參數(shù)回歸模型的貝葉斯半?yún)?shù)方法, 并將其應(yīng)用在通貨膨脹問題的分析中.Chow等[12]在非參數(shù)部分服從Dirichlet 過程時利用數(shù)值方法計算出貝葉斯估計值.Kim等[13]構(gòu)造了正態(tài)先驗下參數(shù)與非參數(shù)部分的一種貝葉斯估計.

本工作針對建立的半?yún)?shù)順序回歸模型(5), 以非參數(shù)部分隨機過程的有限維分布作為先驗, 構(gòu)造了參數(shù)部分與非參數(shù)部分的貝葉斯估計量.基于Kim等[13]提出的方法, 在正態(tài)情形下推導(dǎo)出估計量的解析表達式, 最后通過仿真模擬與實證分析評價貝葉斯估計量的表現(xiàn).

1 貝葉斯估計

根據(jù)樣本{(xi,ui,zi),i=1,2,··· ,n}構(gòu)造模型(5)中α,β,σ2,f(·)的貝葉斯估計, 從而在新狀態(tài){(xi,ui),i=n+1,n+2,··· ,n+t}下估計zi的值.下面分3 步構(gòu)造α,β,σ2,f(·)的貝葉斯估計量: ①構(gòu)造pj(x,u)的貝葉斯估計; ②根據(jù)式(1)計算γj(x,u); ③結(jié)合先驗分布推導(dǎo)α,β,σ2,f(·)的后驗分布, 從而構(gòu)造貝葉斯估計量.

本工作僅就正態(tài)情形給出推導(dǎo)過程, 但估計方法適用于一般分布.在正態(tài)情形下, 可寫出估計量的解析表達式; 在非正態(tài)情形下, 可利用蒙特卡羅-馬爾科夫(Monte Carol-Markov chain, MCMC)等方法進行數(shù)值計算.

當(dāng)(xi,ui)處的觀察值較少時,pij=pj(xi,ui)的極大似然估計值可能會取不合理的0 或1,因此首先對pij進行貝葉斯估計.取pi= (pi1,pi2,··· ,piK)的先驗為Dirichlet 分布, 響應(yīng)變量服從多項分布, 即

則后驗

因此,pij的貝葉斯估計[14]為

在(xi,ui)下, 模型(5)滿足因此, 在不計一個常數(shù)的情況下,{αj,1 ≤j

當(dāng)n=2,t=1,K=3,k=1 時, 模型(5)的樣本模型為

T=為n階單位矩陣的每一行重復(fù)K?1 行,f(u)=(f(u1),f(u2),··· ,f(un+t))′,ε=(ε11,··· ,ε1,K?1,ε21,··· ,εn,K?1)′, 則

為對t個新樣本預(yù)測其分類, 需估計參數(shù)β?和函數(shù)f(·).但f(·)是一個隨機過程, 有無窮多個參數(shù), 因此無法直接估計函數(shù)f(·), 但可估計f(·)在un+1,un+2,··· ,un+t處的函數(shù)值下面對正態(tài)情形構(gòu)造參數(shù)部分與非參數(shù)部分的貝葉斯估計量.

在擾動項ε服從正態(tài)分布的情況下, 有

為了構(gòu)造參數(shù)和非參數(shù)部分的貝葉斯估計量, 設(shè)參數(shù)部分β?的先驗服從多維正態(tài)分布,函數(shù)f(·)服從高斯過程(Gaussian process, GP), 則f(·)在u處的函數(shù)值服從高斯過程的有限維分布N(·,·).σ的先驗分布取Jeffreys 先驗, 即

以下推導(dǎo)過程中將f(u),f0(u)分別簡記為f,f0.

定理1記的后驗分布為正態(tài)-逆伽馬分布.

證明 根據(jù)貝葉斯理論, 在先驗(15), (16)和(17)滿足的情況下, (β?,f,σ)的后驗分布為

式中:

式(18)的推導(dǎo)主要用到矩陣二次型的化簡.

注意到, 對任意v ?=0∈RK?1+k+n+t,

式(23)中的“>”號由v1,v2不同時為0 保證.因此, 矩陣A正定, 式(19)中的A?1存在.

由式(18)可知,θ,σ|Y,x,u服從正態(tài)-逆伽馬分布(normal-inverse-gamma distribution,NIG), 即分別對σ2,θ積分, 可求得θ,σ2的邊際分布.

推論1后驗分布

在二次損失下,θ的貝葉斯估計為后驗期望估計即可得^θ?,f.

推論2后驗分布σ2|Y,x,u服從逆伽馬分布(inverse-gamma distribution, IGa), 即

在二次損失下,σ2的貝葉斯估計為后驗期望估計

對新的樣本(x,u), 根據(jù)

預(yù)測

2 數(shù)值模擬

觀察在不同先驗分布、樣本點處的觀察次數(shù)下^β?,f,^σ2的表現(xiàn).設(shè)回歸模型樣本量n= 70, 預(yù)測t= 30 步, 研究不同先驗分布, 不同觀察次數(shù)下的表現(xiàn).先驗設(shè)置中, 有

表1 N =1 000 時的參數(shù)估計Table 1 Estimation of parameters when N =1 000

(2) 為探究先驗分布對估計量效果的影響, 以f(u)為例, 取先驗均值f0(u)=16u(0.5?u),重復(fù)觀察次數(shù)m=3(隨機種子為123), 由推論1 估計得到的與真實的f(u) = sin(2πu)的比較如圖1(b)所示.

圖1 非參數(shù)部分f 與估計的比較Fig.1 Comparing of nonparametric f and estimation

圖1(a)表明, 在f(u)的先驗分布的均值f0(u) = 0 時,會對先驗分布有一定修正, 估計值趨向于真實的f(u).每個樣本點處的重復(fù)觀察次數(shù)m越多,的表現(xiàn)越好.表1 中與真實的α2,α3差異不大,與真實的β,σ2有一定差異, 但會隨著重復(fù)觀察次數(shù)m的增大越來越接近真實值.從,的1 000 次模擬的方差來看, 每個樣本點處的重復(fù)觀察次數(shù)m越多,,,,的估計值越穩(wěn)定.

對比圖1(b)所示的f0(u)=16u(0.5?u)與f(u)=sin(2πu)的曲線,在區(qū)間[0,0.7]內(nèi)f0(u)接近f(u),估計得到的接近真實值;在區(qū)間[0.7,1]內(nèi)f0(u)偏離f(u)很大,端點處f0(1)=?8,估計得到的很大程度上修正了先驗分布與真實值之間的偏離,端點處接近?2.3.因此在先驗分布接近真實值時,的估計效果較理想, 當(dāng)先驗分布偏離真實值時,也有不錯的表現(xiàn).同理, 可探究先驗分布的選取對的影響.

3 實證分析

對收入等級預(yù)測問題建立一個半?yún)?shù)順序回歸模型, 以食品、衣著、居住、家庭設(shè)備及用品、交通通信、文教娛樂、醫(yī)療保健和其他共8 項主要生活性消費支出[15]的占比作為解釋變量, 根據(jù)消費結(jié)構(gòu)預(yù)測收入等級.在很多經(jīng)濟調(diào)查中, 真實收入的收集非常困難, 本工作根據(jù)消費習(xí)慣來預(yù)測收入等級, 有助于基于收入進行數(shù)據(jù)分析.

根據(jù)國家統(tǒng)計局2002—2012 年人均消費支出占比的統(tǒng)計數(shù)據(jù)(http://data.stats.gov.cn/easyquery.htm?cn=C01), 對不同收入等級居民的消費結(jié)構(gòu)進行實證分析.8 項消費支出占比用每項的消費支出/總消費支出計算得到.農(nóng)村居民家庭的收入等級采取國家統(tǒng)計局的五等份劃分法, 分為低收入、中等偏下收入、中等收入、中等偏上收入、高收入家庭(分別用1, 2, 3, 4,5 表示); 城鎮(zhèn)居民家庭的收入等級在國家統(tǒng)計局的劃分法基礎(chǔ)上稍作改動, 將低收入與較低收入家庭劃分為低收入家庭, 較高收入與最高收入家庭劃分為高收入家庭, 這樣城鎮(zhèn)居民家庭由原來的7 個收入等級改為5 個收入等級.

選取2002—2011 年的數(shù)據(jù)(共120 組)作為訓(xùn)練樣本, 用半?yún)?shù)順序回歸模型(5)擬合8 項人均消費支出占比與家庭收入等級之間的關(guān)系.令u為8 項消費支出占比,

則由模型

從表2 的實驗結(jié)果可以看到, 對收入五等級問題的12 組樣本預(yù)測中, 預(yù)測準(zhǔn)確率為58.33%.預(yù)測錯誤的5 組樣本中, 有4 組樣本的預(yù)測等級與實際等級只相差1 個等級.

表2 半?yún)?shù)順序回歸模型的外推收入等級Table 2 Extrapolation income level of semiparanetric ordinal regression mode

已有研究表明, 在諸多影響消費需求的因素中, 收入水平始終是影響消費需求的最重要因素.而本工作實證分析的結(jié)果表明, 以居民家庭的消費結(jié)構(gòu)為解釋變量可以較準(zhǔn)確地預(yù)測收入等級, 因此消費結(jié)構(gòu)反過來也反映了家庭的收入情況.以8 項生活性消費支出考慮家庭的收入等級的方式, 比采用恩格爾系數(shù)更為全面.綜上可知, 在小樣本情況下, 貝葉斯估計由于利用了先驗信息, 往往有更好的估計效果.

4 結(jié)束語

本工作在正態(tài)情形下構(gòu)造了半?yún)?shù)回歸模型中參數(shù)與非參數(shù)部分的貝葉斯估計, 多次模擬結(jié)果表明, 在先驗均值均取0 時仍然有不錯的估計效果, 在先驗分布接近真實值時, 估計效果會更理想.

相比模型(3), 本工作建立的模型是基于比例優(yōu)勢模型和半?yún)?shù)模型(4)所形成的更一般的半?yún)?shù)順序回歸模型.模型(3)只考慮了被解釋變量與非參數(shù)部分變量u的關(guān)系, 無法反映參數(shù)部分變量x的影響, 而本模型建立了解釋變量與被解釋變量的半?yún)?shù)關(guān)系, 適用范圍更廣泛.例如, 當(dāng)x為外生變量時, 本模型可反映不同情況下u對被解釋變量的影響.在根據(jù)消費結(jié)構(gòu)預(yù)測家庭收入等級的實例中, 本工作將線性部分取為啞變量, 使模型對城鎮(zhèn)居民和農(nóng)村居民的收入等級預(yù)測問題都適用, 同時線性部分的參數(shù)反映了相同收入等級的城鎮(zhèn)居民家庭與農(nóng)村居民家庭的消費結(jié)構(gòu)差異.以2002—2011 年間的數(shù)據(jù)作為樣本, 對2012 年的數(shù)據(jù)作預(yù)測,即在外推情形下, 對五等級預(yù)測問題的預(yù)測準(zhǔn)確率達到58.33%, 預(yù)測錯誤的樣本中, 預(yù)測等級與實際等級大多只相差1 個等級.