产品的失效过程是一个复杂的过程，从失效形式上来看，可分为突发型失效与退化型失效两种[1]．如果一个产品同时具有这两种失效模式，即产品在性能退化过程中也可能发生突发失效，那么该产品的失效是这两种失效模式间竞争的结果，很多产品具有竞争失效模式，例如：对电阻而言，其失效可能是高压击穿、开路或阻值漂移中的一种，前两种失效模式是突发失效，阻值漂移则属于退化失效[2-3]．与考虑单一失效的模型相比，竞争失效模型更能准确刻画产品的失效，更加灵活精确地表征部件性能变化对系统性能和可靠性的影响[4]．文献[5]基于维纳过程构建了竞争失效可靠性模型，运用极大似然估计(MLE)给出参数的估计．文献[6]基于伽马过程构建了竞争失效可靠性模型，分析了相关参数对系统可靠性的影响．目前，基于竞争失效模型的研究已取得很多很好的成果，但是仍存在两方面的问题：一是建立的退化模型主要是维纳过程、伽马过程、逆高斯过程等特定的退化模型，缺少一种适用范围广泛的退化模型；二是对于竞争失效模型的参数估计主要基于频率方法和主观贝叶斯方法，这两种方法在大样本、有历史信息的情况下能取得很好的效果，但针对小样本、无历史信息的情况往往效果不佳[7-8]．针对问题一，目前常用的退化模型主要为维纳过程、伽马过程和逆高斯过程等[9-10]．文献[11]提出指数扩散(ED)过程也可用于刻画产品的性能退化规律，常用三种退化过程都是其特例，相比之下ED过程适用范围更广，效果更好．Tweedie指数扩散(TED)模型是一类重要的ED模型，文献[12]针对TED过程进行了加速退化试验和带有测量误差的模型分析，并通过实例证明TED模型在产品性能退化研究中的有效性．鉴于此，本研究引入TED过程刻画产品的退化失效过程．针对问题二，目前关于竞争失效模型的参数估计大部分是基于频率方法和主观贝叶斯方法．对于小样本情况且没有历史信息借鉴的情况，基于频率方法和主观贝叶斯方法所获得的参数估计效果往往很一般[13]．然而，客观贝叶斯估计在无历史信息的小样本情形下能有相对更好的参数估计效果[14]．Jeffreys先验及Reference先验是两种常用的客观先验，也具有一些较好的性质，如参数一对一变换下的不变性等．鉴于此，本研究选取客观贝叶斯方法对竞争失效模型参数进行估计．1 模型介绍1.1　TED过程退化失效模型令Y(t)表示性能退化量，t表示时间，满足以下条件的随机过程{Y(t),t≥0}称之为ED过程：a．Y(0)=0；b．{Y(t),t≥0}具有平稳独立增量，即对于0≤t1≤t2≤t3≤t4，Y(t2)-Y(t1)和Y(t4)-Y(t3)是独立的；c．{Y(t),t≥0}的增量服从ED分布，即Y(t2)-Y(t1)~FED(μ(t2-t1),λ)，(1)式中FED(μt,λ)为ED分布．FED概率密度函数为f(y;μ,λ)=c(y|λ,t)exp{λ[yθ(μ)-th(θ(μ))]}, (2)式中：λ为扩散参数；c(⋅)为正则化函数，保证式(2)的积分为1；h(⋅)为累积函数，当λ=1时h(⋅)导数表示ED分布的连续累积量；μ为漂移参数，满足μ=h'(θ)τ(θ)，h'(θ)为h(θ)关于θ的一阶导数，θ=τ-1(μ)，τ-1(⋅)为τ(⋅)的反函数．根据式(1)和(2)，可以得到ED过程的均值和方差分别为E(Y(t))=μt；Var(Y(t))=V(μ)t/λ，式中：Var(⋅)为方差；V(μ)=h″(τ-1(μ))=h″(θ)为方差函数，其中h″(θ)为h(θ)关于θ的二阶导数．TED模型的方差函数与均值间是幂函数关系，有V(μ)=μp,p∈(-∞,0]⋃[1,∞)，其中p为类别参数，不同的p值表征了不同的退化模型．对TED模型，除一些特殊情况，函数c(⋅)没有具体表达式．根据现有研究[9]，鞍点近似方法为TED模型提供了高精度的近似，可表示为：f(y;μ,λ)≅λ/(2πt1-pyp)∙exp(-λtd(y;μ)/2); (3)d(y;μ)=2ytlnyμt-(yt-μ)(p=1),2lnμty+yμt-1(p=2),     2y2-ptp-2(1-p)(1-p)-yμ1-p(1-p)t+μ2-p2-p(p≠1,2). (4)设ω为产品的失效阈值，TED为Y(t)首次达到失效阈值ω的时刻，则此时产品发生失效，有TED=inft|Y(t)≥ω,t0，式中inf{⋅}为集合的下确界．1.2　突发失效模型突发失效模型的失效时间TW服从韦布尔分布，则TW的概率密度函数为fW(t)=(βtβ-1/kβ)exp[-(t/k)β]，(5)TW的累积分布函数为FW(t)=1-exp[-(t/k)β]，(6)式中：k为尺度参数；β为形状参数．因此产品的失效时间为T=min(TED,TW)．1.3　极大似然估计假定从一批产品中抽取n个进行性能退化试验，记Y(tij)为第i个产品在时刻tij的性能退化量，i=1,2,…,n,j=1,2,…,mi，增量Δyij=Y(ti,j)-Y(ti,j-1)，相互独立且服从FED(μΔtij,λ)，其中Y(ti,0)=0,Δtij=ti,j-ti,j-1，而第i(i=1,2,…,n)个产品的突发失效时间记为Ti1．令Θ=(μ,λ,k,β)．基于观测数据，其对数似然函数可表示为logL(Θ)=∑i=1n∑j=1mi12[logλ2π-(1-p)logΔtij-plogΔyij-λΔtijd]+∑i=1n[log β+(β-1)logTi1-βlogk-(Ti1/k)β]. (7)在式(7)中对韦布尔分布的参数(k,β)求一阶偏导，即可获得MLE值(k^M,β^M)．对TED过程的参数(μ,λ)求一阶偏导，并令其为零，可得μ^M=B/C,λ^M=(1-p)(2-p)N/(2U)，(8)式中：N=∑i=1nmi；U=A-B2-pCp-1；A=∑i=1n∑j=1mi(Δyij2-pΔtijp-1)；B=∑i=1nyimi；C=∑i=1ntimi．再将(k^M,β^M,μ^M,λ^M)代入式(7)中，求解偏对数似然函数的最大值，即可获得p的MLE值p^M．2 客观先验先验信息是贝叶斯估计的关键要素之一，选取合理的先验信息对评估结果有重要的影响．但某些产品样本量很少，且很难获取先验信息，因此无信息先验的问题不可避免存在．目前，常用的两类无信息先验分布是Jeffreys先验和Reference先验．2.1　Fisher信息矩阵由于Jeffreys先验及Reference先验的推导都要用到Fisher信息矩阵，这里先给出Fisher信息矩阵．对于TED模型，根据文献[15]，当p值变化时，Fisher信息矩阵可能是奇异的，从而无法计算出其逆矩阵．另外，由于TED模型关于p的一阶偏导及二阶偏导太复杂，因此无法获得相应的期望值，从而导致后续矩阵无法获得显示表达式．鉴于此，对于TED模型只考虑参数λ和μ的客观先验．定理1 基于TED过程与韦布尔分布的竞争失效模型，参数的Fisher信息矩阵为I(Θ)=diag(I(λ,μ),I(k,β))，(9)式中：I(λ,μ)=λ-2/200Cλμ-p；I(k,β)=nβ2k-2-(1+γ1)/k-(1+γ1)/k(1+2γ1+γ2)/β2．证明 对似然函数(7)关于参数Θ求二阶偏导数，再添加负号，求期望，即可得到Fisher信息矩阵．这里不再给出详细求导步骤．2.2　Jeffreys先验及Reference先验Jeffreys先验是基于Fisher信息矩阵行列式的平方根所求得的[14]，该先验由于其在参数一对一变换下的不变性而被广泛使用．定理2 基于竞争失效模型，参数Θ的Jeffreys先验πJ(Θ)表示为πJ(Θ)∝λ-1/2μ-p/2k-1．(10)证明 由Fisher信息矩阵(9)，可得I(Θ)=I(λ,μ)I(k,β)=Cn2(γ2-γ12)/(2λμpk2).再根据Jeffreys先验的定义，定理即可得证．Reference先验通过最大化后验分布和先验分布间的相对熵得到[16]．定理3 关于序为{μ,λ,k,β}和{k,β,μ，λ}的Reference先验为πR(Θ)∝λ-1μ-p/2k-1β-1．(11)证明 根据文献[16]给出的多参数Reference先验计算公式及式(9)可得S=I-1(Θ)=diag(|I-1(μ,λ)|,|I-1(k,β)|)．令：W1=2λ2;W2=|I-1(μ,λ)|；W3=diag(|I-1(μ,λ)|,β-2Zk2)；W4=diag(|I-1(μ,λ)|,|I-1(k,β)|)，式中Z=(1+2γ1+γ2)/[n(γ2-γ12)]．令Hi=Wi-1(i=1，2，3，4)，hi为矩阵Hi的右下第i行第i列元素．对Θ选择紧集Ωi=[a1i,b1i][a2i,b2i][a3i,b3i]⋅[a4i,b4i]，并令limi→∞a2x+1i=0，limi→∞b2x+1i=0,x=0,1，可得π4i(β|μ,λ,k)=R4il[a4i,b4i]/β，其中lΩ表示集合Ω的指示函数．可知πji(Θ[∼(j-1)]|Θ[j-1])=πj+1i(Θ[∼(j)]|Θ[j])exp12Eji(log(hj)|Θ[j])l[aji,bji]∫ajibjiexp12Eji(log(hj)|Θ[j])dΘj,式中：j=1,2,3；Eji(log(hj)|Θ[j])=∫log(hj)πj+1i(Θ[∼j]|Θ[j])dΘ[∼j]；Θ[j]=(Θ1,Θ2,…,Θj)；Θ[∼j]=(Θj+1,Θj+2,…,Θ4)．由此可得π3i(k,β|μ,λ)的表达式，以此类推得π1i(Θ)∝λ-1μ-p/2k-1β-1R1iR2iR3iR4il[a1i,b1i]l[a2i,b2i]l[a3i,b3i]l[a4i,b4i]，式中Rji=(logbji-logaji)-1为不依赖于Θ的常数．最终得πR(Θ)=limi→∞[π1i(Θ)/π1i(1,1,⋯,1)]=λ-1μ-p2k-1β-1．类似地，在基于序为k,β,μ,λ的情况下，可证明Reference先验也为πR(Θ)．2.3　概率匹配先验基于概率匹配先验得到的参数的后验分布具有良好的频率性质，但值得注意的是：绝大多数情况下参数的精确概率匹配先验不易求得，故一般是研究渐进概率匹配先验．假设参数向量θ=(θ1,θ2)，其中θ1为感兴趣的参数，X1,X2,...,Xn是来自分布密度函数为f(x|θ)的样本．若关于样本分布f(x|θ)有Px|θ(θ1≤θ1α)=α+O(n-i/2)，(12)式中θ1α为θ1的后验分布的α分位数．若Pθ|x(θ1≤θ1α|X1,X2,...,Xn)=α，则称π(θ)是i阶(渐近)概率匹配先验；若Px|θ(θ1≤θ1α)=α，则称π(θ)是i阶精确概率匹配先验．为分析Jeffreys先验和Reference先验是否为概率匹配先验，先给出须满足的两个方程[17]．a．Fisher信息矩阵为I(θ1,θ2)=I11I12I21I22，基于θ1的概率匹配先验π(θ1,θ2)是以下微分方程的解，即∂∂θ2I12π(θ1,θ2)I22W-∂∂θ1π(θ1,θ2)W=0，(13)式中W=I11-I122/I22．b．假如θ1和θ2正交，即I12=I21=0，则上述微分方程可化简为∂∂θ1π(θ1,θ2)W=0．(14)化简式(14)，可得参数θ1的概率匹配先验为π(θ1,θ2)=g(θ2)I11，其中g(θ2)为θ2的一个连续函数．根据以上结论，可以给出定理4．定理4 Jefffreys先验πJ(Θ)是关于参数μ的概率匹配先验，但不是λ,k,β的概率匹配先验；Reference先验πR(Θ)是关于参数λ,μ,k,β的联合概率匹配先验．证明 根据式(13)和(14)，可得出关于Jeffreys先验πJ(Θ)的概率匹配先验须满足微分方程：∂∂λ2λ2πJ=∂∂λ2λμ-pk≠0;∂∂μμpcλπJ=∂∂μ1Cλk=0;    ∂∂kkβ1+2γ1+γ2n(γ2-γ12)πJ+∂∂ββ(1+γ1)1+2γ1+γ2⋅1+2γ1+γ2n(γ2-γ12)πJ≠0;    ∂∂kk(1+γ1)βn(γ2-γ12)πJ+∂∂ββ1n(γ2-γ12)πJ≠0.可见，对于Jeffreys先验πJ(Θ)，只有参数μ满足微分方程，所以πJ(Θ)是关于参数μ的概率匹配先验，并不是参数λ,k,β的概率匹配先验．同理，Reference先验πR(Θ)是λ,μ,k,β的概率匹配先验．3 后验分析3.1　后验分布性质基于Jeffreys无信息先验分布，推导相应的后验分布．根据贝叶斯公式得参数的联合后验分布为πJ(Θ)∝λ-12μ-p2∏i=1n∏j=1miλ2πΔtij1-pΔyijp⋅exp-λΔtijd(Δyij;μ)21k∏i=1nβTi1β-1kβexp-Ti1kβ. (15)定理5 基于Jeffreys先验πJ(Θ)的联合后验分布πJ(Θ)是正常的．证明 令：R1=∫0+∞∫0+∞1λ12μp2∏i=1n∏j=1miλ2πΔtij1-pΔyijp⋅exp(-λΔtijd/2)dλdμ;R2=∫0+∞∫0+∞1k∏i=1nβTi1β-1kβexp-Ti1kβdkdβ，则根据式(15)，后验分布πJ(Θ)具有如下关系：∫ΘπJ(Θ)dΘ∝R1R2．所以，要证明后验分布πJ(Θ)是正常的，只须证明R1和R2是有界的．对于R1，有R1∝μ-p2λN+12-1exp-λ2∑i=1n∑j=1miΔtijd(Δyij;μ)=∫0BCμ-pΓ[(N+1)/2]A(1-p)(2-p)-Bμ1-p1-p+Cμ2-p2-p(N+1)/2dμ+∫BC+∞μ-pΓ[(N+1)/2]A(1-p)(2-p)-Bμ1-p1-p+Cμ2-p2-p(N+1)/2dμ.当p∈(-∞,0⋃(1,2)时，有μ-pA(1-p)(2-p)-Bμ1-p1-p+Cμ2-p2-pN+12=Ο[μp-2N2-1](μ→+∞);当p∈(2,+∞)时，有μpn/2Γ[(N+1)/2]Aμp(1-p)(2-p)-Bμ1-p+Cμ22-pN+12=Ο(μ-p/2)(μ→+∞).故可证R1+∞．对于R2，有R2Γ(n)/∏i=1nTi1∫0+∞βn-1Q-βdβ+∞，式中Q=∏i=1n(TM/Ti1)1，其中TM为(T11,T21,...，Tn1)中最大值．因此R1和R2均有界，则证明πJ(Θ)是正常的．定理6 基于Reference先验πR(Θ)的后验分布πR(Θ)是正常的．证明 证明过程同定理5类似．在Jeffreys先验πJ(Θ)下，后验分布为πJ(Θ)，为计算方便，对后验分布化简，令δ=kβ，可得参数λ,μ,δ,β的满条件边缘后验分布有：πJ(λ|μ,δ,β)∝GaN+12,A(1-p)(2-p)-Bμ1-p1-p+Cμ2-p2-p;    πJ(μ|λ,δ,β)∝μ-p/2exp-λA(1-p)(2-p)-Bμ1-p1-p+Cμ2-p2-p;    πJ(δ|λ,μ,β)∝δ-(n+1)exp-∑i=1nTi1β/δ~IGan,∑i=1nTi1β;πJ(β|λ,μ,δ)∝βnk-nβ∏i=1nTi1βexp-k-β∑i=1nTi1β.在Reference先验πR(Θ)下，参数λ,μ,δ,β的满条件后验分布有：    πR(λ|μ,δ,β)∝GaN2,A(1-p)(2-p)-Bμ1-p1-p+Cμ2-p2-p;    πR(μ|λ,δ,β)∝μ-p/2exp-λA(1-p)(2-p)-Bμ1-p1-p+Cμ2-p2-p;    πR(δ|λ,μ,β)∝δ-(n+1)exp-∑i=1nTi1β/δ~IGan,∑i=1nTi1β;πR(β|λ,μ,δ)∝βn-1k-nβ∏i=1nTi1βexp-k-β∑i=1nTi1β.以上的参数λ的满条件分布是伽马分布，记为Ga(⋅)，参数δ的满条件分布是逆伽马分布，记为IGa(⋅)，都可以轻松获得其后验分布的均值或众数等，即可获得参数λ,δ的客观贝叶斯估计，再对δ化简，可得λ,k的客观贝叶斯估计．但对于μ及β，并不是常规的分布函数，很难获得后验均值或众数．那么就可以考虑下面的Gibbs抽样算法．3.2　后验分布的Gibbs抽样算法基于Jeffreys先验下Gibbs抽样算法步骤如下．步骤1 选择初始值Θ(0)=(λ(0),μ(0),k(0),β(0))，δ(0)=(k(0))β(0)，并令i=1．步骤2 对于Jeffreys先验，从伽马分布中提取λ(i)，并从逆伽马分布中提取δ(i)．步骤3 使用Metropolis-Hastings(MH)算法从πJ(μ(i-1)|λ(i),δ(i),β(i-1))生成μ(i)，提议分布选择N(μ(i-1),σμ2)，具体操作如下：a．从提议分布N(μ(i-1),σμ2)中产生提议样本μ*；b．从标准正态分布U(0,1)中产生随机样本Uμ；c．计算接受概率rμ=min1,πJ(μ*|λ(i),δ(i),β(i-1))πJ(μ(i-1)|λ(i),δ(i),β(i-1))；d．根据μ(i)=μ*(Uμ≤rμ)μ(i-1)(Uμrμ)，获取样本μ(i)；e．令i=i+1．步骤4 用MH算法生成β(i)，同步骤3．步骤5 重复D次步骤2~4，得到λ(i),μ(i),δ(i)，β(i),将其化简可得(λ(i),μ(i),k(i),β(i))，其中i=1,2,…,D．步骤6 为保证收敛性并消除对初始值选择的影响，舍弃前D0个抽取的样本，剩余的D-D0个样本为λ(i),μ(i),k(i),β(i)(i=D0+1,D0+2,…,D)，基于剩余样本对参数进行贝叶斯估计．因此，基于Jeffreys先验，λ,μ,k和β的客观贝叶斯估计为：λ^BJ=1D-D0∑i=D0+1Dλ(i);μ^BJ=1D-D0∑i=D0+1Dμ(i);k^BJ=1D-D0∑i=D0+1Dk(i);β^BJ=1D-D0∑i=D0+1Dβ(i).类似可以获得Reference先验下的贝叶斯估计．4 数值模拟本研究运用模拟仿真方法来验证客观贝叶斯估计方法的有效性．随机选取参数p=2.5,λ=5,μ=1．初始时刻Y(0)=0，每个产品检测次数m=30，检测间隔时间Δt=1，失效阈值ω=40．突发失效模型假定服从于k=10，β=3的韦布尔分布．为了对比分析，将客观贝叶斯估计方法与MLE及主观先验下的贝叶斯估计进行比较．选取TED模型的参数λ,μ的先验分布分别是正态分布和逆伽马分布．韦布尔分布参数k和β先验分布分别为正态分布和伽马分布．选取两组先验分布参数值，分别如下：第一个主观先验πX1为λ~N(5，0.32)，u~IGa(3,2)，k~N(10,22)，β~Ga(3,1)，这个先验分布的均值即为参数的真值；第二个主观先验πX2为λ~N(4,0.42)，u~IGa(2,2)，k~N(8，2.32)，β~Ga(4,1)，这个先验分布的均值不等于参数的真值．4.1　收敛速度比较对TED退化过程竞争失效模型进行参数估计时，由于公式复杂，须要进行迭代计算．下面对客观贝叶斯估计、MLE及主观先验下的贝叶斯方法进行比较．设θ为待估参数，θ^i为θ的第i步迭代的估计值，给定误差精度ε=1×10-3，令θ^i+1-θ^iε来衡量迭代的收敛性，并通过Nnum{θ^i+1-θ^iε}/N来计算收敛速度，其中Nnum{⋅}为迭代满足θ^i+1-θ^iε的总次数，总体比值越高，代表收敛效果越好．基于上述假设，模拟不同样本下的竞争失效数据，使用Gibbs抽样来计算不同参数估计方式下的收敛比值，每次模拟迭代次数为8 000．结果如表1所示．10.13245/j.hust.239299.T001表1基于不同估计下的收敛速度比较n参数πJπRMLEπX1πX220λ0.6790.6840.6290.6250.612μ0.6740.6900.5320.5280.524k0.7630.7800.6950.6880.676β0.7910.7900.5320.5300.51830λ0.7320.7340.7280.7160.708μ0.7180.7200.6080.6020.596k0.8350.8230.7930.7600.772β0.8710.8780.6480.6650.65450λ0.7500.7520.7460.7380.736μ0.8460.8520.8440.8400.832k0.9270.9350.8890.8810.879β0.9240.9260.8550.8620.858由表1可知：对于相同的样本量，在大多数情况下，客观贝叶斯估计的收敛速度比MLE和主观贝叶斯估计的收敛速度快．特别是在小样本情况下，客观贝叶斯估计的收敛速度比其他方法好得多．其中Reference先验下的贝叶斯估计收敛速度最快，主观先验分布πX2下的收敛速度最慢，这可能是由于该先验分布的均值并不等于参数真值．基于主观先验分布πX1的贝叶斯估计的收敛速度介于MLE及主观先验分布πX2之间，虽然主观先验分布πX1均值等于参数真值，但是先验分布函数的形式也具有信息量，而先验分布πX1的选取具有主观性，可能没有完全正确刻画参数的先验信息，从而导致收敛速度比MLE慢，但是比主观先验分布πX2下的收敛速度更快．因此，如果在并无先验信息的情况下对参数进行贝叶斯估计，建议采用无先验信息的客观贝叶斯估计，尤其是可以考虑Reference先验．随着样本量增加，不同估计方法的收敛速度越来越接近，即参数估计迭代的收敛速度与样本量的增加成正比，这是因为随着样本量增加，提供的信息量越多，能尽早迭代到满足精度要求的解．4.2　均值与均方误差比较为进一步对不同的估计方法进行比较，基于参数真值p=2.5,λ=5,μ=1,k=10,β=3随机模拟产生一组TED退化数据及韦布尔分布突发失效数据，并运用客观贝叶斯方法、极大似然方法及主观贝叶斯方法给出参数估计，重复此过程3 000次，获得参数估计的均值及均方误差，所得结果列于表2中．10.13245/j.hust.239299.T002表2不同估计方法下的均值和均方误差nπJπRMLEπX1πX2均值均方误差均值均方误差均值均方误差均值均方误差均值均方误差20λ5.022 76.67×10-25.016 56.43×10-25.054 46.78×10-24.597 30.615 04.567 30.663 0μ0.999 86.12×10-60.999 86.01×10-60.999 76.13×10-60.991 74.13×10-50.970 58.26×10-5k10.087 20.470 09.951 50.417 09.914 40.498 09.907 10.703 09.902 60.886 0β3.001 81.75×10-42.999 31.56×10-43.099 54.49×10-23.081 74.25×10-23.108 27.95×10-230λ5.027 77.20×10-25.016 06.39×10-25.028 16.17×10-24.798 20.498 04.761 00.654 0μ0.999 73.59×10-60.999 83.52×10-60.999 73.56×10-60.989 52.26×10-50.972 83.35×10-5k10.052 40.309 010.031 10.305 010.082 50.256 010.102 40.606 010.120 70.839 0β3.001 01.58×10-42.999 51.26×10-43.098 84.47×10-22.964 73.51×10-22.959 14.36×10-250λ4.983 85.61×10-24.987 04.57×10-64.983 65.62×10-24.833 94.99×10-24.825 86.15×10-2μ1.000 02.62×10-61.000 02.58×10-61.000 02.64×10-60.991 21.12×10-50.981 92.27×10-5k10.018 00.225 010.025 10.244 09.988 80.219 010.051 50.327 010.072 60.517 0β2.999 11.29×10-42.999 51.07×10-23.076 41.89×10-22.983 42.49×10-22.970 23.72×10-2从表2可以看出：a．随着样本量增加，无信息先验下的客观贝叶斯估计、MLE及主观先验下的贝叶斯估计越来越接近真实值，均方误差也越来越小，且当样本量较大时，不同估计的均值及均方误差整体上越来越接近；b．当样本量相同时，πR的客观贝叶斯估计效果最好，πX2下的贝叶斯估计效果最差，所以当进行参数估计时，如果无任何先验信息，建议采用Reference先验下的客观贝叶斯估计方法；c．对于客观贝叶斯估计，样本量n=50的估计值和均方误差略优于n=20下的相应值，但这种差异性整体并不是很大，所以当样本量较小时，客观贝叶斯估计的效果也很好，估计比较稳健．5 实例分析将所提竞争失效模型及相应的客观贝叶斯方法应用于文献[18]记录的GaAs激光器实例数据中进行分析，验证模型及方法的有效性．共有22个样本试验，每隔250 h测试一次，至4 000 h为止．假设产品的失效阈值为10，其中15个受试产品性能逐渐退化，另外7个产品发生突发失效．5.1　退化失效模型与突发失效模型的选择根据数据可知：突发失效伴随产品整个退化过程，失效时间无规律性，故退化程度与突发失效无关．先对15个退化产品进行分析，分别用TED、维纳、伽马和逆高斯过程对其建模，用客观贝叶斯估计给出参数估计，并计算AIC(赤池信息量准则)值和最大对数似然函数值，见表3．AIC值为VAIC=-2ln L^+2k，其中：ln L^为最大对数似然函数值；k为未知参数个数．具有最大对数似然函数值和最小AIC值的模型被认为是拟合效果最好的模型．10.13245/j.hust.239299.T003表3不同模型下GaAs激光器数据的AIC和对数似然对比模型TED维纳伽马逆高斯p2.850.002.003.00λ1.40×10-46.24×10-30.035.44×10-5μ2.06×10-32.02×10-32.03×10-32.03×10-3VAIC-144.345 0-85.135 4-130.222 8-142.211 6lnL̂75.172 544.567 767.111 473.105 8由表3可看出：当p=2.85时，TED模型具有最大的对数似然函数值和最小的AIC值，因此TED模型可以看作是GaAs激光数据的最佳退化模型．对于突发失效模型的选择，本研究通过对突发数据进行韦布尔概率图拟合，从而发现突发失效时间数据服从韦布尔分布．5.2　基于竞争失效模型客观贝叶斯估计根据前面的分析，应用竞争失效模型对该组数据进行建模，通过客观贝叶斯方法给出参数的估计，无信息先验分布分别选取Jeffreys先验及Reference先验分布，所得参数估计结果列于表4中．10.13245/j.hust.239299.T004表4竞争失效模型的客观贝叶斯估计先验分布λ/10-4μ/10-3k/103βVAICπJ1.4002.0272.6871.876-24.687πR1.3962.0582.6771.874-22.523当进行客观贝叶斯估计时，由于后验分布复杂，须要运用Gibbs抽样，选取抽样次数为D=1.6×104，将前面D0=1 000个随机数作为老练过程舍弃，贝叶斯估计是基于后面1.5×104个样本获得．通过绘制后验样本的迹线图、遍历平均图和自相关图来分析抽样过程的收敛性．图1~3给出了基于Jeffreys先验参数λ的后验样本迹线图、遍历平均图和自相关图，其余参数及基于Reference先验的各个图类似，由于篇幅限制不再展示，图中It为迭代次数．图3中红色点代表延迟阶数，红色竖线代表自相关系数大小，黑色横线代表相关系数为0，蓝色的上下两条横线分别代表自相关系数的上下界，超出边界的部分表示存在相关关系．可以看出：后验样本之间的自相关不显著，各种先验下Gibbs抽样都收敛了，因此其结果可用．10.13245/j.hust.239299.F001图1基于Jeffreys先验λ的后验样本迹线图10.13245/j.hust.239299.F002图2基于Jeffreys先验λ的后验样本遍历平均图10.13245/j.hust.239299.F003图3基于Jeffreys先验λ的后验样本自相关图由迹线图、遍历平均图和自相关图可知：参数迭代趋于稳定，则认定后验样本的生成均达到收敛．从估计结果可看出：基于不同方法的贝叶斯估计相对接近，从AIC的角度可认为基于Jeffrey先验的客观贝叶斯估计模型更适合拟合此组数据．6 结论本研究提出一种基于TED过程的竞争失效模型的客观贝叶斯分析方法．首先给出了所对应的竞争失效模型，根据其对数似然函数推导了两种无信息先验，并证明了两种后验的合理性及概率匹配性质．然后在此基础上考虑到一些参数后验分布不是典型的分布函数类型，提出了后验估计的抽样算法，即MH算法的Gibbs抽样．在数值算例中，将客观贝叶斯估计与MLE、主观贝叶斯估计进行比较，结果表明所提出的两种客观贝叶斯估计结果都优于MLE和主观贝叶斯估计．最后，利用一个实例进一步证明了所提出方法及所建模型的有效性．本研究的主要工作为：a．提出了适用范围更广的TED过程用于性能退化数据的建模，从而使退化与突发的竞争失效模型更实用，竞争失效模型的适用范围更广；b．针对竞争失效模型，提出了客观贝叶斯估计方法，通过模拟实例分析验证了所提方法的有效性及优势，尤其是当小样本无历史信息时．在工程实际中，先验信息的得到往往并不容易，如果还是小样本情形，更难刻画产品的先验信息，那么会使得传统的MLE及主观贝叶斯估计方法有一定困难．此时，本研究所提的客观贝叶斯估计可以作为工程实际的一个参考．关于更进一步的研究，由于在实际的产品退化试验中通常须要对试验样本施加应力进行加速处理，因此对加速试验下的竞争失效模型的分析可以进行深入研究．另外，由于在高可靠长寿命的产品中更易发生多阶段退化过程，因此研究多阶段退化过程与突发失效的竞争失效模型的参数估计问题也可作为未来的研究方向之一．