截断删失数据下泊松分布参数的点估计

何朝兵 ，柴士改，刘华文

(1.安阳师范学院数学与统计学院，安阳455000;2.山东大学数学学院，济南250100)

泊松分布是一种很重要的离散型分布，适于描述单位时间或空间内随机事件发生的次数，如某一服务设施在一定时间内到达的人数，电话交换机接到呼叫的次数等.文献［1］－［5］讨论了泊松分布的参数估计. 近年来，对截断删失数据的研究较多［6－11］，但缺乏涉及泊松分布的报道. 本文主要利用EM 算法和MCMC 方法对截断删失数据下泊松分布寿命参数的点估计进行了研究. 利用逆变换法和舍选法得到了产品的完全数据和参数的EM 迭代公式.把Gibbs 样本的算术平均值作为参数的MCMC 估计.随机模拟的估计效果较好，估计值比较稳定，且精度也较高.

1 截断删失数据模型

设(X，Y，T)是离散型随机变量，X 的分布函数和分布律分别为F(x，λ)=P(X≤x)、f(x，λ);Y 的分布函数和分布律分别为G(y)、g(y);T 的分布函数和分布律分别为H(t)、h(t)，假设X、Y、T 独立.对于n个受试样品，截断删失数据模型是:Zi≥Ti时可观察到数据(Zi，Ti，δi)，而Zi＜Ti时无法得到观察值，其中

下面求样本的似然函数.P(无样本观察值)=P(Zi＜Ti)=1－P(Xi≥Ti，Yi≥Ti)=

为方便研究，引入示性变量νi=I(min(Xi，Yi)≥Ti)(i=1，2，…，n).则基于数据{(zi，ti，δi):vi=1，1≤i≤n}的似然函数为

2 泊松分布参数的点估计

2.1 完全数据似然函数

截断删失数据模型中，假设产品寿命X 服从参数为λ 的泊松分布，即X ～P(λ).

由式(1)易知似然函数很复杂. 为了得到完全数据，不妨添加缺损寿命数据.具体如下:

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

νi=0 时，添加Z2i=Xi=z2i.Z2i的分布律为

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

所以可以利用舍选法随机产生Z2i的值z2i.

令δ，ν，z，u1，u2分别表示由δi，νi，zi，z1i，z2i组成的向量，则完全数据似然函数为

2.2 EM 算法

EM 算法［12］是一种求后验分布众数的迭代方法，它每一步迭代由E 步(求期望)和M 步(极大化)组成.下面求λ 的EM 迭代公式.

取λ 的先验分布为共轭先验分布伽玛分布Ga(α，β)，α，β 已知，即

则λ 的添加后验分布为

假设在第m+1 次迭代开始时的估计值为λ(m).令U1、U2分 别 表 示Z1i、Z2i组 成 的 向 量. 简 记

E 步:

当νi=1，δi=0 时，

当νi=0 时，利用Monte Carlo 方法计算Z2i的期望，称为MCEM方法，它将E 步变为如下2 步:

(E1)利用舍选法从ψ2(k，λ(m))抽取l个随机数k1，k2，…，kl;

(E2)计算1］，其中

M 步:

式(2)给出了由EM 算法得到的参数λ 的迭代公式.

2.3 贝叶斯估计

当νi=1，δi=0 时，ψ1(z1i，λ)，其中z－1i={z1j:j≠i}.

当νi=0 时λ)，其中z－2i={z2j:j≠i}.

参数λ 的满条件分布为

对z1i，z2i，λ 的满条件分布分别进行Gibbs 抽样［13］.设λ(j)(j=1，2，…，M1，…，M2)为λ 的容量为M2的Gibbs 样本，第M1次以后抽样收敛，把后面M2－M1个样本的均值作为λ 的估计，即

EM 估计是后验分布的众数(MLE)，MCMC 估计是后验分布的期望，它们差别不大.

3 随机模拟

设X ～P(λ)，Y ～P(λ1)，T ～P(λ2)，取λ 的先验分布为Ga(α，β)，对于超参数α、β 的确定应该充分利用各种先验信息，而利用先验矩确定超参数是一种常用的方法.虽然进行随机模拟试验时没有先验信息，选取α、β 时，却可以遵循“利用先验矩确定超参数”这一方法的思想，尽可能使先验期望α/β与E(X)= λ 相差不要太大，例如若λ =12，不妨取α=7，β =0.6，此时α/β≈11.67. 样本容量分别取n=30，50，100，200，300，500，800. 对于(λ，λ1，λ2)分别取3 组不同的参数进行随机模拟试验，具体如下:

(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. 参数估计结果见表1.

n=300 时的迭代过程见图1 ～图6.

表1 参数λ 的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)，截断删失模型的原理使得此组参数下的X 的缺失数据与其他2 组相比明显增多，导致与其他2 组相比估计的精度稍低.样本量对估计的影响很小，所以估计比较稳定.综上所述，随机模拟试验的估计效果较好，估计值比较稳定，并且精度也较高.

图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］何朝兵，刘华文. 左截断右删失数据下二项分布参数多变点的贝叶斯估计［J］.华南师范大学学报:自然科学版，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.