近年来已有一些文献研究了竞争失效产品加速寿命试验的统计分析与优化设计。文献[1]研究了广义指数分布竞争失效产品步进应力加速(步加)寿命试验的统计分析,给出了分布参数的点估计和置信区间。文献[2]在失效原因服从对数正态分布的情形下,利用EM算法研究了竞争失效产品加速寿命试验的极大似然分析。文献[3]基于Copula函数建立了恒定应力加速(恒加)寿命试验中相依竞争失效模型,给出了二维和多维失效模式下指数分布参数的极大似然估计和最小方差估计。文献[4]在竞争失效原因相依的情形下应用Copula函数连接了边缘分布和联合分布函数,给出了未知参数的极大似然估计,并简单介绍了正确建立Copula连接函数的方法。文献[5]在逐步Ⅱ型混合截尾恒加寿命试验下,给出了指数分布竞争失效产品未知参数及可靠性指标的极大似然估计和贝叶斯估计。文献[6]在逐步Ⅰ型混合截尾下,建立了威布尔产品恒加试验竞争失效模型,给出了模型参数的极大似然估计、渐近置信区间、贝叶斯估计和近似可信区间。

加速寿命试验的特点是把产品放在高于正常应力水平下进行试验,在获得失效数据后,利用产品寿命和应力水平之间的关系(加速模型),对产品在正常应力水平下的寿命特征进行统计推断。但在工程实际中,产品寿命和应力水平之间的关系不一定总是已知的,例如测试对象是新产品。在这种情形下,可考虑部分加速寿命试验(PALT)。PALT的统计分析不需要加速模型,它通过引入加速因子,将高应力水平下的寿命特征外推到正常应力水平下,进而对产品在正常应力水平下的寿命特征进行统计推断^[7]。近年来PALT已被应用于产品的可靠性分析,见文献[7-10]。但上述PALT的统计推断方法仅适用产品只有一个失效机理的情况,而实际上很多产品包含多个失效机理。当产品具有多个竞争失效机理时,如何利用PALT研究其寿命特征的统计分析具有重要的理论意义和应用价值。鉴于此,本文将逐步混合截尾方式应用于加速寿命试验,在逐步Ⅰ型混合截尾下,讨论Pareto分布竞争失效产品恒定应力部分加速寿命试验的统计推断问题。 文中给出了模型参数及加速因子的极大估计(MLE)和贝叶斯估计(BE)。利用Bootstrap方法及贝叶斯理论分别获得了模型参数的Bootstrap置信区间、最高后验概率密度(HPD)置信区间。最后利用Monte Carlo方法给出数值例子,比较了参数的MLE和BE的优劣性,并讨论了样本量和分配比例对估计精度的影响。

1 基本假设和模型描述

选择2个应力水平S₀＜S₁,其中S₀为正常应力水平,S₁为加速应力水平。基本假设如下:

1) 产品失效由且仅由2个失效机理之一引起。T₁和T₂分别代表机理1和机理2的发生时间,且它们相互独立,产品的寿命T=min(T₁,T₂)。

2) 在应力水平S₀下,机理j(j=1,2)的发生时间服从Pareto分布,其可靠度函数和失效率函数分别为

(1)

(2)

式中，θ_j(θ_j>0)是尺度参数,α_j(α_j>0)是形状参数且已知。

3) 在应力水平S₁下,机理j(j=1,2)的失效率函数为

(3)

式中，λ(λ>1)是加速因子。

根据失效率函数和可靠度函数的关系,可得应力S₁下机理j(j=1,2)的可靠度函数为

(4)

逐步Ⅰ型混合截尾恒定应力部分加速寿命试验模型可描述如下:从n个具有上述失效机理的产品中,选取n₀=n(1-r),0＜r＜1)个产品放在S₀下进行寿命试验,剩余的n₁=nr个产品放在S₁下进行试验,r称为分配比例。在试验过程中,应力水平保持不变。预先设定失效数m＜min(n₀,n₁)和时间τ,在应力水平S_i(i=0,1)下同时进行如下试验:观测到第一个失效时刻t_i1时,从剩下的n_i-1个未失效产品中随机移走R_i1个。类似地,在第二个失效时刻t_i2,从剩余的n_i-2-R_i1个未失效产品中随机移除R_i2个,试验继续。若第m个失效时刻t_im出现在时刻τ之前,则试验在时刻t_im终止,并将剩下的R_im个产品全部移离试验。反之,如果在时间τ之前只观测到了J_i (J_i＜m)个产品失效,则试验在τ时刻停止,并将剩下的R_iJi^*=n_i-J_i-个未失效产品全部移除。因此,在应力水平S_i(i=0,1)下,试验的停止时间可表示为T_i*=min｛τ,t_im｝,得到的截尾样本有2种情形:

情形1 t_i1＜t_i2＜…＜t_im＜τ；

情形2 t_i1＜t_i2＜…＜t_iJi＜τ＜t_i(Ji+1)＜t_i(Ji+2)＜…＜t_im。

令D_i代表应力水平S_i(i=0,1)下时刻τ之前的总失效数。那么,在情形1下T_i*=t_im,D_i=m,在情形2下T_i*=τ,D_i=J_i。在时刻T_i*移走的产品数为R_i*=n_i-D_i-。

若观测样本记为(t_i1,c_i1),(t_i2,c_i2),…,(t_iDi,c_iDi),其中c_il∈｛1,2｝,l=1,2,…,D_i为应力水平S_i下引起第l个失效的机理编号,则在应力水平S_i(i=0,1)下由机理j(j=1,2)引起的失效样本的似然函数为

(5)

式中，δ_j(c_il)=1,当c_il=j,δ_j(c_il)=0,其他,j=1,2。R_ik(i=0,1;k=1,2,…,D_i)为应力水平S_i下第k次失效后移走的产品数。将(1)~(4)式分别代入(5)式,可得到应力水平S₀和S₁下机理j(j=1,2)引起的失效样本似然函数为。其中n_ij=。因此,机理j(j=1,2)引起的总失效样本的似然函数为

由于机理1和机理2的发生时间相互独立,所以全似然函数为

(6)

式中，。

2 极大似然估计和Bootstrap置信区间

首先讨论参数的极大似然估计。对(6)式两边取对数得

(7)

进一步对(7)式分别关于λ和θ_j(j=1,2)求导并令其等于零,可得

(8)

(9)

将(8)、(9)式联立,可得θ_j及λ的极大似然估计分别为

(10)

(11)

式中，t_j1^*是由机理j(j=1,2)引起的第一个失效时间。

下面利用Bootstrap方法来构造参数的Bootstrap置信区间,具体步骤如下:

1) 利用原始的逐步Ⅰ型混合截尾Pareto分布竞争失效产品恒定应力部分加速寿命试验的数据,计算λ,θ₁,θ₂的极大似然估计值,。

2) 将,作为参数产生样本,并根据新样本,利用(9)和(10)式计算3个参数的极大似然估计值,记为。

3) 重复步骤2)N次,将得到的所有按升序排列,即为λ,θ₁,θ₂的Bootstrap样本:

θ_j的置信度为1-γ的Bootstrap P置信区间是所有形如的区间中长度最短的区间。由θ_j(j=1,2)的定义域和(10)式可知。故θ_j(j=1,2)无法落入形如的任何一个区间内。基于此,本文考虑θ_j(j=1,2)的另一种Bootstrap置信区间即Studentized-t(Stud-t)置信区间。

根据文献[11],分别针对λ,θ_j(j=1,2),基于步骤3)得到Bootstrap样本和统计量,其中为参数θ的极大似然估计,是极大似然估计的渐近方差,k=1,2,…N。计算次序统计量如下

则λ,θ_j(j=1,2)的置信度为1-γ的Stu-t置信区间分别为

3 贝叶斯估计和HPD置信区间

本节将利用贝叶斯方法来估计未知参数和加速因子。设λ,θ₁,θ₂的联合先验为无信息先验π(λ,θ₁,θ₂)∝1/λθ₁θ₂。结合(6)式,计算λ,θ₁,θ₂的联合后验密度函数为

(12)

式中，Θ=｛(λ,θ₁,θ₂)|λ>1,0＜θ_j＜t_j1^* (j=1,2)｝。

将(12)式关于θ₁,θ₂积分得到λ的边际后验密度函数具有如下形式

(13)

因此,给定λ时，θ₁,θ₂的条件后验分布是幂律分布,密度函数分别为

显然,无法根据(13)式求得平方损失函数下λ的贝叶斯估计。因此,采用自适应拒绝抽样法来解决该问题。我们首先证明λ的完全条件后验密度函数π(λ|θ₁,θ₂)是对数凹函数。给定θ₁,θ₂,对λ的完全条件后验密度函数取对数得

上式关于λ求二阶导数有∂²lnπ(λ|θ₁,θ₂)/∂λ²=-(D₁-1)/λ²。显然上式右端小于0,故λ的完全条件后验密度函数是对数凹函数。下面通过自适应拒绝抽样法获得平方损失函数下λ,θ_j(j=1,2)的贝叶斯估计,并构造对应的HPD置信区间:

1) 给定θ₁和θ₂,利用自适应拒绝抽样法产生λ。

2) 利用生成的λ,分别产生服从后验分布密度函数π(θ₁|λ,t)和π(θ₂|λ,t)的θ₁和θ₂。

3) 重复步骤1)和步骤2) 1 000次,得到样本(λ^(k),θ₁^(k),θ₂^(k)),k=1,2,…,1 000。则在平方损失函数下,λ,θ_j(j=1,2)的贝叶斯估计分别为。

4) 将得到的所有λ^(k),θ₁^(k),θ₂^(k))(k=1,2,…,1 000)分别按升序排列为(λ^*[1]＜λ^*[2]＜…＜λ^{*[1 000]})和(θ_j^*[1]＜θ_j^*[2]＜…＜θ_j^{*[1 000]})(j=1,2)。计算所有可能的置信区间[λ^*［p],λ^{*[p+1 000(1-γ)]}］和[θ_j^*［p],θ_j^{*[p+1 000(1-γ)]}］(j=1,2;p=1,2,…,1 000-1 000(1-γ))。

5) 计算步骤4)中所得区间的宽度,其中宽度最小的区间即为参数的HPD置信区间。

4 数值模拟

给定参数n、m、τ、α_i、θ_i(i=1,2)、r和移走方案R_k1=R_k2=…R_k(m-1)=1,R_km=n-2m+1(j=0,1)。根据以下步骤生成逐步Ⅰ型混合截尾Pareto分布竞争失效数据:

1) 分别产生4组逐步Ⅰ型截尾样本U₁₁⁽⁰⁾,U₁₂⁽⁰⁾,…,U_1m⁽⁰⁾;U₂₁⁽⁰⁾,U₂₂⁽⁰⁾,…,U_2m⁽⁰⁾; W₁₁⁽¹⁾,W₁₂⁽¹⁾,…,W_1m⁽¹⁾; W₂₁⁽¹⁾,W₂₂⁽¹⁾,…,W_2m⁽¹⁾,它们都来自均匀分布总体U(0,1)。

2) 对i=1,2和1≤j≤m,通过U_ij⁽⁰⁾=1-θ_i^α_i t_0ij^-α_i 计算t_01j和t_02j,并且把min(t_01j,t_02j)记为t_0j^*。若t_0j^*=t_01j,记δ_0j^*=1; 反之记δ_0j^*=2。

3) 对i=1,2和1≤j≤m,由W_ij⁽¹⁾=1-θ_i^λα_i t_1ij^-λα_i 计算t_11j和t_12j,记min(t_11j,t_12j)为t_1j^*。若t_1j^*=t_11j,记δ_1j^*=1;反之记δ_1j^*=2。

4) 若t_im^*＜τ,则D_i=m; 否则查找满足条件t_iDi^*＜τ＜t_i(Di^*+1) 的D_i(i=0,1)。令t_j⁽ⁱ⁾=t_ij^*,δ_j⁽ⁱ⁾=δ_ij^*,R⁽ⁱ⁾=(R_i1,R_i2,…,R_iDi),(1≤k≤D_i)。

(t₁⁽⁰⁾,t₂⁽⁰⁾,…,t_D0⁽⁰⁾)和(t₁⁽¹⁾,t₂⁽¹⁾,…,t_D1⁽¹⁾)分别为S₀和S₁下的逐步Ⅰ型混合截尾失效数据,对应的失效机理编号和移走方案分别为(δ₁⁽⁰⁾,δ₂⁽⁰⁾,…,δ_D0⁽⁰⁾),(R₀⁽⁰⁾,R^*)和(δ₁⁽¹⁾,δ₂⁽¹⁾,…,δ_D1⁽¹⁾),(R₁⁽¹⁾,R^*)。

为了研究样本总数和样本分配比例对估计的影响,给定初始值α₁=0.2,α₂=0.3,θ₁=2.4,θ₂=3,λ=1.8,τ=8,m=0.2n。根据上述算法,分别在r=0.4,0.5和0.6时模拟产生1 000组数据,并基于数据计算λ,θ_j(j=1,2)的(MLE)和(BE)。将2种估计的平均相对误差(ARE)和均方误差(MSE)列于表 1~表 3中,其中ARE和MSE的计算公式分别为。

从表 1~表 3可以看出,参数的MLE和BE的ARE和MSE均随样本量增大而变小;而贝叶斯估计的ARE和MSE均比MLE的ARE和MSE小,可见参数的BE比MLE更精确。另一方面,样本分配比例r与2种估计所对应的ARE和MSE呈负相关关系。为进一步比较区间估计的效果,从n=80时模拟产生的1 000组样本中选取一组列于表 4。根据前文假设,每个应力水平下的预期失效数m=16。

基于表 4的数据,分别计算λ,θ₁,θ₂的MLE、BE以及2种区间估计的上、下限列于表 5。由于θ₁和θ₂的MLE分别是由机理1和机理2引起的第一个失效时间,所以对应的Bootstrap样本值都大于初始值并且相对集中。

由表 5可知,参数的BE更接近初始值。由于参数λ的本身是正数,所以表 5中λ的Stud-t置信区间的下限取为0。比较参数的Stud-t置信区间和HPD置信区间的长度,发现HPD置信区间长度较小,Stud-t置信区间长度较大。这说明在相同的置信度95% 情形下,HPD置信区间优于Stud-t置信区间。

5 结论

基于逐步Ⅰ型混合截尾恒定应力部分加速寿命下Pareto分布竞争失效产品的寿命数据,利用极大似然法和贝叶斯理论获得了分布参数和加速因子的MLE和BE。分别运用Bootstrap法和自适应拒绝抽样法得到了参数和加速因子的Stud-t置信区间HPD置信区间。最后,利用Monte Carlo方法给出数值例子,比较了参数的MLE和BE的优劣性,并讨论了样本量和分配比例对估计精度的影响。模拟计算表明,贝叶斯点估计比MLE更精确;参数及加速因子的MLE和BE的精度随着样本量和分配比例的增大而提高;在相同的置信度下,HPD置信区间优于Stud-t置信区间。