双边定时截尾下BS分布的统计分析

马怡舟,徐晓岭,顾蓓青

(上海对外经贸大学 a.国际经贸学院; b.统计与信息学院，上海 201620)

在可靠性寿命试验中，产品失效有各种原因，其中有一种是由于在某种周期应力作用下，产品产生裂缝，如果裂缝长度达到或者超过某一值时，则产品失效。两参数Birnbaum-Saunders模型是概率物理方法中的一个重要失效分布模型，Birnbaum和Sauders 1969年在文献[1]中研究主因裂纹扩展导致的材料失效过程中推导出来一个两参数疲劳寿命分布，这一失效分布模式在机械产品可靠性研究中应用广泛，主要应用于疲劳失效研究，对于电子产品性能退化失效分析也有重要应用。该两参数疲劳寿命分布的分布函数、密度函数分别为

t≥0,α>0,β>0

(1)

t≥0,α>0,β>0

(2)

其中，Φ(x),φ(x)分别为标准正态分布的分布函数和密度函数。参数α和β分别称为形状参数和刻度参数，记此分布为BS(α,β)。

值得指出的是在国内，文献[2]根据混凝土在疲劳荷载下的损伤机理，提出以应变做为表征其损伤程度的量度，认为混凝土的疲劳损伤过程具有马尔可夫性，通过求解FOKKER-PLANCK方程，导出了在指定时间下疲劳损伤分布服从两参数Birnbaum Saunders疲劳寿命分布。利用文献[3]所给出的在直接拉伸作用下混凝土抗拉疲劳寿命的试验数据资料，认为用两参数BS疲劳寿命分布BS(α,β)拟合比用对数正态分布更合适。

由于Birnbaum-Saunders疲劳寿命分布是从疲劳过程的基本特征出发导出的，因此它比常用寿命分布如威布尔分布,对数正态分布更适合描述某些由于疲劳而引起失效的产品寿命规律。此分布已经成为可靠性统计分析中的常用分布之一。关于两参数BS疲劳寿命分布BS(α,β)的统计分析已有众多的文献研究。文献[4-5]说明了在较文献[1]更弱的条件下产品的寿命分布仍服从BS(α,β)，因此用BS(α,β)分布拟合这类产品的寿命分布比常用的Weibull、对数正态分布更合适。在完全样本情形，文献[6]讨论了BS(α,β)分布参数的点估计，导出了参数的极大似然估计。文献[7]利用模拟方法和极大似然估计的渐近正态性讨论了参数的区间估计和假设检验。文献[8]研究了BS(α,β)分布的对数线性模型，导出了参数的极大似然估计和最小二乘估计，利用极大似然估计的渐近正态性给出了参数的近似置信区间。文献[9]讨论了可靠度函数的区间估计方法。文献[10]讨论了BS(α,β)分布在截尾试验情形下的统计分析，给出了参数的拟最小二乘估计和极大似然估计，另外还给出了刻度参数β的区间估计。文献[11]研究了缺失数据场合下两参数疲劳寿命分布BS(α,β)刻度参数β的区间估计。文献[12]提出了一个新的枢轴量来构造刻度参数β的经典置信区间，并且具有较为简单的显式表达式。文献[13]在失效机理保持不变的条件下，讨论两参数BS疲劳寿命分布环境因子的区间估计问题。文献[14]研究了形状参数以及均值、分位数、可靠度等可靠性指标的广义区间估计。文献[15]通过Monte Carlo模拟说明文献[11-12]所给出的两种方法可能无法得到两参数疲劳寿命分布BS(α,β)刻度参数β区间估计，同时指出文献[14]在利用广义枢轴量法给出刻度参数β以及参数函数的置信区间过程中存在错误，并用反例进行了说明，同时给出了正确的证明。文献[16]应用回归分析方法给出了两参数疲劳寿命分布BS(α,β)参数的最小二乘估计和形状参数α的简单易算的区间估计方法，并用计算机随机模拟方法研究了区间估计的效果，模拟结果显示效果非常好。文献[17]针对两参数疲劳寿命分布BS(α,β)利用形状参数α的区间估计给出了分布变异系数的区间估计和假设检验方法。文献[18]针对逐步增加的II型截尾数据给出了参数的近似极大似然估计。文献[19]通过对数变换给出了求两参数疲劳寿命分布BS(α,β)在全样本场合下参数的对数矩估计，并通过大量Monte Carlo模拟比较了各种点估计方法的精度，基于对数变换通过一阶泰勒展开，将两参数疲劳寿命分布BS(α,β)近似看作两参数对数正态分布，由此得到了刻度参数β和形状参数α的近似区间估计，且区间估计都存在，通过Monte Carlo模拟比较发现所给出的近似方法比原有方法更精确，最后通过若干实例说明方法的可行性。

本文针对两参数Birnbaum-Saunders疲劳寿命分布产品在双边定时截尾寿命试验数据场合下，利用“半个产品失效”的想法，并利用Balakrishan提出的一阶泰勒展开，将似然方程作适当的近似，得到了两个参数的近似极大似然估计，最后通过一实际案例说明方法的可行性。

1 双边定时截尾寿命试验场合下参数的点估计

设产品的寿命T服从两参数BS(α,β)分布，如果做双边定时截尾寿命试验，定时截尾时间设为τ1,τ2其中τ1<τ2，在τ1前共有r1-1个产品失效，在τ2前共有r2个产品失效，次序失效时间为：τ1≤t(r1)≤t(r1+1)≤…≤t(r2)≤τ2，也就是说在时间区间[τ1,τ2]内共有k=r2-r1+1个产品失效。此时，似然函数为(其中C+为正常数)：

L(α,β)=C+[F(τ1;α,β)]r1-1[1-

lnL(α,β)=lnC++(r1-1)ln[F(τ1;α,β)]+

(3)

(4)

式(3)与式(4)为含参数α,β的二元超越方程组，要从中求解得参数α,β的极大似然估计是一个非常复杂的问题。为解决这一问题，下面通过近似处理，通过求解仅含参数β的一元超越方程得到β的近似极大似然估计，进而得到参数α的近似极大似然估计。

似然函数为(其中C+为正常数)：

(5)

lnL(α,β)=lnC++(r1-1)ln[Φ(Z1)]+

(n-r2)ln[1-Φ(Z2)]+

lnL(α,β)=lnC+-klnα-klnβ+(r1-1)ln[Φ(Z1)]+

(n-r2)ln[1-Φ(Z2)]+

(6)

(7)

(8)

(9)

(10)

观察双边定时截尾寿命试验的数据形式，在时刻τ1前有r1-1个产品失效，而在时间t(r1-1)与τ1之间没有产品失效，在此将其视作有“半个产品失效”。在时刻τ2前有r2个产品失效，而在时间t(r2)与τ2之间没有产品失效，在此将其视作有“半个产品失效”。

利用上述泰勒展开对式(9)化简得：

(r1-1)Z1(a(1)-b(1)Z1)-k=0

即

也即

(11)

其中，

利用上述泰勒展开对式(10)化简得：

即

也即

(12)

(13)

(14)

特别地：当r1=1时，为一般定时截尾寿命试验场合，即产品试验做到时间τ0为止，在τ0前有r个产品失效，此时τ1=0,τ2=τ0,r2=r。类似在可以给出一般定时截尾寿命试验参数的近似极大似然估计。

2 实例

本实例[21]中的数据集为6061-T6铝合金的疲劳寿命。铝合金的切取位置应平行于轧制方向，振荡频率为18赫兹。该数据集包括101个观测值，最大应力为31 000帕。数据如下：70， 90，96， 97， 99，100，103，104，104，105，107，108，108，108，109，109，112，112，113，114，114，114，116，119，120，120，120，121，121，123，124，124，124，124，124，128，128，129，129，130，130，130，131，131，131，131，131，132，132，132，133，134，134，134，134，134，136，136，137，138，138，138，139，139，141，141，142，142，142，142，142，142，144，144，145，146，148，148，149，151，151，152，155，156，157，157，157，157，158，159，162，163，163，164，166，166，168，170，174，196，212

文献[19]在全样本场合下给出了14种点估计，如表1所示。

表1 全样本场合下参数的点估计

从双边定时截尾寿命试验与一般定时截尾寿命试验的计算结果看与表1中的结果近似。