混合非參數(shù)回歸的貝葉斯推斷

李道揚(yáng), 何幼樺

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

混合回歸模型是研究在樣本服從的總體分布是不同分布混合時(shí), 解釋變量與被解釋變量關(guān)系的統(tǒng)計(jì)模型, 由Goldfeld 等[1]在1973 年提出. 之后, 很多學(xué)者對(duì)混合回歸的問(wèn)題展開(kāi)了研究. Wedel 等[2]于1995 年將混合線性回歸模型推廣到混合廣義線性回歸模型; 1999 年,Gaffney 等[3]使用期望最大化(expectation maximisation, EM)算法解決了混合回歸中的一般參數(shù)估計(jì)問(wèn)題; Song 等[4]在2014 年提出了基于拉普拉斯分布的混合穩(wěn)健回歸模型等. 對(duì)于混合非參數(shù)回歸問(wèn)題, 也有一些學(xué)者進(jìn)行了研究. Huang 等[5]在2013 年使用基于樣本局部似然函數(shù)的EM 算法解決了混合非參數(shù)回歸問(wèn)題; Wu 等[6]在2016 年提出了基于半?yún)?shù)模型的混合分位數(shù)回歸模型; 胡燁[7]于2017 年在混合泊松回歸模型的基礎(chǔ)上進(jìn)一步研究了半?yún)?shù)混合泊松回歸的問(wèn)題等. 在混合回歸模型的應(yīng)用方面, Fraley 等[8]在1998 年將混合回歸模型應(yīng)用于糖尿病患者與健康人的分類和密度估計(jì)問(wèn)題; Gaffney 等[3]于1999 年將混合回歸模型應(yīng)用于視頻中的動(dòng)作識(shí)別問(wèn)題等.

本工作針對(duì)一般的混合非參數(shù)回歸問(wèn)題, 給出了一種基于貝葉斯框架的方法. 首先, 對(duì)回歸中的每個(gè)非參數(shù)混合成分用一個(gè)隨機(jī)過(guò)程的有限維分布族作為先驗(yàn), 同時(shí)分別構(gòu)造混合比例、隨機(jī)誤差的方差和非參數(shù)混合成分的貝葉斯估計(jì), 并通過(guò)馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo, MCMC)法抽樣來(lái)進(jìn)行后驗(yàn)推斷; 然后, 通過(guò)數(shù)值模擬來(lái)分析混合非參數(shù)回歸模型的貝葉斯推斷方法的有效性; 最后, 應(yīng)用該推斷方法解決了蚜蟲(chóng)數(shù)量與受感染煙草植物數(shù)量關(guān)系的實(shí)際問(wèn)題.

1 混合非參數(shù)回歸的貝葉斯估計(jì)

考慮有K個(gè)非參數(shù)回歸的混合, 回歸的解釋變量為x, 被解釋變量為y, 第k個(gè)非參數(shù)混合成分gk(x)的比例為αk, 滿足αk ∈[0,1],αk= 1,隨機(jī)誤差εk ～N(0,σ2k),則混合非參數(shù)回歸模型(以概率αk產(chǎn)生)為

在混合非參數(shù)回歸問(wèn)題中, 對(duì)任意給定的x, 需要給出其每一個(gè)混合成分下回歸的估計(jì)值k(x)、混合比例的估計(jì)值k以及隨機(jī)誤差εk方差的估計(jì)值.

混合非參數(shù)回歸模型中被解釋變量y的分布函數(shù)一般為

式中:Φ(u) 為標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù), 記標(biāo)準(zhǔn)正態(tài)分布的概率密度函數(shù)為φ(u);g(x) =(g1(x),g2(x),··· ,gK(x)),α= (α1,α2,··· ,αK),σ= (σ1,σ2,··· ,σK). 在本模型中, 假定gk(x) 是非參數(shù)的, 而αk和εk與x無(wú)關(guān). 如果gk(x) 是線性模型, 那么模型(1) 就退化為Goldfeld 等[1]提出的混合線性模型, 也就是說(shuō), 本工作中的混合非參數(shù)回歸模型是混合線性回歸模型中的線性回歸在非參數(shù)回歸上的一種推廣.

設(shè){(xi,yi),i= 1,2,··· ,n}是總體的一組隨機(jī)樣本,U= (xn+1,xn+2,··· ,xn+t)′是待預(yù)測(cè)點(diǎn)的向量. 記X=(x1,x2,··· ,xn)′,Y=(y1,y2,··· ,yn)′,T=[X′,U′]′. 那么可以得到樣本的聯(lián)合密度:

式中:φσk(u)=).

在本工作中, 引入中間變量zik來(lái)簡(jiǎn)化式(2)的復(fù)雜度. 如果第i個(gè)樣本(xi,yi)屬于第k類, 那么就記zik=1; 否則zik=0, 也就是說(shuō)zik ∈{0,1}, 且有

記z=(zik)n×K, 則樣本的聯(lián)合密度函數(shù)可化簡(jiǎn)為

為了構(gòu)造非參數(shù)成分gk(x)、混合比例αk和隨機(jī)誤差方差σ2k的貝葉斯估計(jì), 需要分別給定其先驗(yàn)分布. 由于gk(x)是一個(gè)函數(shù), 因此其先驗(yàn)需要使用一個(gè)隨機(jī)過(guò)程的有限維分布族來(lái)刻畫(huà). 假設(shè)其服從高斯過(guò)程GP(gk0(·),·), 由于只關(guān)心gk(x)在T上的估計(jì), 因此gk(x)在T上的先驗(yàn)只需用高斯過(guò)程的有限維分布來(lái)描述, 即N(gk0(T),τ-1Σk), 其中Σk是階數(shù)為(n+t) 的協(xié)方差矩陣,τ是衡量先驗(yàn)分布信任程度的超參數(shù); 考慮到混合比例αk的先驗(yàn)需滿足αk0∈[0,1],αk0= 1, 取Dirichlet 分布作為混合比例聯(lián)合密度的先驗(yàn); 對(duì)于隨機(jī)擾動(dòng)的方差, 則將參數(shù)為θk和βk的逆伽馬分布作為其先驗(yàn), 即

式中:τ,g10(T),g20(T),··· ,gK0(T),Σ1,Σ2,··· ,ΣK,α10,α9,··· ,αK0,θ1,θ2,··· ,θK,β1,β2,··· ,βK均為先驗(yàn)超參數(shù).

定理 在先驗(yàn)滿足式(3)～(5) 的情況下, 關(guān)于α,gk(T)和σk的后驗(yàn)分布則有如下結(jié)論.

結(jié)論1gk(T)|g-k,α,σ,z,X,Y,T ～N(μ*,Σ*),k=1,2,··· ,K, 其中

結(jié)論2α|z,g,σ,X,Y,T ～Dirichlet(α*), 其中

結(jié)論3σ2k|z,g,α,σ,X,Y,T ～I(xiàn)G(σ2k;θ*k,β*k),k=1,2,··· ,K, 其中

證明 根據(jù)貝葉斯定理可以得到各個(gè)參數(shù)的聯(lián)合后驗(yàn)分布.

式中:g=(g1(T),g2(T),··· ,gK(T));Zk=diag(zik),i=1,2,··· ,n,k=1,2,··· ,K.

根據(jù)聯(lián)合密度函數(shù)可以得到在給定其他參數(shù)的情況下gk(T)的后驗(yàn)分布:

引入輔助矩陣E=[In×n On×t], 那么上式可化簡(jiǎn)為

整理可得

結(jié)論1 得證.

同理對(duì)α可以得到后驗(yàn)分布:

從而α|z,g,σ,X,Y,T服從參數(shù)為α*的Dirichlet 分布. 結(jié)論2 得證.

同樣對(duì)于σ可以計(jì)算后驗(yàn)分布:

于是,π(|z,g,α,X,Y,T)～I(xiàn)G(),k=1,2,··· ,K, 其中

結(jié)論3 得證.

從定理中可以得出,g,α,σ的后驗(yàn)分布均是結(jié)合先驗(yàn)信息和樣本信息的結(jié)果. 特別地, 對(duì)于具體的gk(T), 其后驗(yàn)分布的均值可以看作先驗(yàn)信息與樣本信息的加權(quán)平均, 也就是說(shuō), 在同樣的樣本信息的情況下, 較好的先驗(yàn)分布能得到更準(zhǔn)確的后驗(yàn)估計(jì).

由于聯(lián)合后驗(yàn)分布的密度函數(shù)的結(jié)構(gòu)較為復(fù)雜, 因此使用MCMC 方法中Gibbs 抽樣分別得到g,α,σ的后驗(yàn)分布. 利用Gibbs 抽樣構(gòu)造MCMC 鏈的關(guān)鍵在于分別得到g,α,σ的滿條件分布, 然后根據(jù)滿條件分布依次輪換抽樣得到相應(yīng)的馬氏鏈. 而馬氏鏈的平穩(wěn)性, 使得把這些鏈的實(shí)現(xiàn)作為樣本來(lái)推斷所需的聯(lián)合分布得以實(shí)現(xiàn). 結(jié)合定理的結(jié)果, 可以得到各參數(shù)滿條件分布, 從而給出MCMC 抽樣算法.

(1) 設(shè)定先驗(yàn)分布的超參數(shù)τ,g10(T),g20(T),··· ,gK0(T),Σ1,Σ2,··· ,ΣK,α10,α20,··· ,αK0,θ1,θ2, ··· ,θK,β1,β2,··· ,βK和g,α,σ的初值g(0),α(0),σ(0), 以及MCMC 抽樣次數(shù).

(2) 抽樣第j步, 根據(jù)P(zik= 1|g(j-1),α(j-1),σ(j-1),X,Y) 分別計(jì)算zik= 1 的后驗(yàn)概率pik,k= 1,2,··· ,K,服從參數(shù)為(pi1,pi2,··· ,piK) 的多項(xiàng)分布, 即從～Multinomial(pi1,pi2,··· ,piK)抽取,i=1,2,··· ,n.

(3)根據(jù)gk(T)的滿條件分布π(gk(T)|,X,Y,T)抽取(T), 其中g(shù)-k=(g1(T),g2(T),··· ,gk-1(T),gk+1(T),··· ,gK(T)),k=1,2,··· ,K.

(4) 根據(jù)α的滿條件分布π(α|z(j),g(j),σ(j-1),X,Y,T)抽取α(j).

(5) 對(duì)σ2k, 根據(jù)滿條件分布π(σ2k|z(j),g(j),α(j),X,Y,T)抽取σ2(j)k ,k=1,2,··· ,K.

(6) 重復(fù)(2)～(5), 直至完成預(yù)設(shè)的抽樣次數(shù).

2 數(shù)值模擬

本工作通過(guò)討論混合非參數(shù)回歸模型的貝葉斯推斷方法在不同情況下的表現(xiàn), 以驗(yàn)證該方法的有效性. 下面分別從樣本量大小、各類回歸線的相對(duì)位置以及多種混合成分的情況這3 個(gè)角度來(lái)討論此推斷方法的效果. 為了更加直觀地顯示貝葉斯推斷方法的有效性, 采用Xiang[9]在2018 年提出的全局期望最大化(global EM, GEM)算法作為對(duì)照方法. GEM 算法是一種用于處理混合半?yún)?shù)混合回歸問(wèn)題的改進(jìn)類EM 算法, 同樣適合處理混合非參數(shù)回歸問(wèn)題.

2.1 樣本量不同的情況

這里, 主要討論考察樣本量大小對(duì)貝葉斯推斷方法估計(jì)精度的影響, 并與對(duì)照方法作比較. 為此, 本工作采取了一種2 個(gè)成分混合的模型, 考慮到產(chǎn)生樣本的模型需要一定的波動(dòng)性來(lái)反映貝葉斯推斷方法和對(duì)照方法的擬合能力, 故采用如下方式設(shè)置產(chǎn)生樣本:

為了更直觀地比較貝葉斯推斷方法和對(duì)照方法的精度, 本工作采用平均根誤差(square root of the average squared error, RASE)作為衡量估計(jì)精度的指標(biāo):

式中:xn+i,i= 1,2,··· ,t為待預(yù)測(cè)樣本點(diǎn);N為實(shí)驗(yàn)的重復(fù)次數(shù);(xi)為第l次實(shí)驗(yàn)中g(shù)k(xi)的估計(jì)值. 采用均方根誤差(root of mean square error, RMSE)作為衡量混合比例α和隨機(jī)誤差項(xiàng)ε方差σ2的估計(jì)精度指標(biāo):

式中:,分別為第l次實(shí)驗(yàn)中α,σ的估計(jì)值. 在貝葉斯推斷方法中, 采用后驗(yàn)均值作為回歸和參數(shù)的估計(jì)值. 在貝葉斯推斷方法中設(shè)置先驗(yàn)及超參數(shù):

式中:g10(x) =(2;g20(x) )=-2;τ-1= 0.1;α10= 0.5;α20= 0.5;σ0= 0.1;θ=-1.5;β=0;(Σ)ij=exp;h為經(jīng)驗(yàn)帶寬.

本工作使用貝葉斯推斷方法和GEM 算法分別在樣本量n= 50,100,150,500 時(shí)計(jì)算100個(gè)預(yù)測(cè)點(diǎn)的結(jié)果, 并對(duì)各種情況重復(fù)100 次實(shí)驗(yàn), 計(jì)算預(yù)測(cè)點(diǎn)的誤差結(jié)果如表1 所示.

表1 樣本量不同情況下的估計(jì)誤差Table 1 Estimation error in different sample sizes

觀察表1, 可以得到2 個(gè)結(jié)論: ①隨著樣本量的增加, 貝葉斯推斷方法和GEM 算法的精度都得到了提高; ②對(duì)比貝葉斯推斷方法和GEM 算法的結(jié)果可知, 在對(duì)α的估計(jì)方面, 貝葉斯推斷方法略好于GEM 算法. 而在回歸的精度和對(duì)σ的估計(jì)方面, 貝葉斯推斷方法顯著優(yōu)于GEM 算法, 這表明貝葉斯推斷方法能夠利用先驗(yàn)信息來(lái)增強(qiáng)模型的擬合和預(yù)測(cè)能力, 并且貝葉斯推斷方法即使在樣本量不大的情況下, 也能得到較為精確的結(jié)果

2.2 各類回歸線的相對(duì)位置

下面考察在各類回歸線的不同相對(duì)位置情況下混合非參數(shù)回歸的貝葉斯推斷方法的表現(xiàn). 當(dāng)樣本的類間距離較大時(shí), 這類問(wèn)題較容易解決, 因此本工作著重考慮類間距離較小時(shí)的情況.

當(dāng)?shù)玫饺鐖D1(a)所示的樣本時(shí), 2 條回歸線的形狀存在2 種可能: ①混合成分之間存在交叉的情況; ②2 個(gè)混合成分在x= 0 處幾乎相切但未交叉. 為討論上述情況下模型的表現(xiàn), 應(yīng)用貝葉斯推斷方法和GEM 算法下的所得結(jié)果如圖1(c)～(d)所示.

此時(shí), GEM 算法的混合非參數(shù)回歸只能得到一種情況, 即默認(rèn)曲線不交叉. 而通過(guò)貝葉斯推斷方法可以結(jié)合先驗(yàn)信息對(duì)模型有一個(gè)初步判斷, 并通過(guò)先驗(yàn)反映到模型上, 這樣貝葉斯推斷方法就可以根據(jù)不同的具體情況擬合2 種形狀(見(jiàn)圖1(c), (d)), 這一點(diǎn)是傳統(tǒng)的類EM 方法不能做到的. 從這個(gè)角度看, 貝葉斯推斷方法能夠滿足更一般的情況, 這也是貝葉斯推斷方法顯著的優(yōu)勢(shì).

圖1 樣本及不同方法下的回歸Fig.1 Sample and regression in different methods

從回歸的結(jié)果看, GEM 算法對(duì)樣本的擬合程度不高, 尤其是在兩側(cè)樣本邊緣處, 這也是非參數(shù)模型的共性問(wèn)題; 而在貝葉斯推斷方法的結(jié)果中, 無(wú)論是取交叉的先驗(yàn)還是不交叉的先驗(yàn), 總體回歸效果較好, 這說(shuō)明貝葉斯推斷方法能夠在樣本較少的地方結(jié)合先驗(yàn)信息來(lái)優(yōu)化回歸的結(jié)果, 這從側(cè)面印證了貝葉斯推斷方法的精度要高于GEM 算法的結(jié)論.

2.3 多種混合成分的情況

以3 種混合成分問(wèn)題為例, 驗(yàn)證模型在多種混合成分情況下的有效性(見(jiàn)圖2). 為測(cè)試貝葉斯推斷方法和GEM 算法在復(fù)雜情況下的適應(yīng)性樣本分別由線性函數(shù)、三角函數(shù)和二次函數(shù)產(chǎn)生, 然后混合而成, 其實(shí)驗(yàn)結(jié)果如圖2 所示. 圖2 中的樣本從上至下分別由線性函數(shù)、三角函數(shù)、二次函數(shù)產(chǎn)生.

圖2 3 種非參數(shù)混合成分下的回歸Fig.2 Mixtures of three nonparametric regressions

從圖2 中可以看出, 2 種方法都能處理多種成分混合回歸問(wèn)題, 而不是僅限于2 種混合成分的非參數(shù)回歸問(wèn)題. 在波動(dòng)性較強(qiáng)的回歸曲線中, GEM方法的表現(xiàn)較差, 因?yàn)閭鹘y(tǒng)的非參數(shù)方法只能通過(guò)核函數(shù)加權(quán)利用附近點(diǎn)的樣本信息; 而貝葉斯推斷方法不僅可以利用樣本信息,還可以結(jié)合先驗(yàn)信息, 這就使得貝葉斯方法具有更強(qiáng)的擬合和預(yù)測(cè)能力. 對(duì)本例來(lái)說(shuō), 給定的先驗(yàn)信息只包含了3 條曲線的區(qū)分度信息, 而并未包含曲線的形狀和趨勢(shì)信息. 結(jié)果表明, 即便如此貝葉斯推斷方法也能得到較好的結(jié)果.

3 實(shí)證分析

Boiteau 等[10]給出了一個(gè)關(guān)于蚜蟲(chóng)和感染煙草植物關(guān)系的生物學(xué)實(shí)驗(yàn)數(shù)據(jù), 該數(shù)據(jù)包含51 個(gè)獨(dú)立實(shí)驗(yàn)的結(jié)果. 在每個(gè)實(shí)驗(yàn)環(huán)境中設(shè)置12 株感染煙草植物和69 株健康煙草植物, 并在環(huán)境中釋放不同數(shù)量的蚜蟲(chóng). 在每組實(shí)驗(yàn)進(jìn)行24 h 后, 檢測(cè)先前健康植物的感染情況并記錄受感染的植物數(shù)量. 所記錄的釋放蚜蟲(chóng)數(shù)量與受感染植物數(shù)量樣本散點(diǎn)圖如圖3(a)所示.

隨著蚜蟲(chóng)數(shù)量的增加, 感染植物的數(shù)量與其不是單純的正向相關(guān), 而是圍繞在2 條回歸線周圍. 其中一條隨著蚜蟲(chóng)數(shù)量的增加而增加, 且增長(zhǎng)速度明顯大于線性增長(zhǎng); 另一條受感染煙草植物數(shù)量極少, 且對(duì)蚜蟲(chóng)數(shù)量的增加并不敏感. 這種情況是傳統(tǒng)的一條回歸曲線所不能處理的, 需要用到混合回歸. Boiteau 等[10]使用了混合線性回歸來(lái)處理此問(wèn)題, 擬合結(jié)果并不理想;Gr¨un 等[11]在2008 年用混合廣義線性回歸來(lái)擬合結(jié)果, 所得結(jié)果有一定的改進(jìn), 但因其對(duì)模型做了預(yù)設(shè)使得模型不具有一般性.

使用混合非參數(shù)回歸的貝葉斯推斷方法結(jié)果如圖3(b)所示. 模型同時(shí)估計(jì)了混合比例和隨機(jī)誤差的標(biāo)準(zhǔn)差的后驗(yàn)均值. 在位于上方的曲線中, 混合比例α1的后驗(yàn)均值為0.621 5, 標(biāo)準(zhǔn)差σ1為5.412 0; 在位于下方的曲線中, 混合比例α2的后驗(yàn)均值為0.378 5, 標(biāo)準(zhǔn)差σ2為1.111 5. 總的來(lái)說(shuō), 在所有樣本中約有62%的樣本圍繞上方的曲線波動(dòng), 此時(shí)感染植物的數(shù)量隨蚜蟲(chóng)數(shù)量增加而增長(zhǎng)較快, 隨機(jī)誤差的方差也較大; 在大約有38%的樣本中, 受感染植物的數(shù)量幾乎沒(méi)有, 且與釋放蚜蟲(chóng)數(shù)量的關(guān)系并不明顯, 隨機(jī)誤差的方差也較小. 可見(jiàn)這2 條回歸線表明此推斷方法可以較好地捕捉樣本的性質(zhì)和結(jié)構(gòu).

圖3 蚜蟲(chóng)數(shù)量與受感染植物數(shù)量的關(guān)系Fig.3 Relationship between the number of aphids and the number of infected plants

4 結(jié)束語(yǔ)

本工作基于貝葉斯的框架給出了一種混合非參數(shù)回歸的統(tǒng)計(jì)推斷方法. 首先, 分別給出非參數(shù)成分gk(x)、混合比例αk以及隨機(jī)誤差的方差σ2k的先驗(yàn)分布(其中對(duì)非參數(shù)成分gk(x)的先驗(yàn)以一個(gè)隨機(jī)過(guò)程的有限維分布族的形式給出); 然后, 引入中間變量來(lái)簡(jiǎn)化樣本的聯(lián)合密度函數(shù); 最后, 結(jié)合先驗(yàn)信息和樣本信息給出后驗(yàn)分布. 數(shù)值結(jié)果表明, 由于混合非參數(shù)回歸的貝葉斯推斷方法能夠利用先驗(yàn)信息, 故在適應(yīng)性、一般性和精度上都優(yōu)于GEM 算法, 并且即使在小樣本的情況下也能得到較為精確的結(jié)果. 最后蚜蟲(chóng)數(shù)量與受感染植物數(shù)量關(guān)系的實(shí)證分析表明, 混合非參數(shù)回歸的貝葉斯推斷方法能夠同時(shí)有效地處理聚類與回歸問(wèn)題, 且方法更具一般性.