论文部分内容阅读
摘 要：对于设计寿命之前很少失效，设计寿命之后失效比例大幅增加的一类存储产品，其寿命变量用ZZ分布来描述较为合适。基于逐次定数截尾样本，首先利用泰勒展式将似然函数中的非线性部分转化为线性表达，使似然方程可解，进而得到参数的近似极大似然估计。其次，运用最佳线性无偏估计法，给出了参数无偏估计的一般解析表达式，同时针对两种特殊逐次定数截尾样本，进一步化简了估计的表达形式，且为便于使用，给出了计算所需数表的构造公式。最后对两种估计做了数值模拟比较，结果表明ZZ分布的最佳线性无偏估计优于近似极大似然估计。
关键词：ZZ分布;极大似然估计;近似极大似然估计;最佳线性无偏估计
DOI：10.15938/j.jhust.2021.03.023
中图分类号： O231
文献标志码： A
文章编号： 1007-2683（2021）03-0153-07
Parameter Estimation of ZZ Distribution under
Progressive Type-II Censored Data
ZHAO Hong-kai， ZHANG Guo-zhi， WANG Ping
（School of Science， Harbin University of Science and Technology，Harbin 150080，China）
Abstract：ZZ distribution is suitable to describe the life variable of a type of storage products that rarely fail before a given design life， and the proportion of failures has increased significantly after the design life. Based on a progressive type-II censored sample， firstly using Taylor′s formula to transform the the nonlinear part of the likelihood function into a linear expression， the likelihood equation is solvable， and the approximate maximum likelihood estimation（AMLE） of the parameter is given. Secondly， using the best linear unbiased estimation（BLUE） method， a general analytical expression of the parameter is given. At the same time， the expression form of the estimation is further simplified for two kinds of special progressive type-II censored sample and the required number table is constructed for the convenience of use. Finally， a numerical comparison of the two methods is performed， and the results show that the best linear unbiased estimation of ZZ distribution is better than the approximate maximum likelihood estimation.
Keywords：ZZ distribution; maximum likelihood estimation; approximate maximum likelihood estimation; best linear unbiased estimation
0 引 言
傳统的定时和定数截尾数据有一个共同的特点：在截尾时间点和截尾数上缺少灵活性，即在寿命试验中不允许有样品在试验中途退出试验。然而有些试验样品非常昂贵或者在研究耐久性试验时，如果试验人员为了减少试验时间希望在一些时间点将一些未失效的产品移出试验，传统的截尾样本数据就不适用了。随着基础理论的发展和实际应用中的需要，逐次定数截尾数据近年来受到广泛关注。
有关逐次定数截尾数据的研究，在国外，最早研究这类数据问题的文献可以追溯到1963年，文[1-3]详细的讨论了逐次截尾数据下针对常见几种分布的推断问题。1995年后，Sandhu等[4]给出了获得逐次定数截尾样本数据的抽样方法。Balasooryia等[5]和Tse等[6]分别在两参数指数分布和韦布尔分布下基于逐次定数截尾数据推导出了诸可靠性指标。Balakrishnan等[7]的著作对逐次定数截尾寿命试验进行了深入研究，给出了相关模型较为详细的阐述。Balakrishnan等[8-12]讨论了逐次定数截尾样本下正态分布、指数分布和极值分布参数的极大似然估计、区间估计和假设检验问题。Soliman[13]研究了服从Burr-XII分布的极大似然估计和Bayes估计。Mahmoud[14]推导出了逐次定数截尾样本下来自韦布尔-伽玛分布的顺序统计量的近似矩，并利用近似矩得到了参数的最佳线性无偏估计和极大似然估计。Singh等[15]研究了对数正态分布模型的统计分析问题，得到了参数的极大似然估计和Bayes估计。Mahmoud等[16]得到指数-帕累托参数的极大似然估计，并利用近似矩求得参数的最佳线性无偏估计。 国内对逐次定数截尾数据的研究始于上世纪九十年代，贺国芳等[17]说明了在实际问题中经常遇到这类数据，因此有必要讨论这类数据的统计分析方法。徐晓玲等[18]和王炳兴[19]讨论了Weibull分布基于逐次定数截尾寿命数据的参数估计，分别得到了参数的极大似然估计、逆矩估计以及区间估计。杨君慧等[20]给出了广义指数分布参数和可靠度函数的极大似然估计，并在熵损失和加权平方损失函数下，给出参数和可靠度函数的Bayes估计。杨敏[21]基于平方损失函数和LINEX损失函数，讨论了双参数Rayleigh分布参数的贝叶斯估计。罗嘉成等[22]探讨了Lomax分布的形状参数和可靠性指标的Bayes估计。
由于上述给出的研究均基于常见的几种寿命分布，然而有些存储产品，在给定的设计寿命之前很少失效，过了设计寿命之后失效的比例大幅增加，再用常见一些分布进行统计推断得到的结果与实际明显不符。张宁[23]给出一个较好解决上述存储产品寿命的ZZ分布，并对这类产品的可靠性指标在截尾数据情况下进行了相关的统计分析。
此外，上述对逐次定数截尾数据的研究，所获得的大多是参数极大似然估计和贝叶斯估计的数值解，并未给出具体的解析表达式，且研究多基于几种常见分布，并不能涵盖所有产品的寿命类型，因此本文在逐次定数截尾样本下，试给出ZZ分布参数估计的解析表达式，并配以数表以便使用。
1 ZZ分布与逐次定数截尾试验
1.1 ZZ分布
ZZ分布的分布函数F（t），密度函数f（t），失效率函数分别是λ（t）[23]：
F（t）=1-exp1-etηm，t>0
f（t）=mtm-1ηm·e（tη）m·e1-e（tη）m，t>0
λ（t）=mtm-1ηm·e（tη）m，t>0
它含有兩个参数η>0与m>0，记为ZZ（m，η），其中η称为特征寿命，m称为形状参数。
1.2 逐次定数截尾试验
产品实施逐次增加定数截尾试验的操作如下：从一批产品中随机抽取n个独立同分布产品在相同环境和应力条件下同时进行寿命试验，当观察到第1个样品失效时刻XR1∶r∶n时，在剩下的n-1个未失效的样品中任意抽取R1个撤离试验，还有n-1-R1个未失效的产品留下继续试验;当观察到第2个样品失效时刻XR2∶r∶n时，在剩下的n-2-R1个未失效样品中任意抽取R2个撤离试验，将余下n-2-R1-R2个未失效样品留下继续进行试验;如此下去，直至观察到r个样品失效时刻XRr∶r∶n时结束试验，此时余下的Rr=n-r-R1-R2-…-Rr-1个样品退出试验。则XR1∶r∶n<XR2∶r∶n<…<XRr∶r∶n可以看作是移出数据为R=[R1，R2，…，Rr]的逐次定数截尾试验数据，其中1<r≤n，且n=r+∑ri=1Ri。
考虑到在上述失效时刻XRj∶r∶n，j=1，2，…，r，都会随机移出Rj个元件退出试验，故上述观察到的失效时间不再是容量为n的样本的前r个次序观测值，而是n个次序统计量中的某r个顺序观测值。实际上，可以将XR1∶r∶n看成容量为n的寿命为t1<t2<…<tn的第i1（i1=1）个次序统计量ti1，将XR2∶r∶n看成容量为n的寿命为t1<t2<…<tn的第i2个次序统计量ti2，直到XRr∶r∶n为容量为n的寿命为t1<t2<…<tn的第ir个次序统计量tir，其中1=i1≤i2≤…≤ir≤n。
设X=[XR1∶r∶n，XR2∶r∶n，…，XRr∶r∶n]是来自ZZ分布逐次定数截尾的样本，R=[R1，R2，…，Rr]是试验中相应被移出的产品数，为书写方便，简记Xj≡XRj∶r∶n，j=1，2，…，r，则逐次定数截尾下ZZ分布的极大似然函数为
L=c∏rj=1f（xj）[1-F（xj）]Rj=
cmrηmr∏rj=1xm-1j·exp∑rj=1（Rj+1）1-exjηm+∑rj=1xmjηm（1）
其中
c=n（n-1-R1）（n-2-R1-R2）…（n-r-R1-R2-…-Rr+1）是与未知参数无关的常数因子。
2 近似极大似然估计
针对ZZ分布在逐次定数截尾样本下的极大似然估计没有解析解的情况，借鉴文[9]的方法，考虑近似似然函数，并求解显示结果。
设随机变量X～ZZ（m，η），则Y=lnX服从EZ分布[23]并且它的概率密度函数为
fy（μ，σ）=1σey-μσ·expey-μσ·exp1-expey-μσ
其中μ=lnη，σ=1/m。
μ=0，σ=1时标准的EZ分布的密度函数和分布函数有以下的形式：
fy（0，1）=ey·exp{ey}exp{1-exp{ey}}
F（y）=1-exp{1-exp{ey}}
X=[XR1∶r∶n，XR2∶r∶n，…，XRr∶r∶n]是来自ZZ分布逐次定数截尾的随机样本，令YRj∶r∶n=lnXRj∶r∶n，zRj∶r∶n=YRj∶r∶n-μσ，则z=[zR1∶r∶n，zR2∶r∶n，…，zRr∶r∶n]可看成来自标准EZ分布逐次定数截尾的随机样本，简记z≡[z1，z2，…，zr]，那么由式（1）得到的似然函数可表示成以下式子：
L（μ，σ）=C1σr∏rj=1w（zj）（W（zj））Rj
其中：w（z）=exp{z+ez+[1-exp（ez）]}，W（z）=exp{1-exp（ez）}。
从而，对数似然函数为
lnL=-rlnσ+∑rj=1ln（w（zj））+∑rj=1Rj·ln（W（zj）） 上式對μ，σ分别求偏导数，有
lnLμ=-1σ∑rj=1w′（zj）w（zj）-∑rj=1RjσW′（zj）W（zj）=
-1σ∑rj=1w′（zj）w（zj）+∑rj=1Rjσw（zj）W（zj）
lnLσ=-mσ-∑rj=1zjσw′（zj）w（zj）+∑rj=1Rjzjσw（zj）W（zj）
得到下列的似然方程组：
-∑rj=1w′（zj）w（zj）+∑rj=1Rjw（zj）W（zj）=0（2）
-m-∑rj=1zjw′（zj）w（zj）+∑rj=1Rjzjw（zj）W（zj）=0（3）
由式（2）和式（3）组成的似然方程组不能解出和的明显表达式。由于式（2）和式（3）不能直接解出是w′（zj）/w（zj）和w（zj）/W（zj）引起的，根据文[9]中给出的方法，利用泰勒展式在点cj=E（zj）处将w′（zj）/w（zj）和w（zj）/W（zj）近似展开。其中，zj=F-1（Uj），Uj，j=1，2，…，r是来自于总体为U（0，1）均匀分布的逐次定数截尾样本。由于cj=E（zj）直接求解比较困难，而E（Uj）在文[4]中已有结论，故有如下的近似式
cj=E（zj）≈F-1（E（Uj））
其中F-1（Uj）=lnln（1-ln（1-Uj）），且
pj=E（Uj）=1-∏rl=r-j+1l+Rr+Rr-1+…Rr-l+11+l+Rr+Rr-1+…Rr-l+1，qj=1-pj，j=1，2，…，r。如此，能得到如下的展开式：
w′（zj）w（zj）≈aj-bjzj
w（zj）W（zj）≈1-aj+bjzj
其中，
aj=w′（cj）w（cj）-cjw″（cj）w（cj）-w′（cj）w（cj）2=
1+ecj[1-cj-（1+cj）exp（ecj）-cjexp（ecj）ecj]=
1+ln（1-lnqj）[1-lnln（1-lnqj）-
（1+lnln（1-lnqj））（1-lnqj）-
lnln（1-lnqj）ln（1-lnqj）（1-lnqj）]
bj=w′（cj）w（cj）2-w″（cj）w（cj）=
ecj[exp（ecj）·ecj+exp（ecj）-1]=
ln（1-lnqj）[（1-lnqj）·ln（1-lnqj）-lnqj]
那么极大似然方程组式（2）和式（3）能近似表示如下
-∑rj=1（aj-bjzj）+∑rj=1Rj（1-aj+bjzj）=0（4）
-r-∑rj=1（aj-bjzj）zj+∑rj=1Rj（1-aj+bjzj）zj=0（5）
将zj=yj-μσ代入式（4），化简得：
=∑rj=1（1+Rj）bjyj∑rj=1（1+Rj）bj-∑rj=1（1+Rj）aj-∑rj=1Rj∑rj=1（1+Rj）bj（6）
记A=∑rj=1（1+Rj）bjyj∑rj=1（1+Rj）bj，B=∑rj=1（1+Rj）aj-∑rj=1Rj∑rj=1（1+Rj）bj，
则式（6）可以表示为
=A-B（7）
然后，再将zj=yj-μσ及式（7）代入式（5），得到关于的二次方程：
D2+E-H=0（8）
其中，
D=r+B∑rj=1aj-B∑rj=1Rj（1-aj）-B2∑rj=1bj-B2∑rj=1Rjbj=
r+B∑rj=1（1+Rj）aj-∑rj=1Rj-∑rj=1（1+Rj）aj-∑rj=1Rj=r
E=∑rj=1（aj+Rjaj-Rj）（yj-A）-2B∑rj=1（1+Rj）bj（yj-A）
H=∑rj=1（1+Rj）bj（yj-A）2>0
解方程式（8），得到关于的值：
=-E+E2+4rH2r
将求出的的值代入式（7）即可求出的值。
3 最佳线性无偏估计
设在逐次截尾R=[R1，R2，…，Rr]的试验中，失效元件的寿命记为ti1≤ti2≤…≤tir，其下标看做是随机向量（I1，I2，…，Ir）的一组观测值，该随机向量所有可能取值的集合记为G，设其分布律为P{I1=i1，I2=i2，…，Ir=ir}≡Pi1，i2，…，ir，其中[i1，i2，…，ir]∈G且1=i1≤i2≤…≤ir≤n。则在{I1=i1，I2=i2，…，Ir=ir}条件下，有
F（tij）=1-exp1-etijηm，j=1，2，…，r
移项，并取三重对数后可得
lnln（1-ln（1-F（tij）））=mlntij-mlnη（9）
令Ti1，i2，…，ir=[Ti1，Ti2，…，Tir]T，其中Tij=lnln（1-ln（1-F（tij））），j=1，2，…，r。显然F（ti1）<F（ti2）<…<F（tir）是来自U（0，1）均匀分布的统计量，设
ETi1，i2，…，ir=ai1，i2，…，ir=[ai1，ai2，…，air]T，Σi1，i2，…，ir=E[Ti1，i2，…，ir-ai1，i2，…，ir][Ti1，i2，…，ir-ai1，i2，…，ir]T则Ti1，i2，…，ir均值和方差只与n和r有关，而不依赖于其它参数。那么式（9）可表示为
lntij=lnη+1maij+1mδij，j=1，2，…，r（10）
其中δi1，i2，…，ir=[δi1，δi2，…，δir]T，并且 E（δi1，i2，…，ir）=0，Cov（δi1，i2，…，ir）=Σi1，i2，…，ir。
再令
Yi1，i2，…，ir=[Yi1，Yi2，…，Yir]T=[lnti1，lnti2，…，lntir]T
Xi1，i2，…，ir=1ai11ai21ar
θ=（μ，σ）T，其中μ=lnη，σ=1m
εi1，i2，…，ir=1mδi1，i2，…，ir
将式（10）转化为如下矩阵形式
Yi1，i2，…，ir=Xi1，i2，…，irθ+εi1，i2，…，ir，
Cov（εi1，i2，…，ir）=σ2Σi1，i2，…，ir。
根据高斯-马尔科夫定理，用加权最小二乘法求得参数的最佳线性无偏估计（BLUE）
i1，i2，…，ir=i1，i2，…，ir
i1，i2，…，ir=
[X′i1，i2，…，ir-1i1，i2，…，irXi1，i2，…，ir]-1
X′i1，i2，…，ir×-1i1，i2，…，irYi1，i2，…，ir
当n，r及i1，i2，…，ir给定时，矩阵Xi1，i2，…，ir与Σi1，i2，…，ir都是已知的，所以为了便于使用，将参数的最佳线性无偏估计简化如下形式
（i1，i2，…，ir）=∑rj=1D0（i1，i2，…，ir，n，r，j）lntij
（i1，i2，…，ir）=∑rj=1D1（i1，i2，…，ir，n，r，j）lntij
其中：D0（i1，i2，…，ir，n，r，j）称为μ（i1，i2，…，ir）的最佳线性无偏估计系数;D1（i1，i2，…，ir，n，r，j）称为σ（i1，i2，…，ir）的最佳线性无偏估计系数。
因此，在R=（R1，R2，…，Rr）下，ZZ分布参数最佳线性无偏估计的一般表达形式为
=∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤nPi1，i2，…，ir·（i1，i2，…，ir）（11）
=∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤nPi1，i2，…，ir·（i1，i2，…，ir）（12）
上面的表达式涉及到系数D0（i1，i2，…，ir，n，r，j），D1（i1，i2，…，ir，n，r，j）和概率P{I1=i1，I2=i2，…，Ir=ir}≡Pi1，i2，…，ir的计算，本文仅对下面两种情况，给出参数估计的简化形式，并构造相应的数表计算公式，便于使用。
3.1 R=[n-r，0，0，…，0]1×r时的参数估计
当移出数据为R=[n-r，0，0，…，0]1×r时，n个元件参加试验，在第一个元件失效时刻一次随机移出n-r个未失效元件，剩下的r-1个元件继续参加试验直至全部失效。ti1<ti2<…<tir可以看成来自容量为n的次序统计量t1<t2<…<tn中固定第一个次序统计量后，从余下n-1元件中随机抽取r-1个的子样本，且该子样本有Cr-1n-1种等可能。所以此时有
P{I1=i1，I2=i2，…，Ir=ir}≡Pi1，i2，…，ir=1/Cr-1n-1
由式（11）和式（12）得ZZ分布参数估计可表示为如下形式
R=（n-r，0，0，…，0）=1Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤n（i1，i2，…，ir）
R=（n-r，0，0，…，0）=1Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤n（i1，i2，…，ir）
令
1Cr-1n-1∑（i1，i2，…，ir）∈G1=ii<i2<…<ir≤nD0（ii，i2，…，ir，n，r，j）=D′0，R=（n-r，0，0，…，0）（n，r，j）
1Cr-1n-1∑（i1，i2，…，ir）∈G1=ii<i2<…<ir≤nD1（ii，i2，…，ir，n，r，j）=D′1，R=（n-r，0，0，…，0）（n，r，j）
则R=[n-r，0，0，…，0]1×r时ZZ分布的最佳线性无偏估计如下
R=（n-r，0，0，…，0）=1Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤n∑rj=1D0（i1，i2，…，ir，n，r，j）lntij=
∑rj=11Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤nD0（i1，i2，…，ir，n，r，j）lntij=
∑rj=1D′0，R=（n-r，0，0，…，0）（n，r，j）lntij
R=（n-r，0，0，…，0）=1Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤n∑rj=1D1（i1，i2，…，ir，n，r，j）lntij=
∑rj=11Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤nD1（i1，i2，…，ir，n，r，j）lntij=
∑rj=1D′1，R=（n-r，0，0，…，0）（n，r，j）lntij
為了便于计算参数的估计，其系数D′0，R=（n-r，0，0，…，0）（n，r，j）和D′1，R=（n-r，0，0，…，0）（n，r，j）在n=5到n=20时已由计算机给出。
3.2 R=[1，1，…，1，0，0，…，0]1×r时的参数估计
在每次失效时刻XRj∶r∶n，j=1，2，…，r分别移出1个未失效产品的情况：n个元件参加试验，设有r个失效，且在每次失效时刻移出1个未失效元件。在R=[1，1，…1，0，0，…，0]1×r时，每种情况可能出现的概率P{I1=i1，I2=i2，…，Ir=ir}≡Pi1，i2，…，ir不再相等，但其计算可借助计算机完成，最终体现在后面的数表中。 由式（11）和式（12），移出數据为R=[1，1，…，1，0，…，0]1×r的ZZ分布的最佳线性无偏估计表示如下
R=（1，…，1，0，…，0）=∑rj=1D′0，R=（1，…，1，0，…，0）（n，r，j）lntij
R=（1，…，1，0，…，0）=∑rj=1D′1，R=（1，…，1，0，…，0）（n，r，j）lntij
其中，
D′0，R=（1，…，1，0，…，0）（n，r，j）=∑（i1，i2，…，ir）∈G1=ii<i2<…<ir≤nPi1，i2，…，ir×
D0（ii，i2，…，ir，n，r，j）
D′1，R=（1，…，1，0，…，0）（n，r，j）=∑（i1，i2，…，ir）∈G1=ii<i2<…<ir≤nPi1，i2，…，ir×
D1（ii，i2，…，ir，n，r，j）
为了便于计算参数的估计，其系数D′0，R=（1，1，…，1，0，…，0）（n，r，j）和D′1，R=（1，1，…，1，0，…，0）（n，r，j）在n=5到n=20时已由计算机给出。
4 数值模拟
基于Monte Carlo模拟，给出以上两种方法下ZZ分布参数的估计结果。根据文[4]给出的算法，产生ZZ分布下逐次定数截尾数据的步骤如下：
1）产生r个来自均匀分布U（0，1）的独立随机样本W1，W2，…，Wr;
2）在预先设定的逐次定数截尾移出数据
R=[R1，R2，…，Rr]下，令
Vj=W1/（j+Rr+Rr-1+…Rr-j+1）j，j=1，2，…，r;
3）再令Uj=1-VrVr-1…Vm-i+1，j=1，2，…，r，则U1，U2，…，Ur是来自均匀分布U（0，1）的逐次定数截尾数据;
4）令Xj=F-1（Uj），j=1，2，…，r，得到来自ZZ分布的逐次定数截尾试验数据，其中F-1（·）是ZZ分布函数的反函数，
Xj=（β·ln（1-ln（1-Uj）））1/m，j=1，2，…，r。
设定两组不同样本量和两组不同的移出数据，根据3.1和3.2节方法，构造最佳线性无偏估计的系数值分别如表1和表2所示。
应用上述数据，模拟得参数估计的结果如表3所示。
从表3中可以看出，最佳线性无偏估计结果和近似极大似然估计结果都较接近参数真值，且最佳线性无偏估计结果较优。因此，二者均可应用于寿命服从ZZ分布、逐次定数截尾寿命试验下的参数估计。
5 结 论
本文基于逐次定数截尾寿命试验模型，运用近似极大似然法和最佳线性无偏估计法，对产品寿命服从ZZ分布的参数进行了统计分析。主要构建了r维不可观测但可求分布的随机向量，据此将收集到的逐次定数截尾样本数据转化为一定概率分布条件下的顺序统计量，进而给出ZZ分布参数的最佳线性无偏估计，并可依据构造出的数表再进行简单计算即得到参数的估计值。且此法具有的解析表达形式、运用数表计算方便等优点，易可应用到Weibull分布或其他指数型寿命分布中。
参 考 文 献：
[1] COHEN， CLIFFORD A. Progressively Censored Samples in Life Testing[J].Technometrics， 1963， 5（3）：327.
[2] COHEN A C. Progressively Censored Sampling in the Three Parameter Log-normal Distribution[J]. Technometrics， 1976， 18（1）：99.
[3] NORGAARD C N J. Progressively Censored Sampling in the Three-parameter Gamma Distribution[J]. Technometrics， 1977， 19（3）：333.
[4] SANDHU， BALAKRISHNANR N. A Simple Simulational Algorithm for Generating Progressive Type-II Censored Samples[J]. The American Statistician， 1995， 49（2）：229.
[5] BALASOORIYA U， SAW S L C. Reliability Sampling Plans for the Two-parameter Exponential Distribution under Progressive Censoring[J]. Journal of Applied Statistics， 1998， 25（5）：707.
[6] TSE S K， YANG C. Reliability Sampling Plans for the Weibull Distribution under Type-II Progressive Censoring with Binomial Removals[J]. Journal of Applied Statistics， 2003， 30（6）：709.
[7] BALAKRISHNAN N， AGGARWALA R. Progressive Censoring[M]. Birkhauser， 2000.
[8] BALAKRISHNAN N， JIE Mi. Existence and Uniqueness of the MLEs for Normal Distribution Based on General Progressively Type-II Censored Samples[J]. Statistics and Probability Letters， 2003， 64（4）：407. [9] BALAKRISHNAN N， KANNAN N， LIN C T. Point and Interval Estimation for Gaussian Distribution， Based on Progressively Type-II Censored Samples[J]. IEEE Transactions on Reliability， 2003， 52（1）：90.
[10]BALAKRISHNAN N， LIN C T. On the Distribution of a Test for Exponentiality Based on Progressively Type-II Right Censored Spacings[J]. Journal of Statistical Computation & Simulation， 2003， 73（4）：277.
[11]BALAKRISHNAN N， NG H K T， KANNAN N. A Test of Exponentiality Based on Spacings for Progressively Type-II Censored Data[J]. Goodness-of-it Tests and Model Validity， 2002：89.
[12]BALAKRISHNAN N， KANNAN N， LIN C T， et al. Inference for the Extreme Value Distribution under Progressive Type-II Censoring[J]. Journal of Statistical Computation and Simulation， 2004， 74（1）：25.
[13]SOLIMAN， A. A. Estimation of Parameters of Life from Progressively Censored Data Using Burr-XII Model[J]. IEEE Transactions on Reliability， 2005， 54（1）：34.
[14]MAHMOUD M A W， MOSHREF M， YHIEA N M， et al. Progressively Censored Data from the Weibull Gamma Distribution Moments and Estimations[J]. J Appl Probab Stat， 2014， 3 （1）：45.
[15]SINGH S， TRIPATHI Y M， WU S J. On Estimating Parameters of a Progressively Censored Log-normal Distribution[J]. Journal of Statistical Computation and Simulation， 2015， 85（6）：1071.
[16]ALIM A M， JAHEEN Z F. Estimation of Parameters for the Exponentiated Pareto Distribution Based on Progressively Type-II Right Censored Data[J]. Journal of the Egyptian Mathematical Society， 2016， 24（3）.
[17]賀国芳， 瞿荣贞编写. 可靠性数据的收集与分析[M]. 北京：国防工业出版社， 1995.41.
[18]徐晓玲， 王蓉华. Weibull 分布逐步增加的Ⅱ型截尾试验的统计分析[J].强度与环境， 2003， 30（2）：31.
XU Xiaoling， WANG Ronghua. Statistical Analysis of Type Ⅱ Censoring Test with Increasing Weibull Distribution[J]. Strength and Environment， 2003， 30 （2）： 31.
[19]王炳兴. Weibull分布基于定数逐次截尾寿命数据的统计分析[J].科技通报， 2004（6）：488.
WANG Bingxing. Statistical Analysis of Weibull Distribution Based on Constant Successive Censored Life Data[J]. Science and Technology Bulletin， 2004 （6）： 488.
[20]杨君慧， 师义民， 曹弘毅. 逐步增加Ⅱ型截尾试验下广义指数分布的统计分析[J]. 统计与决策， 2014（16）：28.
YANG Junhui， SHI Yimin， CAO Hongyi. Statistical Analysis of Generalized Exponential Distribution under Gradually Increasing Type II Censoring Test[J]. Statistics and Decision Making， 2014 （16）： 28.
[21]杨敏. 删失数据场合双参数Rayleigh分布的参数估计[D]. 广西师范学院， 2015：12.
[22]罗嘉成， 陈勇明. 逐次定数截尾下Lomax分布的贝叶斯分析[J]. 成都信息工程大学学报， 2016， 31（3）：323.
LUO Jiacheng， CHEN Yongming. Bayesian Analysis for Lomax Distribution under Progressively Type-Ⅱ Censored Data[J]. Journal of Chengdu University of Information Technology， 2016，31（3）：323.
[23]张宁. 截尾数据下ZZ分布参数估计及性质[D]. 哈尔滨：哈尔滨理工大学， 2019.
（编辑：温泽宇）
关键词：ZZ分布;极大似然估计;近似极大似然估计;最佳线性无偏估计
DOI：10.15938/j.jhust.2021.03.023
中图分类号： O231
文献标志码： A
文章编号： 1007-2683（2021）03-0153-07
Parameter Estimation of ZZ Distribution under
Progressive Type-II Censored Data
ZHAO Hong-kai， ZHANG Guo-zhi， WANG Ping
（School of Science， Harbin University of Science and Technology，Harbin 150080，China）
Abstract：ZZ distribution is suitable to describe the life variable of a type of storage products that rarely fail before a given design life， and the proportion of failures has increased significantly after the design life. Based on a progressive type-II censored sample， firstly using Taylor′s formula to transform the the nonlinear part of the likelihood function into a linear expression， the likelihood equation is solvable， and the approximate maximum likelihood estimation（AMLE） of the parameter is given. Secondly， using the best linear unbiased estimation（BLUE） method， a general analytical expression of the parameter is given. At the same time， the expression form of the estimation is further simplified for two kinds of special progressive type-II censored sample and the required number table is constructed for the convenience of use. Finally， a numerical comparison of the two methods is performed， and the results show that the best linear unbiased estimation of ZZ distribution is better than the approximate maximum likelihood estimation.
Keywords：ZZ distribution; maximum likelihood estimation; approximate maximum likelihood estimation; best linear unbiased estimation
0 引 言
傳统的定时和定数截尾数据有一个共同的特点：在截尾时间点和截尾数上缺少灵活性，即在寿命试验中不允许有样品在试验中途退出试验。然而有些试验样品非常昂贵或者在研究耐久性试验时，如果试验人员为了减少试验时间希望在一些时间点将一些未失效的产品移出试验，传统的截尾样本数据就不适用了。随着基础理论的发展和实际应用中的需要，逐次定数截尾数据近年来受到广泛关注。
有关逐次定数截尾数据的研究，在国外，最早研究这类数据问题的文献可以追溯到1963年，文[1-3]详细的讨论了逐次截尾数据下针对常见几种分布的推断问题。1995年后，Sandhu等[4]给出了获得逐次定数截尾样本数据的抽样方法。Balasooryia等[5]和Tse等[6]分别在两参数指数分布和韦布尔分布下基于逐次定数截尾数据推导出了诸可靠性指标。Balakrishnan等[7]的著作对逐次定数截尾寿命试验进行了深入研究，给出了相关模型较为详细的阐述。Balakrishnan等[8-12]讨论了逐次定数截尾样本下正态分布、指数分布和极值分布参数的极大似然估计、区间估计和假设检验问题。Soliman[13]研究了服从Burr-XII分布的极大似然估计和Bayes估计。Mahmoud[14]推导出了逐次定数截尾样本下来自韦布尔-伽玛分布的顺序统计量的近似矩，并利用近似矩得到了参数的最佳线性无偏估计和极大似然估计。Singh等[15]研究了对数正态分布模型的统计分析问题，得到了参数的极大似然估计和Bayes估计。Mahmoud等[16]得到指数-帕累托参数的极大似然估计，并利用近似矩求得参数的最佳线性无偏估计。 国内对逐次定数截尾数据的研究始于上世纪九十年代，贺国芳等[17]说明了在实际问题中经常遇到这类数据，因此有必要讨论这类数据的统计分析方法。徐晓玲等[18]和王炳兴[19]讨论了Weibull分布基于逐次定数截尾寿命数据的参数估计，分别得到了参数的极大似然估计、逆矩估计以及区间估计。杨君慧等[20]给出了广义指数分布参数和可靠度函数的极大似然估计，并在熵损失和加权平方损失函数下，给出参数和可靠度函数的Bayes估计。杨敏[21]基于平方损失函数和LINEX损失函数，讨论了双参数Rayleigh分布参数的贝叶斯估计。罗嘉成等[22]探讨了Lomax分布的形状参数和可靠性指标的Bayes估计。
由于上述给出的研究均基于常见的几种寿命分布，然而有些存储产品，在给定的设计寿命之前很少失效，过了设计寿命之后失效的比例大幅增加，再用常见一些分布进行统计推断得到的结果与实际明显不符。张宁[23]给出一个较好解决上述存储产品寿命的ZZ分布，并对这类产品的可靠性指标在截尾数据情况下进行了相关的统计分析。
此外，上述对逐次定数截尾数据的研究，所获得的大多是参数极大似然估计和贝叶斯估计的数值解，并未给出具体的解析表达式，且研究多基于几种常见分布，并不能涵盖所有产品的寿命类型，因此本文在逐次定数截尾样本下，试给出ZZ分布参数估计的解析表达式，并配以数表以便使用。
1 ZZ分布与逐次定数截尾试验
1.1 ZZ分布
ZZ分布的分布函数F（t），密度函数f（t），失效率函数分别是λ（t）[23]：
F（t）=1-exp1-etηm，t>0
f（t）=mtm-1ηm·e（tη）m·e1-e（tη）m，t>0
λ（t）=mtm-1ηm·e（tη）m，t>0
它含有兩个参数η>0与m>0，记为ZZ（m，η），其中η称为特征寿命，m称为形状参数。
1.2 逐次定数截尾试验
产品实施逐次增加定数截尾试验的操作如下：从一批产品中随机抽取n个独立同分布产品在相同环境和应力条件下同时进行寿命试验，当观察到第1个样品失效时刻XR1∶r∶n时，在剩下的n-1个未失效的样品中任意抽取R1个撤离试验，还有n-1-R1个未失效的产品留下继续试验;当观察到第2个样品失效时刻XR2∶r∶n时，在剩下的n-2-R1个未失效样品中任意抽取R2个撤离试验，将余下n-2-R1-R2个未失效样品留下继续进行试验;如此下去，直至观察到r个样品失效时刻XRr∶r∶n时结束试验，此时余下的Rr=n-r-R1-R2-…-Rr-1个样品退出试验。则XR1∶r∶n<XR2∶r∶n<…<XRr∶r∶n可以看作是移出数据为R=[R1，R2，…，Rr]的逐次定数截尾试验数据，其中1<r≤n，且n=r+∑ri=1Ri。
考虑到在上述失效时刻XRj∶r∶n，j=1，2，…，r，都会随机移出Rj个元件退出试验，故上述观察到的失效时间不再是容量为n的样本的前r个次序观测值，而是n个次序统计量中的某r个顺序观测值。实际上，可以将XR1∶r∶n看成容量为n的寿命为t1<t2<…<tn的第i1（i1=1）个次序统计量ti1，将XR2∶r∶n看成容量为n的寿命为t1<t2<…<tn的第i2个次序统计量ti2，直到XRr∶r∶n为容量为n的寿命为t1<t2<…<tn的第ir个次序统计量tir，其中1=i1≤i2≤…≤ir≤n。
设X=[XR1∶r∶n，XR2∶r∶n，…，XRr∶r∶n]是来自ZZ分布逐次定数截尾的样本，R=[R1，R2，…，Rr]是试验中相应被移出的产品数，为书写方便，简记Xj≡XRj∶r∶n，j=1，2，…，r，则逐次定数截尾下ZZ分布的极大似然函数为
L=c∏rj=1f（xj）[1-F（xj）]Rj=
cmrηmr∏rj=1xm-1j·exp∑rj=1（Rj+1）1-exjηm+∑rj=1xmjηm（1）
其中
c=n（n-1-R1）（n-2-R1-R2）…（n-r-R1-R2-…-Rr+1）是与未知参数无关的常数因子。
2 近似极大似然估计
针对ZZ分布在逐次定数截尾样本下的极大似然估计没有解析解的情况，借鉴文[9]的方法，考虑近似似然函数，并求解显示结果。
设随机变量X～ZZ（m，η），则Y=lnX服从EZ分布[23]并且它的概率密度函数为
fy（μ，σ）=1σey-μσ·expey-μσ·exp1-expey-μσ
其中μ=lnη，σ=1/m。
μ=0，σ=1时标准的EZ分布的密度函数和分布函数有以下的形式：
fy（0，1）=ey·exp{ey}exp{1-exp{ey}}
F（y）=1-exp{1-exp{ey}}
X=[XR1∶r∶n，XR2∶r∶n，…，XRr∶r∶n]是来自ZZ分布逐次定数截尾的随机样本，令YRj∶r∶n=lnXRj∶r∶n，zRj∶r∶n=YRj∶r∶n-μσ，则z=[zR1∶r∶n，zR2∶r∶n，…，zRr∶r∶n]可看成来自标准EZ分布逐次定数截尾的随机样本，简记z≡[z1，z2，…，zr]，那么由式（1）得到的似然函数可表示成以下式子：
L（μ，σ）=C1σr∏rj=1w（zj）（W（zj））Rj
其中：w（z）=exp{z+ez+[1-exp（ez）]}，W（z）=exp{1-exp（ez）}。
从而，对数似然函数为
lnL=-rlnσ+∑rj=1ln（w（zj））+∑rj=1Rj·ln（W（zj）） 上式對μ，σ分别求偏导数，有
lnLμ=-1σ∑rj=1w′（zj）w（zj）-∑rj=1RjσW′（zj）W（zj）=
-1σ∑rj=1w′（zj）w（zj）+∑rj=1Rjσw（zj）W（zj）
lnLσ=-mσ-∑rj=1zjσw′（zj）w（zj）+∑rj=1Rjzjσw（zj）W（zj）
得到下列的似然方程组：
-∑rj=1w′（zj）w（zj）+∑rj=1Rjw（zj）W（zj）=0（2）
-m-∑rj=1zjw′（zj）w（zj）+∑rj=1Rjzjw（zj）W（zj）=0（3）
由式（2）和式（3）组成的似然方程组不能解出和的明显表达式。由于式（2）和式（3）不能直接解出是w′（zj）/w（zj）和w（zj）/W（zj）引起的，根据文[9]中给出的方法，利用泰勒展式在点cj=E（zj）处将w′（zj）/w（zj）和w（zj）/W（zj）近似展开。其中，zj=F-1（Uj），Uj，j=1，2，…，r是来自于总体为U（0，1）均匀分布的逐次定数截尾样本。由于cj=E（zj）直接求解比较困难，而E（Uj）在文[4]中已有结论，故有如下的近似式
cj=E（zj）≈F-1（E（Uj））
其中F-1（Uj）=lnln（1-ln（1-Uj）），且
pj=E（Uj）=1-∏rl=r-j+1l+Rr+Rr-1+…Rr-l+11+l+Rr+Rr-1+…Rr-l+1，qj=1-pj，j=1，2，…，r。如此，能得到如下的展开式：
w′（zj）w（zj）≈aj-bjzj
w（zj）W（zj）≈1-aj+bjzj
其中，
aj=w′（cj）w（cj）-cjw″（cj）w（cj）-w′（cj）w（cj）2=
1+ecj[1-cj-（1+cj）exp（ecj）-cjexp（ecj）ecj]=
1+ln（1-lnqj）[1-lnln（1-lnqj）-
（1+lnln（1-lnqj））（1-lnqj）-
lnln（1-lnqj）ln（1-lnqj）（1-lnqj）]
bj=w′（cj）w（cj）2-w″（cj）w（cj）=
ecj[exp（ecj）·ecj+exp（ecj）-1]=
ln（1-lnqj）[（1-lnqj）·ln（1-lnqj）-lnqj]
那么极大似然方程组式（2）和式（3）能近似表示如下
-∑rj=1（aj-bjzj）+∑rj=1Rj（1-aj+bjzj）=0（4）
-r-∑rj=1（aj-bjzj）zj+∑rj=1Rj（1-aj+bjzj）zj=0（5）
将zj=yj-μσ代入式（4），化简得：
=∑rj=1（1+Rj）bjyj∑rj=1（1+Rj）bj-∑rj=1（1+Rj）aj-∑rj=1Rj∑rj=1（1+Rj）bj（6）
记A=∑rj=1（1+Rj）bjyj∑rj=1（1+Rj）bj，B=∑rj=1（1+Rj）aj-∑rj=1Rj∑rj=1（1+Rj）bj，
则式（6）可以表示为
=A-B（7）
然后，再将zj=yj-μσ及式（7）代入式（5），得到关于的二次方程：
D2+E-H=0（8）
其中，
D=r+B∑rj=1aj-B∑rj=1Rj（1-aj）-B2∑rj=1bj-B2∑rj=1Rjbj=
r+B∑rj=1（1+Rj）aj-∑rj=1Rj-∑rj=1（1+Rj）aj-∑rj=1Rj=r
E=∑rj=1（aj+Rjaj-Rj）（yj-A）-2B∑rj=1（1+Rj）bj（yj-A）
H=∑rj=1（1+Rj）bj（yj-A）2>0
解方程式（8），得到关于的值：
=-E+E2+4rH2r
将求出的的值代入式（7）即可求出的值。
3 最佳线性无偏估计
设在逐次截尾R=[R1，R2，…，Rr]的试验中，失效元件的寿命记为ti1≤ti2≤…≤tir，其下标看做是随机向量（I1，I2，…，Ir）的一组观测值，该随机向量所有可能取值的集合记为G，设其分布律为P{I1=i1，I2=i2，…，Ir=ir}≡Pi1，i2，…，ir，其中[i1，i2，…，ir]∈G且1=i1≤i2≤…≤ir≤n。则在{I1=i1，I2=i2，…，Ir=ir}条件下，有
F（tij）=1-exp1-etijηm，j=1，2，…，r
移项，并取三重对数后可得
lnln（1-ln（1-F（tij）））=mlntij-mlnη（9）
令Ti1，i2，…，ir=[Ti1，Ti2，…，Tir]T，其中Tij=lnln（1-ln（1-F（tij））），j=1，2，…，r。显然F（ti1）<F（ti2）<…<F（tir）是来自U（0，1）均匀分布的统计量，设
ETi1，i2，…，ir=ai1，i2，…，ir=[ai1，ai2，…，air]T，Σi1，i2，…，ir=E[Ti1，i2，…，ir-ai1，i2，…，ir][Ti1，i2，…，ir-ai1，i2，…，ir]T则Ti1，i2，…，ir均值和方差只与n和r有关，而不依赖于其它参数。那么式（9）可表示为
lntij=lnη+1maij+1mδij，j=1，2，…，r（10）
其中δi1，i2，…，ir=[δi1，δi2，…，δir]T，并且 E（δi1，i2，…，ir）=0，Cov（δi1，i2，…，ir）=Σi1，i2，…，ir。
再令
Yi1，i2，…，ir=[Yi1，Yi2，…，Yir]T=[lnti1，lnti2，…，lntir]T
Xi1，i2，…，ir=1ai11ai21ar
θ=（μ，σ）T，其中μ=lnη，σ=1m
εi1，i2，…，ir=1mδi1，i2，…，ir
将式（10）转化为如下矩阵形式
Yi1，i2，…，ir=Xi1，i2，…，irθ+εi1，i2，…，ir，
Cov（εi1，i2，…，ir）=σ2Σi1，i2，…，ir。
根据高斯-马尔科夫定理，用加权最小二乘法求得参数的最佳线性无偏估计（BLUE）
i1，i2，…，ir=i1，i2，…，ir
i1，i2，…，ir=
[X′i1，i2，…，ir-1i1，i2，…，irXi1，i2，…，ir]-1
X′i1，i2，…，ir×-1i1，i2，…，irYi1，i2，…，ir
当n，r及i1，i2，…，ir给定时，矩阵Xi1，i2，…，ir与Σi1，i2，…，ir都是已知的，所以为了便于使用，将参数的最佳线性无偏估计简化如下形式
（i1，i2，…，ir）=∑rj=1D0（i1，i2，…，ir，n，r，j）lntij
（i1，i2，…，ir）=∑rj=1D1（i1，i2，…，ir，n，r，j）lntij
其中：D0（i1，i2，…，ir，n，r，j）称为μ（i1，i2，…，ir）的最佳线性无偏估计系数;D1（i1，i2，…，ir，n，r，j）称为σ（i1，i2，…，ir）的最佳线性无偏估计系数。
因此，在R=（R1，R2，…，Rr）下，ZZ分布参数最佳线性无偏估计的一般表达形式为
=∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤nPi1，i2，…，ir·（i1，i2，…，ir）（11）
=∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤nPi1，i2，…，ir·（i1，i2，…，ir）（12）
上面的表达式涉及到系数D0（i1，i2，…，ir，n，r，j），D1（i1，i2，…，ir，n，r，j）和概率P{I1=i1，I2=i2，…，Ir=ir}≡Pi1，i2，…，ir的计算，本文仅对下面两种情况，给出参数估计的简化形式，并构造相应的数表计算公式，便于使用。
3.1 R=[n-r，0，0，…，0]1×r时的参数估计
当移出数据为R=[n-r，0，0，…，0]1×r时，n个元件参加试验，在第一个元件失效时刻一次随机移出n-r个未失效元件，剩下的r-1个元件继续参加试验直至全部失效。ti1<ti2<…<tir可以看成来自容量为n的次序统计量t1<t2<…<tn中固定第一个次序统计量后，从余下n-1元件中随机抽取r-1个的子样本，且该子样本有Cr-1n-1种等可能。所以此时有
P{I1=i1，I2=i2，…，Ir=ir}≡Pi1，i2，…，ir=1/Cr-1n-1
由式（11）和式（12）得ZZ分布参数估计可表示为如下形式
R=（n-r，0，0，…，0）=1Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤n（i1，i2，…，ir）
R=（n-r，0，0，…，0）=1Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤n（i1，i2，…，ir）
令
1Cr-1n-1∑（i1，i2，…，ir）∈G1=ii<i2<…<ir≤nD0（ii，i2，…，ir，n，r，j）=D′0，R=（n-r，0，0，…，0）（n，r，j）
1Cr-1n-1∑（i1，i2，…，ir）∈G1=ii<i2<…<ir≤nD1（ii，i2，…，ir，n，r，j）=D′1，R=（n-r，0，0，…，0）（n，r，j）
则R=[n-r，0，0，…，0]1×r时ZZ分布的最佳线性无偏估计如下
R=（n-r，0，0，…，0）=1Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤n∑rj=1D0（i1，i2，…，ir，n，r，j）lntij=
∑rj=11Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤nD0（i1，i2，…，ir，n，r，j）lntij=
∑rj=1D′0，R=（n-r，0，0，…，0）（n，r，j）lntij
R=（n-r，0，0，…，0）=1Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤n∑rj=1D1（i1，i2，…，ir，n，r，j）lntij=
∑rj=11Cr-1n-1∑（i1，i2，…，ir）∈G1=i1≤i2≤…≤ir≤nD1（i1，i2，…，ir，n，r，j）lntij=
∑rj=1D′1，R=（n-r，0，0，…，0）（n，r，j）lntij
為了便于计算参数的估计，其系数D′0，R=（n-r，0，0，…，0）（n，r，j）和D′1，R=（n-r，0，0，…，0）（n，r，j）在n=5到n=20时已由计算机给出。
3.2 R=[1，1，…，1，0，0，…，0]1×r时的参数估计
在每次失效时刻XRj∶r∶n，j=1，2，…，r分别移出1个未失效产品的情况：n个元件参加试验，设有r个失效，且在每次失效时刻移出1个未失效元件。在R=[1，1，…1，0，0，…，0]1×r时，每种情况可能出现的概率P{I1=i1，I2=i2，…，Ir=ir}≡Pi1，i2，…，ir不再相等，但其计算可借助计算机完成，最终体现在后面的数表中。 由式（11）和式（12），移出數据为R=[1，1，…，1，0，…，0]1×r的ZZ分布的最佳线性无偏估计表示如下
R=（1，…，1，0，…，0）=∑rj=1D′0，R=（1，…，1，0，…，0）（n，r，j）lntij
R=（1，…，1，0，…，0）=∑rj=1D′1，R=（1，…，1，0，…，0）（n，r，j）lntij
其中，
D′0，R=（1，…，1，0，…，0）（n，r，j）=∑（i1，i2，…，ir）∈G1=ii<i2<…<ir≤nPi1，i2，…，ir×
D0（ii，i2，…，ir，n，r，j）
D′1，R=（1，…，1，0，…，0）（n，r，j）=∑（i1，i2，…，ir）∈G1=ii<i2<…<ir≤nPi1，i2，…，ir×
D1（ii，i2，…，ir，n，r，j）
为了便于计算参数的估计，其系数D′0，R=（1，1，…，1，0，…，0）（n，r，j）和D′1，R=（1，1，…，1，0，…，0）（n，r，j）在n=5到n=20时已由计算机给出。
4 数值模拟
基于Monte Carlo模拟，给出以上两种方法下ZZ分布参数的估计结果。根据文[4]给出的算法，产生ZZ分布下逐次定数截尾数据的步骤如下：
1）产生r个来自均匀分布U（0，1）的独立随机样本W1，W2，…，Wr;
2）在预先设定的逐次定数截尾移出数据
R=[R1，R2，…，Rr]下，令
Vj=W1/（j+Rr+Rr-1+…Rr-j+1）j，j=1，2，…，r;
3）再令Uj=1-VrVr-1…Vm-i+1，j=1，2，…，r，则U1，U2，…，Ur是来自均匀分布U（0，1）的逐次定数截尾数据;
4）令Xj=F-1（Uj），j=1，2，…，r，得到来自ZZ分布的逐次定数截尾试验数据，其中F-1（·）是ZZ分布函数的反函数，
Xj=（β·ln（1-ln（1-Uj）））1/m，j=1，2，…，r。
设定两组不同样本量和两组不同的移出数据，根据3.1和3.2节方法，构造最佳线性无偏估计的系数值分别如表1和表2所示。
应用上述数据，模拟得参数估计的结果如表3所示。
从表3中可以看出，最佳线性无偏估计结果和近似极大似然估计结果都较接近参数真值，且最佳线性无偏估计结果较优。因此，二者均可应用于寿命服从ZZ分布、逐次定数截尾寿命试验下的参数估计。
5 结 论
本文基于逐次定数截尾寿命试验模型，运用近似极大似然法和最佳线性无偏估计法，对产品寿命服从ZZ分布的参数进行了统计分析。主要构建了r维不可观测但可求分布的随机向量，据此将收集到的逐次定数截尾样本数据转化为一定概率分布条件下的顺序统计量，进而给出ZZ分布参数的最佳线性无偏估计，并可依据构造出的数表再进行简单计算即得到参数的估计值。且此法具有的解析表达形式、运用数表计算方便等优点，易可应用到Weibull分布或其他指数型寿命分布中。
参 考 文 献：
[1] COHEN， CLIFFORD A. Progressively Censored Samples in Life Testing[J].Technometrics， 1963， 5（3）：327.
[2] COHEN A C. Progressively Censored Sampling in the Three Parameter Log-normal Distribution[J]. Technometrics， 1976， 18（1）：99.
[3] NORGAARD C N J. Progressively Censored Sampling in the Three-parameter Gamma Distribution[J]. Technometrics， 1977， 19（3）：333.
[4] SANDHU， BALAKRISHNANR N. A Simple Simulational Algorithm for Generating Progressive Type-II Censored Samples[J]. The American Statistician， 1995， 49（2）：229.
[5] BALASOORIYA U， SAW S L C. Reliability Sampling Plans for the Two-parameter Exponential Distribution under Progressive Censoring[J]. Journal of Applied Statistics， 1998， 25（5）：707.
[6] TSE S K， YANG C. Reliability Sampling Plans for the Weibull Distribution under Type-II Progressive Censoring with Binomial Removals[J]. Journal of Applied Statistics， 2003， 30（6）：709.
[7] BALAKRISHNAN N， AGGARWALA R. Progressive Censoring[M]. Birkhauser， 2000.
[8] BALAKRISHNAN N， JIE Mi. Existence and Uniqueness of the MLEs for Normal Distribution Based on General Progressively Type-II Censored Samples[J]. Statistics and Probability Letters， 2003， 64（4）：407. [9] BALAKRISHNAN N， KANNAN N， LIN C T. Point and Interval Estimation for Gaussian Distribution， Based on Progressively Type-II Censored Samples[J]. IEEE Transactions on Reliability， 2003， 52（1）：90.
[10]BALAKRISHNAN N， LIN C T. On the Distribution of a Test for Exponentiality Based on Progressively Type-II Right Censored Spacings[J]. Journal of Statistical Computation & Simulation， 2003， 73（4）：277.
[11]BALAKRISHNAN N， NG H K T， KANNAN N. A Test of Exponentiality Based on Spacings for Progressively Type-II Censored Data[J]. Goodness-of-it Tests and Model Validity， 2002：89.
[12]BALAKRISHNAN N， KANNAN N， LIN C T， et al. Inference for the Extreme Value Distribution under Progressive Type-II Censoring[J]. Journal of Statistical Computation and Simulation， 2004， 74（1）：25.
[13]SOLIMAN， A. A. Estimation of Parameters of Life from Progressively Censored Data Using Burr-XII Model[J]. IEEE Transactions on Reliability， 2005， 54（1）：34.
[14]MAHMOUD M A W， MOSHREF M， YHIEA N M， et al. Progressively Censored Data from the Weibull Gamma Distribution Moments and Estimations[J]. J Appl Probab Stat， 2014， 3 （1）：45.
[15]SINGH S， TRIPATHI Y M， WU S J. On Estimating Parameters of a Progressively Censored Log-normal Distribution[J]. Journal of Statistical Computation and Simulation， 2015， 85（6）：1071.
[16]ALIM A M， JAHEEN Z F. Estimation of Parameters for the Exponentiated Pareto Distribution Based on Progressively Type-II Right Censored Data[J]. Journal of the Egyptian Mathematical Society， 2016， 24（3）.
[17]賀国芳， 瞿荣贞编写. 可靠性数据的收集与分析[M]. 北京：国防工业出版社， 1995.41.
[18]徐晓玲， 王蓉华. Weibull 分布逐步增加的Ⅱ型截尾试验的统计分析[J].强度与环境， 2003， 30（2）：31.
XU Xiaoling， WANG Ronghua. Statistical Analysis of Type Ⅱ Censoring Test with Increasing Weibull Distribution[J]. Strength and Environment， 2003， 30 （2）： 31.
[19]王炳兴. Weibull分布基于定数逐次截尾寿命数据的统计分析[J].科技通报， 2004（6）：488.
WANG Bingxing. Statistical Analysis of Weibull Distribution Based on Constant Successive Censored Life Data[J]. Science and Technology Bulletin， 2004 （6）： 488.
[20]杨君慧， 师义民， 曹弘毅. 逐步增加Ⅱ型截尾试验下广义指数分布的统计分析[J]. 统计与决策， 2014（16）：28.
YANG Junhui， SHI Yimin， CAO Hongyi. Statistical Analysis of Generalized Exponential Distribution under Gradually Increasing Type II Censoring Test[J]. Statistics and Decision Making， 2014 （16）： 28.
[21]杨敏. 删失数据场合双参数Rayleigh分布的参数估计[D]. 广西师范学院， 2015：12.
[22]罗嘉成， 陈勇明. 逐次定数截尾下Lomax分布的贝叶斯分析[J]. 成都信息工程大学学报， 2016， 31（3）：323.
LUO Jiacheng， CHEN Yongming. Bayesian Analysis for Lomax Distribution under Progressively Type-Ⅱ Censored Data[J]. Journal of Chengdu University of Information Technology， 2016，31（3）：323.
[23]张宁. 截尾数据下ZZ分布参数估计及性质[D]. 哈尔滨：哈尔滨理工大学， 2019.
（编辑：温泽宇）