基于逐步Ⅱ型删失恒定应力部分加速寿命试验下BurrⅢ分布的统计分析

杜雪萌，张峰源

（贵州民族大学 数据科学与信息工程学院，贵阳 550025）

由于寿命试验产品的寿命时间越来越长，在进行生存分析与可靠性试验时很难在正常情况下获得这些产品的故障数据。因此，为了缩短测试时间，加速性能退化，对所有或部分测试产品施加比平时更严重的压力条件，其被称为加速寿命测试。如果只将部分产品置于加速条件下，则该试验又称为部分加速寿命试验。近年来，有不少文献对加速寿命试验的统计分析进行了研究。Almarashi[1]基于恒定应力部分加速寿命试验，得到了逐步Ⅱ型删失数据下广义倒指数分布的极大似然估计和Bayes估计；Kumar等[2]在恒定应力加速寿命试验下，利用9种常用估计方法对Lindley分布进行了参数估计；武东等[3]基于定数删失数据，对瑞利分布恒定应力加速寿命试验进行了贝叶斯统计分析。考虑到逐步Ⅱ型删失数据的优良性，本文将基于逐步Ⅱ型删失数据，在恒定应力部分加速寿命试验下，对BurrⅢ型分布的形状参数和加速因子进行参数估计。

1 基本假设和模型描述

逐步Ⅱ型删失下BurrⅢ分布恒定应力部分加速寿命试验的具体试验设计方案如下。

假设有n个寿命服从BurrⅢ分布且相互独立的试验样品，将其随机分成2组，其容量分别为n0，n1。其中，n0个样品被安排到正常应力水平S0下进行寿命试验，n1个样品被安排到加速应力水平S1（S0＜S1）下进行寿命试验。在应力水平Sj（j=0，1）下进行逐步Ⅱ型删失试验。当观测到第一个试验样品失效时，从剩余的nj-1个未失效的样品中随机挑选出Rj1个样品退出试验，当观测到第二个试验样品失效时，再从剩余的nj-Rj1-2个未失效的样品中随机挑选出Rj2个样品退出试验，按照这种方法一直试验，当观测到第mj（mj＜nj）个失效样品时，将剩下的个样品全部撤出试验。共m（m=m0+m1）个失效样本。记在应力水平Sj（j=0，1）观测到的失效时刻为xj=（xj1，xj2，…，xjmj）。本文假设从试验中随机剔除的Rji完全随机，保证Rj1+Rj2+…+Rjm=nj-mj即可。

基本假设如下：

（1）在应力水平Sj（j=0，1）下，试验产品失效机理相同。

（2）在应力水平Sj（j=0，1）下，产品的失效时间xji（j=1，2，i=1，2，…，mj）相互独立。

（3）在应力水平S0下，产品的失效时间x0服从BurrⅢ分布，其分布函数、密度函数、可靠性函数和失效率函数分别为

式中：α＞0，θ＞0分别称为尺度参数和形状参数，本篇文章假设尺度参数α已知。

（4）在应力水平S1下，产品的失效率函数为h1（x）=βh0（x），β＞1，式中：β为加速因子。则应力水平S1下产品失效时间x1的分布函数、密度函数、可靠性函数和失效率函数分别为

2 极大似然估计

基于上述模型描述，通过极大似然方法[4]对BurrⅢ型分布的形状参数θ和加速因子β进行估计，可以得到如下定理。

定理1：按照上述逐步Ⅱ型删失恒定应力部分加速寿命试验，BurrⅢ型分布形状参数θ的极大似然估计由方程（9）确定，并且该方程有唯一解，同时，加速因子β极大似然估计式为式（10）。

证明：假设xj=（xj1，xj2，…，xjmj）是应力水平Sj（j=0，1）下服从BurrⅢ型分布的逐步Ⅱ型删失数据，记x=（x0，x1）为恒定应力寿命试验下服从BurrⅢ型分布的逐步Ⅱ型删失数据。Rj=（Rj1，Rj2，…，Rjmj）为应力水平Sj（j=0，1）下对应的删失模式，参数α已知，则关于x的似然函数为

式中：M

将式（1）、式（2）、式（7）和式（8）带入式（11）可得

对L（x）取对数后分别对θ，β求偏导，并分别令其偏导数为0可得

参数θ显然无法由式（13）得到显式解，为此，探求此方程在（0，+∞）上是否具有唯一解。

显然g（θ）′＜0，同时因此方程（13）在（0，+∞）上有唯一解。

式（13）没有显式表达式，但有唯一解，故考虑应用Newton-Raphson近似法求出近似解。记此解为参数θ的极大似然估计近似值θ^M，将θ^M带入式（14）中可得到参数β的极大似然估计近似值β^M。

3 贝叶斯估计

在一些实际情况下，参数值的信息可以以独立的方式获得。因此，假设参数先验独立，设加速度因子β的先验分布为无信息先验，即

取θ的先验分布为伽马分布，其形状参数与尺度参数分别为a和b，概率密度函数为

因此，参数θ和β的联合先验可表示为

结合式（12）和式（18）可得θ和β的联合后验密度分布函数为

参数θ和β的全条件后验分布函数分别表示为

基于上述全条件后验密度分布函数式（20）和式（21），本文考虑使用Metropolis-Hastings（MH）方法产生Markov Chain Monte Carlo（MCMC）样本并估计参数。选取提议分布为正态分布和具体算法如下：

（1）给出初始值θ(0)，β(0)；

（2）令i=1；

（5）计算θ（i）和β（i）；

（6）令i=i+1；

（7）重复步骤（2）—（6）M次，生成样本θ(1)，θ(2)，…，θ(M)与β(1)，β(2)，…，β(M)。

下面利用基于上述方法所产生的样本θ(1)，θ(2)，…，θ(M)与β(1)，β(2)，…，β(M)对参数θ和加速因子β进行估计。为方便书写，记γ＝（θ，β）。

3.1 平方损失函数下的估计

设δ是参数γ的任一决策函数，则平方损失下L1=（θ-δ）2，γ的贝叶斯估计为

3.2 广义熵损失函数下的估计

3.3 对称熵损失函数下的估计

不同损失函数下，γ估计的后验均方误差为

4 随机模拟

利用Monte-Carlo方法通过R软件产生一个服从BurrⅢ分布的恒定应力部分加速寿命试验逐步Ⅱ型删失样本，具体步骤如下：①产生2个容量分别为n0，n1且服从均匀分布U0（0，1），U0（0，1）的独立同分布样本，并升序排列为U01，U02，…，U0n0和U11，U12，…，U1n1；②当α=2时，设θ=1.5，β=1.2

则X01，…，X0n0，X11，…，X1n1是一个恒定应力部分加速寿命试验下服从参数θ=1.5，加速因子β=1.2的BurrⅢ分布的独立同分布样本；③根据不同样本量n0，n1和观测次数m0，m1，随机生成相应的删失模式R0=（R01，R02，…，R0m0），R1=（R11，R12，…，R1m1）。按所生成的R1，R2进行随机抽样，获得的失效数据X0（1），…，X0（m0），X1（1），…，X1（m1），为服从BurrⅢ分布的恒定应力部分加速寿命试验逐步Ⅱ型删失样本。

基于上述所生成的删失样本，取参数初值θ=1.8，β=1.4，设定精度e=0.000 1，则由式（9）和式（10）可求得参数θ和加速因子β极大似然估计的近似值重复以上过程1 000次，可求得的均值与均方误差。根据上述MH算法步骤，令M=20 000，生成20 000个θ和β，利用所生成的样本θ(1)，…，θ(20000)，β(1)，…，β(20000)代入式（22）—（24），可得到在3种损失函数下的θ和β的估计值，其中，取广义熵损失函数中参数c为0.8和1。利用式（25）可求得上述3种损失函数下估计的后验均方误差（MSE）。

表1和表2的模拟结果表明，在样本量相同的情况下，观测次数不同时，观测次数越多时，参数估计的均方误差越小，估计效果越好；在观测次数相同，样本量n不同时，n越大，参数估计的MSE越小。2种估计方法得到的结果都比较稳健。

表1 参数θ=1.5，β=1.2，不同样本量n，不同观测次数m时的θ的模拟结果

5 结论

本文基于恒定应力部分寿命试验，在逐步Ⅱ型删失数据下，首先运用极大似然法对BurrⅢ分布的形状参数和加速因子进行估计得到了形状参数的估计方程式，并证明该方程式的解是唯一存在的，由于形状参数的估计方程式无显性表达式，采用Newton-Raphson法求出近似解。其次基于MH算法得到了不同损失函数下未知参数的贝叶斯估计式。最后运用Monte-Carlo模拟方法，对在不同情况下的形状参数与加速因子的估计值进行了模拟比较。结果表明，定义的极大似然估计和贝叶斯估计的MSE均随样本量的增加而减小，2种估计效果近似。