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

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

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

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

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

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

在混合非參數(shù)回歸問題中, 對任意給定的x, 需要給出其每一個混合成分下回歸的估計(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無關(guān). 如果gk(x) 是線性模型, 那么模型(1) 就退化為Goldfeld 等[1]提出的混合線性模型, 也就是說, 本工作中的混合非參數(shù)回歸模型是混合線性回歸模型中的線性回歸在非參數(shù)回歸上的一種推廣.

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

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

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

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

為了構(gòu)造非參數(shù)成分gk(x)、混合比例αk和隨機(jī)誤差方差σ2k的貝葉斯估計(jì), 需要分別給定其先驗(yàn)分布. 由于gk(x)是一個函數(shù), 因此其先驗(yàn)需要使用一個隨機(jī)過程的有限維分布族來刻畫. 假設(shè)其服從高斯過程GP(gk0(·),·), 由于只關(guān)心gk(x)在T上的估計(jì), 因此gk(x)在T上的先驗(yàn)只需用高斯過程的有限維分布來描述, 即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); 對于隨機(jī)擾動的方差, 則將參數(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ù)貝葉斯定理可以得到各個參數(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é)論1 得證.

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

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

同樣對于σ可以計(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é)果. 特別地, 對于具體的gk(T), 其后驗(yàn)分布的均值可以看作先驗(yàn)信息與樣本信息的加權(quán)平均, 也就是說, 在同樣的樣本信息的情況下, 較好的先驗(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)作為樣本來推斷所需的聯(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) 對σ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ù)值模擬

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

2.1 樣本量不同的情況

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

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

式中:xn+i,i= 1,2,··· ,t為待預(yù)測樣本點(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 時計(jì)算100個預(yù)測點(diǎn)的結(jié)果, 并對各種情況重復(fù)100 次實(shí)驗(yàn), 計(jì)算預(yù)測點(diǎn)的誤差結(jié)果如表1 所示.

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

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

2.2 各類回歸線的相對位置

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

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

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

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

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

2.3 多種混合成分的情況

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

3 實(shí)證分析

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

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

使用混合非參數(shù)回歸的貝葉斯推斷方法結(jié)果如圖3(b)所示. 模型同時估計(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. 總的來說, 在所有樣本中約有62%的樣本圍繞上方的曲線波動, 此時感染植物的數(shù)量隨蚜蟲數(shù)量增加而增長較快, 隨機(jī)誤差的方差也較大; 在大約有38%的樣本中, 受感染植物的數(shù)量幾乎沒有, 且與釋放蚜蟲數(shù)量的關(guān)系并不明顯, 隨機(jī)誤差的方差也較小. 可見這2 條回歸線表明此推斷方法可以較好地捕捉樣本的性質(zhì)和結(jié)構(gòu).

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

4 結(jié)束語

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