加速寿命数据的贝叶斯建模与分析

汪建均, 杨桂康, 冯泽彪

(南京理工大学经济管理学院, 江苏 南京 210094)

0 引 言

产品质量不仅是企业的生命线,更是企业在全球市场竞争中赢得顾客的关键。近年来,随着各行业在产品质量方面的竞争日益激烈,可靠性作为最重要的质量维度之一,在提高产品寿命、降低保修成本、实现预期产品功能等方面发挥着至关重要的作用。因此,为了减少保修期内的故障,制造商将注意力集中在高质量、高可靠性的产品设计上[1]。可靠性改进是产品质量开发的重要组成部分。许多具有高可靠性的产品,例如汽车、电视机等,通常寿命超过10年以上,这使得在正常操作条件下难以评估产品的可靠性[2]。此时,想要获取产品信息的常规方法是加速寿命试验(accelerated life test,ALT)。ALT是使产品在加速操作条件下更快地失效,从而估计产品在正常运行条件下的寿命。因此,ALT的关键是在试验寿命估计与加速条件之间建立稳健、合适的关系模型,从而在加速条件下拟合数据,并将模型推广到正常操作条件下。然而,在加速条件下,一些测试产品在试验终止时可能不会失效,此时就产生了删失数据[3]。文献[4-6]全面介绍了加速寿命试验建模和数据分析的方法。

以往大多数关于ALT的研究,由于忽略了真实的试验设计条件,寿命试验数据是通过完全随机设计的方法得到的。完全随机设计意味着试验因子组合被随机地应用于每个测试单元。文献[5-7]详细介绍了关于寿命数据的分析。文献[8]举例说明了一个提高金刚石钻头可靠性的Plackett-Burman设计和分析。在这些文献中,假设试验数据是通过完全随机设计得到的。在实际应用生产中,为了节约成本,通常在具有多个样品的试验台上施加应力水平,此时就产生了子抽样结构。在这种情况下,试验台是试验单元,试验样品是观测单元,所获得的寿命试验数据通常不是完全随机设计的。另外,聚类数据、批处理以及裂区试验设计等情况都会导致试验样品未完全随机化[9]。在涉及到非完全随机化的可靠性数据建模时,忽视随机效应可能会产生错误的分析结果,尤其是当试验单元之间的方差大于观测单元间的方差时。许多研究人员已经意识到将随机效应整合到可靠性试验建模与数据分析的重要意义。文献[10]采用了一种两阶段方法,用最大似然估计方法来估计每个所需参数。文献[11]分别比较了I型删失,II型删失和数据未删失3种情况下两阶段方法的性能。由于分成两个阶段分析未知参数,因此针对其未知参数构建共同的似然函数是极其困难的,因此不能对两阶段方法的某些特征(如失效百分位数)进行推断。为此,文献[10, 12-13]提出了一种考虑随机效应的非线性混合模型(nonlinear mixed model, NLMM)方法。该方法不但能够减少未知参数的偏差,并且能够精确计算低分位数的置信区间。在处理包含子抽样结构的样本数据时,文献[14]整合了多重插补方法和异质性模型,该方法进一步提高了置信区间的预测精度。另外,文献[15]提出了一种贝叶斯马尔可夫链蒙特卡罗(Markov chain Monte Carlo,MCMC)建模方法来处理杜邦公司纤维链的随机线轴效应问题。虽然这些模型考虑了随机效应,但是其往往在构建模型前假设尺度参数是固定的,从而忽略了对数位置-尺度分布的尺度参数的检验。

针对加速寿命试验的大多数研究只考虑了尺度参数不相等的情况,而忽略了对于形状参数的分析。文献[2]中假设加速应力通过尺度参数而不是形状参数纳入到模型中。文献[10]所提出的传统两阶段方法都是建立在形状参数相等的前提下。文献[16]在形状参数相等的情况下,给出了分位数置信区间的方法。另外,有一部分学者意识到了形状参数是变量的检验问题。文献[17]将用于设计ALT的最大似然方法扩展到非恒定形状参数的模型中。另外,在指定应力因子的水平下,该设计方案可以最小化最大似然估计的渐近方差。当应力水平不影响形状参数时,文献[18]推导出了II型删失情况下恒加速寿命试验的参数估计量。文献[18]中的方法计算简单,不会出现估计量不存在的问题,同时该方法的偏差也比线性无偏估计的偏差更小。文献[19]提出了一种有效的算法,无论对数寿命数据的标准方差是否与应力因子相关,都可以计算出未知参数的最大似然估计值。文献[20]在构建模型之前首先检验了分组数据的形状参数是否相等,但是并没有考虑形状参数随应力因子的变化而变化的情况。文献[21]假设Weibull寿命分布的尺度参数和形状参数都随应力因子的变化而变化,并且通过自助法得到了低分位数的置信区间。文献[22]认为寿命分布的尺度参数和形状参数都和应力因子有关,并且通过正态近似法得到低分位数的置信区间。文献[2]同时考虑了随机效应和形状参数是变量的问题,并使用Weibull回归模型从ALT中对分位数进行推断。这些文献研究表明,加速应力不仅对位置-尺度分布的位置参数有影响,而且对尺度参数也有一定的影响,但是由于模型的限制,只假设随机效应影响尺度参数。

在ALT分析中同时将随机效应纳入到尺度参数和形状参数中,并且考虑形状参数不相等的研究较少。事实上,不同类型的寿命试验数据(完整或删失数据集、大数据集或小数据集等)可能采用不同的估计方法。例如,对于不同的寿命试验数据,考虑偏差和均方根误差等统计量时,没有一种方法在参数估计方面总是优于其他方法。此外,其他因素也可能会影响参数估计方法的选择,例如计算的简单性[23]。因此探索小样本的分布特性,提高参数对小样本分布特性的可靠性仍然是必要的[24]。NLMM在整合随机效应方面有很大优势,但是NLMM只能处理随机效应服从正态的情况,并且不能计算两个随机效应。另外NLMM模型的积分异常复杂,它的计算精度和积分点选取的数量有关。因此,针对上述问题,本文在贝叶斯建模的框架下,结合NLMM提出一种新的可靠性寿命数据的分析方法。

1 基本模型和假设

目前在寿命数据分析中,Weibull分布和对数正态分布是常用的两种对数位置-尺度族分布。然而,相关文献表明,可靠性工程师倾向于选择Weibull分布和最小极值分布来建模寿命数据。实际上,由于具有形状参数,Weibull分布的数据拟合能力远强于对数正态分布,其能够灵活地对多种类型的失效机制进行建模,因此本文选用Weibull分布来对寿命数据进行建模。对于表示寿命数据的随机变量T,双参数Weibull概率密度函数的常用参数形式:

(1)

累计分布函数:

(2)

式中,t>0是失效时间;β>0是Weibull分布的无量纲形状参数,不同的形状参数β值代表几种不同的失效机制;η>0是尺度参数,和随机变量T有相同的单位,这个尺度参数可以解释为失效时间分布的近似0.632 1分位数(也就是总体的63.21%失效的时间)。另外,寿命时间t的对数,即log(t)服从最小极值分布,它的位置参数为μ=logη,尺度参数为σ=1/β。

2 模型分析

2.1 NLMM分析

本文提出了一种将随机效应纳入失效时间模型的贝叶斯理论模型,假设随机效应是由于子抽样造成的[25],通常通过两种方式将随机效应纳入到模型中[25]:通过平均响应或通过模型参数。广义线性混合模型采用第一种方法,NLMM更灵活,通过模型参数来考虑随机效应。对于Weibull分布,一个普遍的假设是所有模型项都通过与对数尺度参数成线性关系输入,因此本文在贝叶斯理论的基础上,使用NLMM框架来整合子抽样导致的随机效应。在可靠性寿命试验中,有i=1,2,…,m个独立的试验单元,每个试验单元中有j=1,2,…,ni个子样本或者观测单元,可以将具有子抽样的Weibull分布的NLMM指定为

tij|μi,εi～Indep.Weib(βi,ηi)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

式中,假设有q个因子效应,xi是第i个试验设计的q×1的因子效应矩阵;ηi>0,βi>0分别是第i组的尺度参数和形状参数;tij>0是第i个试验单元上第j个样本的失效时间,服从尺度参数为ηi,形状参数为βi的独立Weibull分布;θ和γ是固定因子效应的回归系数矩阵;μi和εi分别为i=1,2，…,m时的随机效应;F(tij|βi,ηi,μi,εi)是第i个试验单元内数据的Weibull累计分布函数;f(μi)和f(εi)是随机效应μi和εi的正态概率密度函数;tip是第i个试验单元的p分位数(0

2.2 先验信息的选择

在使用极大似然估计(maximum likelihood estimation,MLE)方法时,删失数据的似然函数可能是无穷的,因此不存在MLE值,不能通过MLE来计算未知参数的估计值。规避可估计性问题的一种自然方法是使用贝叶斯方法[24]。通过使用合适的先验分布,可以计算出有代表性的后验分布,然后可以使用其均值或中位数来估计未知参数。

当尺度参数与形状参数均未知时,尚未有合适的联合先验分布形式。本文假设θ和γ的先验分布为q维多元正态分布π(θ|μ0,Σ0)和π(γ|μ1,Σ0),其中,

(10)

(11)

Σ0=M-1

(12)

M=XTX+A0

(13)

2.3 删失数据和失效时间数据的Gibbs采样算法

步骤 1设定(θ,γ,σμ,σε)初始值。

(14)

式中,yi=logti。

在此,采用吉布斯抽样方法来获得模型参数的估计值。首先,舍弃一些模型参数的初始抽样值(即燃烧期);然后对经过处理后的模型参数抽样值(即所获得模型参数的马尔可夫链)进行收敛性诊断。若上述模型参数的马尔可夫链围绕某个确定值在一定范围内上下波动,则可以利用该模型参数的马尔可夫链的均值作为未知参数的估计值。

需要特别指出的是,这里的“链”是指对模型参数采用贝叶斯方法估计时,通过MCMC模拟方法所获得的马尔可夫链,并使其平稳分布为待估计参数的后验分布。燃烧期是指采用MCMC或Gibbs抽样方法所获取的模型参数后验抽样值(即模型参数的马尔可夫链)在没有达到平稳分布时所采集到的样本量。

对式(2)和式(3),经过充分的迭代后所生成的链若通过收敛性诊断,则可以说明式(3)中所绘制的样本来自于(θ,γ,σμ,σε)的后验分布。此外,需要特别说明的是,Gibbs抽样是MCMC抽样方法中的一种抽样算法,适用于条件分布比边缘分布更容易采样的多变量分布。利用MCMC或Gibbs抽样方法所获得的模型参数的马尔可夫链,通常可以利用一些判断收敛性的统计量(如潜在尺度缩减因子 potential scale reduction factor, PSRF)或者一些收敛性诊断工具(如可视化的踪迹图、自相关图和密度曲线图)来推断模型参数的后验抽样值是否已经达到平稳状态,是否可能收敛于某个平稳分布[27-28]。

3 基于贝叶斯NLMM的可靠性数据分析

本文提出了一种将随机效应分别纳入尺度参数和形状参数的贝叶斯理论模型。在表1中,分别列出了文献[5],文献[10],文献[6],文献[2]中可靠性分析的模型以及本文所用到的NLMM模型。由于删失数据的存在,无法找到似然函数的共轭先验分布,后验分布没有简单的闭环形式。因此本文根据以上各参数的先验设定方法,通过MCMC方法动态模拟各参数后验分布的马尔可夫链,并根据各参数的后验样本进行收敛性诊断与参数估计。

表1 5种方法的模型假设

所提方法的基本流程如图1所示,具体步骤如下。

图1 所提方法的可靠性分析流程图

步骤 1首先需要检验收集的数据是否服从Weibull分布,在服从Weibull分布的前提下验证形状参数是否是恒定常数。

步骤 2从最初的试验数据中利用最小二乘法计算尺度参数和形状参数的估计值,并根据式(10)～式(13)设置合理的先验分布。

步骤 3根据步骤2的先验分布和未知参数的似然信息就可以利用Gibbs采样迭代地估计参数值。若经过K次迭代后生成的链已经符合收敛性检验,则可以将去掉燃烧期后的链的均值作为模型参数的估计值。若经过K次迭代后生成的链没有收敛,则需要重新设定初始值,重复上述迭代过程。

步骤 4根据步骤3的后验估计值得出分位数的估计值和置信区间。

步骤 5计算均方根误差(root mean square error,RMSE)和相对偏差(relative bias,RB)[29],并且对几种方法进行对比及验证。

4 实例分析

4.1 实例背景

该实例来自文献[9],主要研究电容器的可靠性问题,响应变量为电容器失效的时间,影响可靠性寿命的因子主要包括:温度和电压。其中温度应力为2水平变量,电压应力为4水平变量,实验设置由Zelen描述为“N个组件同时放置在试验单元上”,这是一个子抽样试验设计。该试验一共有8个试验单元,每个单元有8个电容器,并且采用定数截尾的方式来终止试验,本试验预先设定的失效数目为4。本实验中的观测单位是每个试验台上的8个玻璃电容器,每个试验台均采用单独的温度、电压处理组合。试验的设计和寿命数据如表2所示，*代表删失数据。

4.2 数据分析

运用本文提出的方法重新分析这个例子。首先,检验收集到的数据是否遵循Weibull分布。用Minitab软件使用AD(Anderson-Darling)统计量检验所收集到的数据是否遵循Weibull分布,其结果如图2所示。

表2 试验设计及寿命数据

图2 Minitab输出结果

从图2中发现,每个处理组合的AD值都比较小,这意味着寿命数据很好地遵循Weibull分布。其次,利用Minitab软件检验Weibull分布的形状参数是否为常数,结果表明,χ2=19.507 5,自由度df为7,p=0.007,这意味着可在2α=0.05水平下拒绝原假设。另外,可以看出最大的形状参数(26.991)和最小的形状参数(2.153 2)相差较大,不能简单地将形状参数假设为一个恒定常数。

根据上述数据分析的结果,对模型参数设置合理的先验分布,然后运用MCMC方法对模型参数进行估计。在此,对模型参数进行了120 000次迭代抽样,并舍弃了前20 000次的抽样,获得了100 000次的模型参数抽样值。为了消除未知参数抽样值之间的自相关性,在此对抽样获得的后验样本每间隔10步抽样一次,最后得到10 000个有效的模型参数抽样值。然后,对所获得的模型参数抽样值进行收敛性检验,并列出踪迹图和自相关图来帮助判断参数抽样值的收敛性。考虑到篇幅限制,这里只提供θ2和第3组的第一分位数对应的参数后验值的踪迹图(见图3)和自相关图(见图4)。

从图3和图4中可知,在抽样过程中,由于参数的不确定性其后验抽样值会围绕某个确定的值上下波动,并且呈现出稳态分布的特征,因此可以直观地判断出所有的模型参数估计值具有良好的收敛性,可以利用其参数的抽样值进行后续的统计推断和数据分析。

图3 模型参数后验值的踪迹图

图4 模型参数后验值的自相关图

表3列出了表1描述的模型的参数估计值,-为模型中未估计的参数。在应用本文提出的方法时,在考虑随机效应对尺度参数和形状参数都有影响时,产生的结果与之前的模型结果相比稍有不同。

表3 模型Ⅰ～Ⅴ的参数估计

从表3中可以看出,尺度参数的回归模型系数没有太大差别,但是形状参数的回归系数与前4种模型有较大差别,主要是因为本文方法考虑了随机效应对形状参数的影响。另外,形状参数的随机效应(σε=0.400 8)大于尺度参数的随机效应(σμ=0.262 7),因此仅仅考虑随机效应对尺度参数的影响而忽略随机效应对形状参数的影响是不合理的。在模型IV中,σμ=0.184,小于模型V中的0.262 7,这意味着,如果忽略了加速寿命分析中的随机效应对形状参数的影响,那么尺度参数的随机效应就会被低估。为了获得低分位数的精确估计,因此有必要将随机效应同时纳入到尺度参数和形状参数中。表4给出了5种模型尺度参数和形状参数的估计值。

表4 模型Ⅰ～Ⅴ的尺度参数和形状参数估计

在可靠性优化中,产品在其低分位点处的失效时间往往是生产商关注的焦点[29]。因此分位数估计在可靠性分析,尤其是在加速寿命分析中尤为重要。表5总结了模型Ⅰ～Ⅴ的第1、第5、第10和第50分位数及其95%的置信区间。由于模型Ⅲ、模型Ⅳ和模型Ⅴ考虑了非恒定形状参数的假设,3个模型的第1、第5和第10分位数的结果与模型Ⅰ和模型Ⅱ中的稍有不同。 另外,模型Ⅴ和模型Ⅲ、模型Ⅳ的分位数估计以及置信区间有较大不同,产生这个结果的原因是模型Ⅴ考虑了随机效应对形状参数的影响。

从表5可以看出,模型V的95%置信区间范围大于前4个模型的范围,虽然本文方法在一定程度上增加了分位数估计的不确定性,但是改进了未知参数的估计精度。另外,表3的结果表明不能为了消除这种不确定性而忽略随机效应对形状参数的影响。

表5 模型Ⅰ～Ⅴ的分位数估计以及置信区间

续表5

为了说明本文所提方法的有效性,比较了所提方法与其他方法对分位数的估计精度。为了确保评价估计方法的精度,本文采用RB和RMSE两个指标[30]。

(15)

(16)

RB1、RB2、RB3、RB4分别为文献[5],文献[10],文献[2]以及本文方法的RB;RMSE1、RMSE2、RMSE3、RMSE4分别代表文献[5],文献[10],文献[2]以及本文方法的RMSE。

RB和RMSE的分析比较结果分别列于图5和图6中。

从图5中可以看出,RB4在所有方法中是最小的,这意味着本文提出的方法精度是比较好的。在图6中可知,文献[10]以及本文所提方法所获得RMSE的值小于文献[5]以及文献[2]中的值。在第1、第5、第10分位数中,虽然RB3的值是较小的,但是RMSE3的值是大于其他3种方法的,由此表明虽然文献[2]所提方法正确地考虑了随机效应对尺度参数的影响,但是却增大了分位数预测值的波动。由于本文所提方法考虑了随机效应对形状参数的影响,同时利用贝叶斯随机抽样方法考虑了模型参数的不确定性以及随机误差对模型参数估计的影响,因此在对低分位数的估计方面是更为稳健。

图5 分位数估计的RB

图6 分位数估计的RMSE

5 方法讨论

针对可靠性寿命的低分位数估计的偏差和存在删失情况的问题,为了改善分析结果的可靠性,本文研究加速寿命试验数据的分析与处理方法,结合贝叶斯抽样方法与NLMM,提出一种改进方法。首先,利用Weibull分布拟合寿命试验数据,并检验各组数据的形状参数是否恒定。其次,考虑随机效应对尺度参数和形状参数的影响,建立相应的贝叶斯模型。此外,结合电容器试验案例,对不同优化方法的研究结果进行了比较。根据本文的研究结果可得结论如下。

(1) 在ALT中,因研究对象的复杂程度、试验因子数量[31]、试验条件存在很大差异,加上样本数量有限、存在混合删失现象等原因,使得试验数据的分析、建模存在较大难度。同时,因为不存在通用的模型和方法,需要针对实际的样本数据开展针对性分析和研究,提出合适的模型和算法。由于本文选用的案例数据属于II型删失,因此本文所提方法只是初步验证了算法在II型删失中的有效性,但还不足以证明该方法的通用性和适用范围,例如在I型删失数据中,每组删失个数不同,很难建立一个通用的模型来拟合数据。

(2) 本文在以往相关文献的基础上进一步考虑了尺度参数的随机效应。以往文献通常假设非随机化设计只影响尺度参数,然而,尺度参数和形状参数通常是在建模过程中结合具体的试验数据共同拟合估计出来。因此不能在建模过程中忽视形状参数的随机效应对研究结果的影响。正如文献[2]所指出的那样“在未来研究中应该考虑将尺度参数和形状参数的随机效应纳入统一模型框架中。然而,由于NLMM的数值积分问题非常复杂,因此有待进一步研究。” 本文利用贝叶斯方法对所构建的新模型进行参数估计,不仅在以往研究基础上考虑形状参数的随机效应,而且还利用贝叶斯抽样方法考虑模型参数不确定性以及随机误差对模型参数估计结果的影响。根据研究结果可知,在建模过程中考虑形状参数的随机效应会在一定程度上增加分位数估计结果的不确定性,但同时又能提高低分位数估计结果的稳健性。然而,表3中的研究结果表明:形状参数的随机效应σε大于尺度参数σμ的随机效应。因此,在建模过程中忽视形状参数的随机效应可能会导致一些与实际不太吻合的研究结论,甚至是错误结果。

(3) 若直接用数值积分方法来计算上述的NLMM往往会非常复杂,因为其模型参数的后验分布往往异常复杂,在很多情况下往往无法获得其封闭的概率密度函数。因此,借助贝叶斯抽样方法(如Gibbs抽样)来获得模型参数的后验抽样值,然后据此来计算模型参数后验抽样值的均值,从而借助贝叶斯方法实现对NLMM的估计,解决了NLMM模型的计算精度和积分点选取的问题。

6 结 论

在加速寿命设计的可靠性分析中,往往需要考虑未完全随机化设计对可靠性寿命的影响。本文结合贝叶斯抽样方法与NLMM方法提出了一种可靠性寿命分析的新方法。该方法不仅考虑了随机效应对模型尺度参数和形状参数的影响,而且还解决了两阶段方法中似然函数不统一和NLMM中使用Gauss-Hermite法中计算精度和积分点选取的问题。另外,通过实例证明该方法在存在高波动情况下对低分位数估计也是稳健的。

需要特别指出的是,本文中先验分布的选择具有一定的主观性,下一步需要进一步地验证先验信息对可靠性寿命估计的灵敏度。此外,在可靠性分析时,尤其对于ALT来说,在低应力水平下样本失效数据少,甚至出现一组数据全部删失的情况。若某组试验数据全部删失,此时可利用的信息较少,无论样本的失效机制是否相同,都不能计算该组样品的特征参数。文献[32]针对失效数据少的问题,提出了一种分层贝叶斯方法,从其他分组中共享信息。在处理某组样本数据全部是删失情况时,文献[33]只是简单的删除此组数据,可能会剔除掉一些重要信息。因此,在某组寿命数据全部删失的情况下,未来需要考虑的是如何借助贝叶斯方法获取更多的信息以便进行下一步分析。另外,本文仅考虑试验数据服从威布尔分布的情况。若试验数据不服从Weibull分布时,则需要进一步判断试验数据是否满足常见分布如对数正态分布或指数分布。若试验数据满足常见分布,则可以参考本文所提出的贝叶斯建模方法对相关模型进行参数估计与可靠性数据分析。若试验数据不满足常见分布,本文所提方法将不再适用,则有待未来对此类问题开展更为深入的研究。