零一膨脹二項(xiàng)混合回歸模型的改進(jìn)EM算法

呂敏紅,閆奕榮,吳成晶

(1.西安航空學(xué)院理學(xué)院,陜西 西安 710077; 2.西安交通大學(xué)經(jīng)濟(jì)與金融學(xué)院,陜西 西安 710049)

0 引 言

在現(xiàn)實(shí)生活中，數(shù)據(jù)來源于單一總體的假設(shè)往往是不成立的，如果仍然用同一個(gè)模型對這些“非同質(zhì)”的數(shù)據(jù)進(jìn)行處理，就會(huì)造成數(shù)據(jù)的散度偏大?；旌匣貧w模型是指多個(gè)混合成分按照一定的比例進(jìn)行混合而成的回歸模型，這樣就可以將具有混合分布的數(shù)據(jù)有效地模型化。近年來，各類混合模型在基因研究、金融和貿(mào)易等眾多研究領(lǐng)域迅猛發(fā)展，對混合模型的研究受到了越來越多統(tǒng)計(jì)學(xué)者的重視。Tchetgen等人[1]對廣義線性混合效應(yīng)模型的診斷檢驗(yàn)問題做了研究，Li等人[2]基于EM算法考慮了正態(tài)混合模型的假設(shè)檢驗(yàn)問題，Zhu等人[3]基于似然估計(jì)的方法研究了混合模型的漸進(jìn)理論，Lin等人[4]研究了偏斜正態(tài)分布下有限混合模型的參數(shù)估計(jì)問題，Xiang等人[5]研究了兩成分泊松混合回歸模型的影響診斷問題，陳家驊等人[6]研究了有限混合Von Mises模型的極大似然估計(jì)。

計(jì)數(shù)數(shù)據(jù)廣泛存在于醫(yī)療、保險(xiǎn)和金融等多個(gè)研究領(lǐng)域，最早這類數(shù)據(jù)通常會(huì)用泊松分布或二項(xiàng)分布等經(jīng)典的分布進(jìn)行模擬，但是在實(shí)際問題中，觀測數(shù)據(jù)可能存在過多的0或1，這時(shí)如果還用經(jīng)典分布就會(huì)造成數(shù)據(jù)散度偏大。例如在車險(xiǎn)理賠問題中，由于本年度理賠的次數(shù)會(huì)影響到下年度的保險(xiǎn)費(fèi)用，所以一般在較小的事故發(fā)生時(shí)，車主都會(huì)自掏腰包進(jìn)行維修，所以從保險(xiǎn)公司獲得的車輛發(fā)生交通事故次數(shù)的數(shù)據(jù)就會(huì)出現(xiàn)零觀測值過多的情況，這時(shí)就需要考慮零膨脹回歸模型。自從Cohen[7]第一次考慮了零膨脹數(shù)據(jù)，提出調(diào)整的泊松回歸模型后，對具有零膨脹的計(jì)數(shù)數(shù)據(jù)已經(jīng)有了很多研究，Lambert[8]正式提出了零膨脹泊松回歸模型，F(xiàn)ahrmeir等人[9]對零膨脹的可加模型進(jìn)行了研究，Ridout等人[10]對零膨脹回歸的貝葉斯方法進(jìn)行了深入的研究，Jansakul等人[11]研究了零膨脹模型的score檢驗(yàn)問題，Melkersson等人[12]在對看牙醫(yī)次數(shù)的數(shù)據(jù)進(jìn)行研究時(shí)，提出了零一膨脹泊松分布，Wang等人[13]對零一膨脹分布的性質(zhì)進(jìn)行了研究。盡管目前對于混合回歸模型和零一膨脹模型已經(jīng)有了很多方面的研究，但仍有較大的探索空間，特別是零一膨脹的混合模型還未涉及。

本文首先對零一膨脹二項(xiàng)混合回歸模型的參數(shù)建立了極大似然估計(jì)，然后由于傳統(tǒng)的EM算法只能使得估計(jì)收斂到局部極值，所以本文提出了MCEM算法對EM算法進(jìn)行改進(jìn)，使得零一膨脹二項(xiàng)混合回歸模型能夠找到全局最優(yōu)解。最后再通過一個(gè)模擬研究驗(yàn)證該算法的有效性。

1 混合回歸模型

混合回歸模型[4]是指一個(gè)回歸模型是由幾個(gè)混合成分按照一定比例混合而組成的。也就是說，混合回歸模型是將因變量服從混合分布的情形進(jìn)行模型化，從而對總體的“非同質(zhì)性”進(jìn)行描述。在實(shí)際中，混合成分的個(gè)數(shù)大多是有限的，把這類模型稱之為有限混合回歸模型，具體形式如下：

(1)

于是相應(yīng)的對數(shù)似然函數(shù)為：

(2)

傳統(tǒng)的極大似然估計(jì)是關(guān)于θ對式(2)求極大值，但是由于其結(jié)構(gòu)比較復(fù)雜，直接求解是相當(dāng)困難的。受數(shù)據(jù)添加思想的啟示[13]，首先引入指示向量Wi=(wi1,wi2,…,wip)T，若yi來自第j個(gè)混合成分，記wij=1，否則wij=0。這樣就有：

p(Wi=wi)=π1wi1π2wi2…πpwip

進(jìn)一步可得：

這樣基于完全數(shù)據(jù)的對數(shù)似然函數(shù)為：

(3)

2 零一膨脹二項(xiàng)混合回歸模型

零一膨脹回歸模型是計(jì)數(shù)數(shù)據(jù)中0和1都過多的計(jì)數(shù)模型，假設(shè)φ0表示數(shù)據(jù)中過多的0所占數(shù)據(jù)的比例，φ1表示過多的1所占的比例，而其他數(shù)據(jù)的取值服從某種離散分布，它們是按照一定的比例進(jìn)行混合：

(4)

其中，f(y)表示來自某種離散分布，如二項(xiàng)分布、泊松分布等,φ2=1-φ0-φ1。

從式(4)可以看出數(shù)據(jù)集中的0來自2個(gè)部分，即第一部分的0和第三部分的0，1也是同樣的道理來自2個(gè)部分。當(dāng)φ1=0,φ2=0時(shí)，數(shù)據(jù)完全來自某種離散分布；當(dāng)φ1=0,φ2≠0時(shí)，數(shù)據(jù)只在0處發(fā)生了膨脹，式(4)就簡化為零膨脹模型。

如果式(4)中的某種離散分布為二項(xiàng)分布時(shí)，即f(y)～B(m,y)，模型便為零一膨脹二項(xiàng)分布(ZOIB),具體形式為：

(5)

其中，φ2=1-φ0-φ1，0來自非泊松分布中的0和二項(xiàng)分布中的0，1也是同樣的道理。下面對模型(5)的參數(shù)部分引入?yún)f(xié)變量X和Z，這樣便得到了ZOIB的具體形式：

(6)

其中，β和γ0,γ1是回歸系數(shù)向量。

若式(1)中的每個(gè)混合成分都是ZOIB，則得到零一膨脹二項(xiàng)有限混合回歸模型。特別地，當(dāng)p=2時(shí)，則說明模型具有2個(gè)混合成分。本文重點(diǎn)研究具有2個(gè)混合成分的ZOIB混合回歸模型，具體形式如下：

(7)

3 ZOIB有限混合回歸模型的極大似然估計(jì)

(8)

其中，f(yi|θk)由公式(5)和公式(6)給出。通過式(8)對混合系數(shù)π求偏導(dǎo)，得到相應(yīng)的估計(jì)方程為：

但是在實(shí)際計(jì)算中上式是很難求解的，即使得到混合系數(shù)π的相應(yīng)估值，也可能超出π要求的范圍0<π<1。并且傳統(tǒng)的極大似然估計(jì)是關(guān)于式(8)求極大值，但是由于其結(jié)構(gòu)比較復(fù)雜，直接求解是相當(dāng)困難的。所以首先引入指示變量wi，若yi來自第一個(gè)混合成分，記wi=1，否則wi=0。這樣就可以給出完全數(shù)據(jù)集Ycom=(Y0,wi)，其中Y0=(yi,Xi,Zi)為觀測數(shù)據(jù)。

基于完全數(shù)據(jù)的對數(shù)似然函數(shù)為：

(1-wi)log f(yi|θ2)]}

(9)

可以看出式(9)中的項(xiàng)數(shù)比式(8)中的項(xiàng)數(shù)有所增加，但是添加的潛在變量是線性的，計(jì)算時(shí)相對容易。本文將基于公式(9)建立參數(shù)極大似然估計(jì)的算法。

4 加速M(fèi)CEM算法

EM算法最初由Dempster等人于1977年首次提出，它是當(dāng)數(shù)據(jù)存在缺失時(shí)常用的一種迭代算法，由于操作方便且穩(wěn)定，所以實(shí)用性很強(qiáng)。但是傳統(tǒng)的EM算法[14-15]只能使得估計(jì)收斂到局部極大值，而MCEM算法會(huì)大大降低收斂速度。下面提出加速M(fèi)CEM算法對傳統(tǒng)的EM算法和MCEM進(jìn)行修正[16]。具體就是將MCEM算法與Newton-Raphson算法結(jié)合，利用蒙特卡羅法解決高維空間的積分和優(yōu)化問題，通過某種實(shí)驗(yàn)的方法來估算隨機(jī)變量的期望。該算法具有二次收斂速度，從而使其保留MCEM算法的優(yōu)點(diǎn)，卻改進(jìn)了MCEM算法的收斂速度。具體算法包括如下3個(gè)步驟：

1)E1步。

從條件分布f(wi|θ(t),Y0)中隨機(jī)地抽取ki個(gè)樣本。其中θ(t)表示第t次迭代后θ估計(jì)的當(dāng)前值，f(wi|θ(t),Y0)表示在給定當(dāng)前θ的估計(jì)值和觀測值Y0下，潛在數(shù)據(jù)wi的條件預(yù)測分布。

2)E2步。

運(yùn)用步驟1中抽取的ki個(gè)樣本，令:

(10)

3)M步。

(11)

5 模擬研究

表1 2種算法下的參數(shù)估計(jì)

paraEMMCEMparaEMMCEMβ100.82140.8207β200.69790.6981β110.19780.1974β210.31790.3175γ1,000.41070.4105γ2,000.23050.2306γ1,010.20780.2075γ2,010.32690.3271γ1,100.29870.2989γ2,100.50420.5039γ1,110.10320.1033γ2,110.31070.3108

表2 2種算法下的迭代速度

算法迭代次數(shù)加速M(fèi)CEM27EM39

從表1計(jì)算結(jié)果中容易看出加速M(fèi)CEM算法在計(jì)算的可行性方面表現(xiàn)良好，說明其在E步中使用Monte Carlo模擬所犧牲的精度，在M步中又被補(bǔ)償了回來。從表2可以看出，在收斂速度方面，加速M(fèi)CEM算法的收斂速度也優(yōu)于EM算法。這表明在零一膨脹二項(xiàng)混合回歸模型的參數(shù)估計(jì)問題上，加速M(fèi)CEM算法無論在計(jì)算的可行性方面，還是在收斂速度方面，都有著優(yōu)于EM算法的表現(xiàn)。

6 結(jié)束語

本文首先對存在過多0和1的觀測數(shù)據(jù)提出了零一膨脹回歸模型，由于在現(xiàn)實(shí)生活中，數(shù)據(jù)來源于單一總體的假設(shè)往往是不成立的，所以進(jìn)一步引入了零一膨脹二項(xiàng)混合回歸模型。針對EM算法通常會(huì)使得估計(jì)收斂到局部最優(yōu)解上的缺陷，提出了加速M(fèi)CEM算法,對具有有限混合成分的ZOIB的參數(shù)進(jìn)行估計(jì)。最后通過模擬研究說明該方法的有效性。但是本文并未對混合比例考慮回歸，這將是今后的研究重點(diǎn)。

[1] Tchetgen E J, Coull B A. A diagnostic test for the mixing distribution in a generalised linear mixed model[J]. Biometrika, 2006,93(4):1003-1010.

[2] Li Pengfei, Chen Jiahua. Testing the oder of a finite mixture[J]. Journal of American Statistical Association, 2010,105(491):1084-1092.

[3] Zhu Hongtu, Zhang Heping. Hypothesis testing in mixture regression models[J]. Journal of Royal Statistical Society Series B, 2004,66(1):3-16.

[4] Lin T I, Lee J C, Yen S Y. Finite mixture modeling using the skew normal distribution[J]. Statistica Sinica, 2007,17(3):909-927.

[5] Xiang Liming, Yau K K W, Lee A H, et al. Influence diagnostics for two-component Poisson mixture regression models: Applications in public health[J]. Statistics in Medicine, 2005(19):3053-3071.

[6] 陳家驊,李鵬飛,譚鮮明. 混合Von Mises模型的參數(shù)估計(jì)[J]. 系統(tǒng)科學(xué)與數(shù)學(xué), 2007,27(1):59-67.

[7] Cohen A. Estimation of the Poisson parameter from truncated samples and from censored samples[J]. Journal of American Statistical Association, 1954,49(265):158-168.

[8] Lambert D. Zero-inflated Poisson regression with an application to defects in manufacturing[J]. Technometrics, 1992,34(1):1-14.

[9] Fahrmeir L, Echavarria L O. Structured additive regression for overdispersed and zero-inflated count data[J]. Applied Stochastic Models in Business and Industry, 2006,22(4):351-369.

[10] Ridout M, Hinde J, Demetrio C G B. A score test for testing zero-inflated Poisson regression model against zero-inflated negative binomial alternatives[J]. Biometrics, 2001,57(1):219-223.

[11] Jansakul N, Hinde J P. Score tests for zero-inflated Poission models[J]. Computational Statistics and Data Analysis， 2002,40(1):75-96．

[12] Melkersson M, Olsson C. Is Visiting the Dentist A Good Habit? Analyzing Count Data with Excess Zeros and Excess Ones[D]. Sweden: Umer University, 1999.

[13] Zhang Chi, Tian Guoliang, Wang Kai. Properties of the zero-and-one inflated Poisson distribution and statistical inference methods[J]. Statistics and Its Interface, 2016,9(1):11-32.

[14] Biernacki C. Initializing EM using the properties of its trajectories in Gaussian mixtures[J]. Statistics and Computing, 2004,14(3):267-279.

[15] 呂敏紅,閆奕榮,楊青. 零一膨脹泊松回歸模型的EM算法改進(jìn)[J]. 河南科學(xué), 2017，35(7):1037-1041.

[16] 羅季. Monte Carlo EM加速算法[J]. 應(yīng)用概率統(tǒng), 2008,24(3):312-318.

[17] 盧玉桂，韋新星，趙麗棉. 多層線性模參數(shù)估計(jì)的MCEM算法[J]. 數(shù)學(xué)的實(shí)踐與認(rèn)識(shí), 2016,46(11):225-230.