混合分组数据下恒加试验的参数估计

混合分组数据下恒加试验的参数估计
俞自由;蔡万科
【摘要】本文研究了恒定心力加速寿命试验在获得混合分组数据情况下的参数估计问题.利用EM算法,扶得了参数极大似然估计的迭代式.该结果推广了分组数据场合下一般恒加试验的参数估计.
【期刊名称】《数学杂志》
【年(卷),期】2010(030)002
【总页数】6页(P367-372)
【关键词】加速寿命试验;EM算法;混合分组数据;指数分布
【作者】俞自由;蔡万科
【作者单位】上海财经人学金融学院,上海,200433;上海财经人学金融学院,上海,200433
【正文语种】中文
【中图分类】O212.1
恒定应力加速寿命试验是产品可靠性评定中常用的试验方法.在该试验中，一定数量的产品先被分组，然后分别在不同应力水平下进行寿命试验，试验在进行到一定时间或观测到一定数量的失效产品而获得截尾数据后停止.对此情形下的恒加寿命试验的研究，文献[1,2]有较为详细的论述.
尽管这种加速寿命试验可以减少试验时间，但要获得截尾数据需要对整个试验进行跟踪并记录产品的失效时间，这必然要消耗很多人力和财力.如果在试验过程中将
随时跟踪改为定时检测给定时间区间内的失效个数的话，将会给试验带来很大的方便，这种情况下得到的数据被称为分组数据.与前面提到的截尾数据不同，分组数据只知道产品在给定时间区间内的失效个数，而不知道产品具体的失效时间，因此也给这类试验的统计分析增加了难度.
国内有些文献研究了指数分布场合下分组数据基础上的步加试验，文献[3]对简单步加试验(即两步步加试验)，在相等观测时间间隔获得分组数据的场合下，给出了参数极大似然估计的精确解；文献[4]进一步考虑分组数据场合下简单步加试验中观测时间间隔不相等的情形，利用数值解方法给出极大似然估计的迭代式.对于分组数据场合下一般的恒加试验的统计分析还未有文献涉及.
本文讨论分组数据场合下一般的恒加试验，并做一些改进：在低应力时，产品具有较大的平均寿命，此时，我们只定时检测产品在给定时间内失效的个数；在高应力时，产品失效得快，此时记录产品失效的具体时间.我们称此类试验下获得的数据为：混合分组数据.这样的试验不但可以减少试验的投入，而且也有利于统计分析.在混合分组数据的基础上，本文利用EM算法给出参数极大似然估计的迭代式，并将之应用于解决一个实际的例子
本文剩下部分安排如下：第2节介绍EM算法；第3节解释相关符号和陈述相关假设；第4节给出估计的表达式，并将之应用于解决一个实际的例子.
EM算法是一种迭代方法,最初由Dempster等提出,并且主要是用来求后验分布的众数（即极大似然估计）,它的每次迭代由两步完成：E步（求期望）和M步（极大化）.一般地,以P(θ|Y)表示θ的基于观测数据Y的后验分布密度函数,称为观测后验分布；P(θ|Y,Z)表示添加数据Z后得到的关于θ的后验分布密度函数,称为添加后验分布；P(Z|θ,Y)表示给定θ和观测数据Y下添加数据Z的条件分布密度函数.我们的目的是计算观测后验分布P(θ|Y)的众数,于是EM算法如下进行：记θ(u)为第u+1次迭代开始时后验众数的估计值,则第u+1次迭代的两步为:
E步将lnP(θ|Y,Z)关于随机变量Z求条件期望,即
M步将Q(θ|θ(u),Y)极大化,即找一点θ(u+1),使
如此形成了一次迭代θ(u)→θ(u+1).将上述E步和M步进行迭代直
至‖θ(u+1)−θ(u)‖或者‖Q(θ(u+1)|θ(u),Y)−Q(θ(u)|θ(u),Y)‖达到要求时停止.
下面两个定理讨论EM算法的收敛性.
定理2.1 EM算法的每一次迭代使P(θ(u+1)|Y)≥P(θ(u)|Y)当且仅当
时等号成立.
定理2.2假设一系列EM迭代θ(u)满足
(1)θ(u)收敛于某一个值θ∗,并设p(z|θ,Y)充分光滑.则有
换句话说,如果迭代θ(u)收敛,它们将收敛于P(θ|Y)的稳定点.
(1)确定正常应力水平s0和k个加速应力水平s1＜s2＜···＜sk.这些应力水平一般应有关系:s0＜s1＜s2＜···＜sk
(2)从该批产品中随机选出n个,并分为k个样本,其样本容量分别为n1,n2,···,nk,记n=n1+n2+···+nk,其中第i个样本将安排在si下进行试验.
(3)将k个应力水平分成两组:前k1个和后k2个,k1+k2=k.记前k1个应力水平(s1,s2,···,sk1)的mi,(i=1,2,···,k1)个检测时刻为
为后k2个应力水平确定决定停止试验的k2个失效个数
在前k1个应力水平下的试验中,每到检测时刻τij便去记录产品在(τi(j−1),τij]中失效的个数；而在后k2个应力水平下的试验中,去跟踪产品的失效时间.
(4)假设我们记录了前k1个应力水平下的试验在时间段(τi(j−1),τij],j=1,···,mi有rij 个产品失效.后k2步试验中第i个应力水平下的试验ri个失效产品的具体失效时间为
λi:第i个应力水平si下产品的失效率.θi=1/λi为平均寿命.
ri:第i个应力水平si下总的失效个数,当i=1,2,···,k1时,记
τi:第i个应力水平si下试验的截止时间,当i=1,2,···,k1时,τi=τimi;当i= k1+1,···,k 时,τi=tiri.
假设A 在应力水平si下产品寿命T服从指数分布,其寿命分布函数为
假定B产品的平均寿命θi与应力水平si之间有如下加速模型
其中a,b是待定参数,φ(si)为应力水平si的已知函数,并记φ(si)=φ(i).
定理4.1在上面的试验数据和假设下,参数a,b的极大似然估计可以经下面的迭代算得
其中
当试验为简单步加试验,即k=2时,迭代式具有更简单的表达式
证添加Zijh, h=1,...,rij为第k组产品在时间段(τij−1,τij]上失效的rij个产品的具体失效时间,i=1,···,k,j=1,···,mi.如果第u+1次迭代开始时参数a,b的估计值为
a(u),b(u),记=e−(a(u)+b(u),则具有如下的分布密度函数
因而有
E步
其中ET由式(4.2)给出.
M步由于r＞0,且Q(a,b|a(u),b(u),T)关于a,b任意阶可导,所以Q(a,b|a(u),b(u),T)最大值存在且为其稳定点,对Q(a,b|a(u),b(u),T)关于a,b求一阶导可得方程组
特别地,当步加试验为简单步加试验的时候,即k=2时,有
由于可选择s1,s2,使得矩阵可逆,所以进而有
一般情况下,我们解方程组得a,b满足
这样我们只要求解b就可以了.利用牛顿迭代法可以求得
其中由式4.1定义.当m→∞时,b(u,m+1)→b(u+1).并且当u→∞时,b(u+1)→b.
可见要得到b的估计值,必须经过两轮迭代,这样计算量相当大,求解变得非常麻烦.但是b的估计值可以由下面的一轮迭代式得出
因为由(4.4)式产生的数列是(4.3)式产生的数列的子列,并将收敛于同一数值b.所以定理成立.
例某产品的平均寿命为5×104h,改进设计后其平均寿命有明显提高.为评定新产品在工作电压10V下的平均寿命,现在组织一次加速寿命试验,该试验的前两步获得分组数据；后两步试验获得定数截尾数据.所选的4个加速电压水平、样本容量、失效数和估计结果都列于表4.1.从表4.1可见,利用本文获得的估计式能获得较好的估计结果
【相关文献】
[1]茆诗松,王玲玲.加速寿命试验[M].北京：科学出版社,1997.
[2]张志华.加速寿命试验及其统计分析[M].北京：北京工业大学出版社,2003.
[3]李艳冰,费鹤良.分组数据下简单步进应力加速寿命试验的极大似然估计[J].上海师范大学学
报,2003, 32(1):15–20.
[4]李艳冰,费鹤良.指数分步场合分组数据下简单步加试验的极大似然估计[J].电子产品可靠性与环境试验,2004,31(4):31–34.
[5]张春华,温熙森,陈循.加速寿命试验技术综述[J].兵工学报,2004,25(4):485–489.
[6]张志华,茆诗松.恒加试验简单线性估计的改进[J].高校应用数学学报A辑(中文版),1997,12(4): 417–424.。