截斷刪失數(shù)據(jù)下泊松分布參數(shù)的點估計

何朝兵 ，柴士改，劉華文

(1.安陽師范學(xué)院數(shù)學(xué)與統(tǒng)計學(xué)院，安陽455000;2.山東大學(xué)數(shù)學(xué)學(xué)院，濟南250100)

泊松分布是一種很重要的離散型分布，適于描述單位時間或空間內(nèi)隨機事件發(fā)生的次數(shù)，如某一服務(wù)設(shè)施在一定時間內(nèi)到達(dá)的人數(shù)，電話交換機接到呼叫的次數(shù)等.文獻(xiàn)［1］－［5］討論了泊松分布的參數(shù)估計. 近年來，對截斷刪失數(shù)據(jù)的研究較多［6－11］，但缺乏涉及泊松分布的報道. 本文主要利用EM 算法和MCMC 方法對截斷刪失數(shù)據(jù)下泊松分布壽命參數(shù)的點估計進(jìn)行了研究. 利用逆變換法和舍選法得到了產(chǎn)品的完全數(shù)據(jù)和參數(shù)的EM 迭代公式.把Gibbs 樣本的算術(shù)平均值作為參數(shù)的MCMC 估計.隨機模擬的估計效果較好，估計值比較穩(wěn)定，且精度也較高.

1 截斷刪失數(shù)據(jù)模型

設(shè)(X，Y，T)是離散型隨機變量，X 的分布函數(shù)和分布律分別為F(x，λ)=P(X≤x)、f(x，λ);Y 的分布函數(shù)和分布律分別為G(y)、g(y);T 的分布函數(shù)和分布律分別為H(t)、h(t)，假設(shè)X、Y、T 獨立.對于n個受試樣品，截斷刪失數(shù)據(jù)模型是:Zi≥Ti時可觀察到數(shù)據(jù)(Zi，Ti，δi)，而Zi＜Ti時無法得到觀察值，其中

下面求樣本的似然函數(shù).P(無樣本觀察值)=P(Zi＜Ti)=1－P(Xi≥Ti，Yi≥Ti)=

為方便研究，引入示性變量νi=I(min(Xi，Yi)≥Ti)(i=1，2，…，n).則基于數(shù)據(jù){(zi，ti，δi):vi=1，1≤i≤n}的似然函數(shù)為

2 泊松分布參數(shù)的點估計

2.1 完全數(shù)據(jù)似然函數(shù)

截斷刪失數(shù)據(jù)模型中，假設(shè)產(chǎn)品壽命X 服從參數(shù)為λ 的泊松分布，即X ～P(λ).

由式(1)易知似然函數(shù)很復(fù)雜. 為了得到完全數(shù)據(jù)，不妨添加缺損壽命數(shù)據(jù).具體如下:

νi=1，δi=0 時，添加Z1i=Xi=z1i.Z1i的分布律為

νi=0 時，添加Z2i=Xi=z2i.Z2i的分布律為

由于f(x，λ)是Xi的分布律，Xi與Z2i有相同的取值，且

所以可以利用舍選法隨機產(chǎn)生Z2i的值z2i.

令δ，ν，z，u1，u2分別表示由δi，νi，zi，z1i，z2i組成的向量，則完全數(shù)據(jù)似然函數(shù)為

2.2 EM 算法

EM 算法［12］是一種求后驗分布眾數(shù)的迭代方法，它每一步迭代由E 步(求期望)和M 步(極大化)組成.下面求λ 的EM 迭代公式.

取λ 的先驗分布為共軛先驗分布伽瑪分布Ga(α，β)，α，β 已知，即

則λ 的添加后驗分布為

假設(shè)在第m+1 次迭代開始時的估計值為λ(m).令U1、U2分 別 表 示Z1i、Z2i組 成 的 向 量. 簡 記

E 步:

當(dāng)νi=1，δi=0 時，

當(dāng)νi=0 時，利用Monte Carlo 方法計算Z2i的期望，稱為MCEM方法，它將E 步變?yōu)槿缦? 步:

(E1)利用舍選法從ψ2(k，λ(m))抽取l個隨機數(shù)k1，k2，…，kl;

(E2)計算1］，其中

M 步:

式(2)給出了由EM 算法得到的參數(shù)λ 的迭代公式.

2.3 貝葉斯估計

當(dāng)νi=1，δi=0 時，ψ1(z1i，λ)，其中z－1i={z1j:j≠i}.

當(dāng)νi=0 時λ)，其中z－2i={z2j:j≠i}.

參數(shù)λ 的滿條件分布為

對z1i，z2i，λ 的滿條件分布分別進(jìn)行Gibbs 抽樣［13］.設(shè)λ(j)(j=1，2，…，M1，…，M2)為λ 的容量為M2的Gibbs 樣本，第M1次以后抽樣收斂，把后面M2－M1個樣本的均值作為λ 的估計，即

EM 估計是后驗分布的眾數(shù)(MLE)，MCMC 估計是后驗分布的期望，它們差別不大.

3 隨機模擬

設(shè)X ～P(λ)，Y ～P(λ1)，T ～P(λ2)，取λ 的先驗分布為Ga(α，β)，對于超參數(shù)α、β 的確定應(yīng)該充分利用各種先驗信息，而利用先驗矩確定超參數(shù)是一種常用的方法.雖然進(jìn)行隨機模擬試驗時沒有先驗信息，選取α、β 時，卻可以遵循“利用先驗矩確定超參數(shù)”這一方法的思想，盡可能使先驗期望α/β與E(X)= λ 相差不要太大，例如若λ =12，不妨取α=7，β =0.6，此時α/β≈11.67. 樣本容量分別取n=30，50，100，200，300，500，800. 對于(λ，λ1，λ2)分別取3 組不同的參數(shù)進(jìn)行隨機模擬試驗，具體如下:

(1)(λ，λ1，λ2)=(9，4，7)，α =14，β =1.5，此時E(Y)＜E(T)＜E(X)，即λ1＜λ2＜λ.

(2)(λ，λ1，λ2)=(12，15，10)，α =7，β =0.6，此時E(T)＜E(X)＜E(Y)，即λ2＜λ ＜λ1.

(3)(λ，λ1，λ2)=(30，40，36)，α =86，β =3，此時E(X)＜E(T)＜E(Y)，即λ ＜λ2＜λ1.

EM 迭代的初始值取為λ(0)= λ /2;Gibbs 抽樣迭代的初始值取為λ(0)= λ +5，取M1=5 000，M2=10 000. 參數(shù)估計結(jié)果見表1.

n=300 時的迭代過程見圖1 ～圖6.

表1 參數(shù)λ 的EM 估計和MCMC 估計Table 1 EM estimation and MCMC estimation of parameter λ

圖1 n=300，λ =9 時λ 的MCEM 迭代過程Figure 1 MCEM iterations of λ with n=300，λ =9

圖2 n=300，λ =9 時λ 的Gibbs 抽樣迭代過程Figure 2 Gibbs sampling iterations of λ with n=300，λ =9

圖3 n=300，λ =12 時λ 的MCEM 迭代過程Figure 3 MCEM iterations of λ with n=300，λ =12

由表1 看出，EM 估計與MCMC 估計的差別很小. λ =9 時估計的相對誤差不超過10%，λ =10，λ =30 時估計的相對誤差大部分不超過4%. λ =9 時，由于E(Y)＜E(T)＜E(X)，截斷刪失模型的原理使得此組參數(shù)下的X 的缺失數(shù)據(jù)與其他2 組相比明顯增多，導(dǎo)致與其他2 組相比估計的精度稍低.樣本量對估計的影響很小，所以估計比較穩(wěn)定.綜上所述，隨機模擬試驗的估計效果較好，估計值比較穩(wěn)定，并且精度也較高.

圖4 n=300，λ =12 時λ 的Gibbs 抽樣迭代過程Figure 4 Gibbs sampling iterations of λ with n=300，λ =12

圖5 n=300，λ =30 時λ 的MCEM 迭代過程Figure 5 MCEM iterations of λ with n=300，λ =30

圖6 n=300，λ =30 時λ 的Gibbs 抽樣迭代過程Figure 6 Gibbs sampling iterations of λ with n=300，λ =30

［1］Lu W，Shi D. A new compounding life distribution:The Weibull-Poisson distribution［J］. Journal of Applied Statistics，2012，39(1):21－38.

［2］Barreto-Souza W，Cribari-Neto F. A generalization of the exponential-Poisson distribution［J］. Statistics ＆ Probability Letters，2009，79(24):2493－2500.

［3］Poisson S D. Recherches sur la probabilité des jugements en matière criminelle et en matière civile，précédées des règles générales du calcul des probabilités［M］. Paris:Bachelier，1837.

［4］Consul P C，Jain G C. A generalization of the Poisson distribution［J］. Technometrics，1973，15(4):791－799.

［5］Cohen A C. Estimating the parameter in a conditional Poisson distribution［J］. Biometrics，1960，16(2):203－211.

［6］何朝兵，劉華文. 左截斷右刪失數(shù)據(jù)下二項分布參數(shù)多變點的貝葉斯估計［J］.華南師范大學(xué)學(xué)報:自然科學(xué)版，2014，46(3):34－38.He C B，Liu H W. Bayesian estimation of parameter of binomial distribution with multiple change points for left truncated and right censored data［J］. Journal of South China Normal University:Natural Science Edition，2014，46(3):34－38.

［7］Colchero F，Clark J S. Bayesian inference on age-specific survival for censored and truncated data［J］. Journal of Animal Ecology，2012，81(1):139－149.

［8］Balakrishnan N，Mitra D. Likelihood inference for lognormal data with left truncation and right censoring with an illustration［J］. Journal of Statistical Planning and Inference，2011，141(11):3536－3553.

［9］Cosslett S R. Efficient semiparametric estimation of censored and truncated regressions via a smoothed self-consistency equation［J］. Econometrica，2004，72 (4):1277－1293.

［10］Ahmadi J，Doostparast M，Parsian A. Estimation with left-truncated and right censored data:A comparison study［J］. Statistics ＆ Probability Letters，2012，82(7):1391－1400.

［11］Gross S T，Lai T L. Nonparametric estimation and regression analysis with left-truncated and right-censored data［J］. Journal of the American Statistical Association，1996，91(435):1166－1180.

［12］Dempster A P，Laird N M，Rubin D B. Maximum likelihood from incomplete data via the EM algorithm［J］.Journal of the Royal Statistical Society:Series B，1977，39:1－38.

［13］Geman S，Geman D. Stochastic relaxation，Gibbs distributions，and the Bayesian restoration of images［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1984 (6):721－741.