李 玥, 冶繼民

(西安電子科技大學 數(shù)學與統(tǒng)計學院, 西安 710126)

在某些實際問題中, 為了解系統(tǒng)的失效原因和隨時間變化系統(tǒng)的退化情況, 需要在試驗的不同階段將部分未失效系統(tǒng)移離試驗, 即逐步截尾試驗. 逐步截尾方案通常分為逐步Ⅰ型截尾方案(progressively Type-Ⅰ censored scheme, Type-Ⅰ PCS)和逐步Ⅱ型截尾方案(progressively Type-Ⅱ censored scheme, Type-Ⅱ PCS). 逐步Ⅰ型截尾是指在預先給定的間隔觀察時刻點, 將預先給定數(shù)量的未失效系統(tǒng)移離試驗, 直到試驗終止時刻結束試驗. 逐步Ⅱ型截尾是指每當觀測到系統(tǒng)失效時, 將預先給定數(shù)量的未失效系統(tǒng)移離試驗, 直到有預先給定數(shù)量的系統(tǒng)失效時結束試驗. 逐步截尾方案比傳統(tǒng)的Ⅰ型和Ⅱ型截尾方案更靈活, 且Type-Ⅱ PCS與Type-Ⅰ PCS相比, 可以保證至少觀測到給定數(shù)量的失效系統(tǒng), 從而提高試驗的準確性[1-2].

系統(tǒng)壽命數(shù)據(jù)包括失效時間和失效原因兩方面. 在工程實踐中, 引起系統(tǒng)失效的原因多種多樣, 但由于測試環(huán)境、 測試時間等因素的限制, 有時可能無法識別導致系統(tǒng)失效的確切部件, 而只能將其歸咎于一些可能導致系統(tǒng)失效的部件集合. 此時, 系統(tǒng)的失效原因被屏蔽, 這種數(shù)據(jù)稱為屏蔽數(shù)據(jù), 相應的系統(tǒng)稱為屏蔽數(shù)據(jù)系統(tǒng). Usher等[3]采用極大似然方法研究了串聯(lián)屏蔽數(shù)據(jù)系統(tǒng)的可靠性; Sarhan[4-5]研究了當部件壽命服從Weibull分布和Pareto分布時, 屏蔽數(shù)據(jù)系統(tǒng)模型的統(tǒng)計推斷; 蔡靜等[6]研究了Burr Ⅻ分布串聯(lián)屏蔽數(shù)據(jù)系統(tǒng)的可靠性; Sarhan等[7]討論了并聯(lián)屏蔽數(shù)據(jù)系統(tǒng)的部件可靠性評估問題.

上述研究均假設系統(tǒng)部件壽命變量相互獨立, 但在工程實踐中, 系統(tǒng)部件之間通常不獨立, 部件失效存在一定的相關關系. 例如, 一架雙引擎飛機, 當其中一個引擎出現(xiàn)故障時, 另一個引擎的工作狀態(tài)會發(fā)生變化, 從而維持飛機能正常飛行. 因此, 研究部件之間的相依關系及其對系統(tǒng)可靠性的影響有一定意義. Feizjavadian等[8]和Bai等[9]利用二元Marshall-Olkin Weibull分布描述失效原因之間的相依結構, 對相依串聯(lián)系統(tǒng)模型進行了統(tǒng)計分析； Cai等[10]基于二元Marshall-Olkin Weibull分布對相依串聯(lián)屏蔽數(shù)據(jù)系統(tǒng)進行了統(tǒng)計分析； 付倩嬈等[11]討論了二元Marshall-Olkin指數(shù)分布下含有屏蔽數(shù)據(jù)的相依串聯(lián)系統(tǒng)的可靠性分析問題. 上述研究均采用多元分布刻畫相依結構, 適用于系統(tǒng)部件壽命變量具有相同分布的結構. 但在復雜結構的情形下, 多元分布對變量之間的不對稱性和尾部相依性的計算較困難. Copula函數(shù)因其靈活性可以很好地解決上述問題, 本文選擇Copula函數(shù)進行相依結構的建模. Hsu等[12]基于Clayton Copula函數(shù)研究了在Ⅰ型截尾試驗下多部件相依串聯(lián)系統(tǒng)的可靠性評估問題； 蔡靜等[13]基于Clayton Copula函數(shù)討論了相依并-串聯(lián)屏蔽數(shù)據(jù)系統(tǒng)的可靠性分析問題.

Burr Ⅻ分布[14]是一種重要的連續(xù)壽命分布. 由于其失效率函數(shù)的非單調形式, 使其適合代表許多產品的壽命, 具有很好的靈活性, 廣泛應用于可靠性工程、 通信工程、 質量控制、 航空航天等領域. Burr Ⅻ分布作為一種失效模型, 關于其應用和統(tǒng)計性質目前已有許多研究結果. 但基于Burr Ⅻ分布對相依屏蔽數(shù)據(jù)系統(tǒng)的研究報道較少, 因此, 本文針對部件壽命變量服從Burr Ⅻ分布的情形, 通過引入Copula函數(shù)建立部件壽命變量之間的相依結構, 并討論在逐步Ⅱ型截尾試驗下相依串聯(lián)屏蔽數(shù)據(jù)系統(tǒng)的可靠性. 首先, 推導出模型參數(shù)及系統(tǒng)可靠度函數(shù)的極大似然估計(maximum likelihood estimations, MLEs)； 然后, 基于漸近正態(tài)性理論和bootstrap抽樣算法分別構造出參數(shù)的漸近置信區(qū)間(asymptotic confidence intervals, ACIs)及偏差校正的百分位bootstrap置信區(qū)間(bias-corrected percentile bootstrap confidence intervals, Boot-BCP CIs)； 最后, 用仿真模擬和真實數(shù)據(jù)分析驗證所提出統(tǒng)計方法的可行性和有效性.

1 Copula理論

1.1 Copula函數(shù)

Copula函數(shù)是描述多元隨機變量間相依結構的一種函數(shù), 它可以將多元隨機變量的聯(lián)合分布函數(shù)與其邊際分布函數(shù)聯(lián)系起來. 下面簡要介紹二元變量情形下Copula函數(shù)的基本理論[15-16].

失效Copula函數(shù)是一個多維聯(lián)合分布函數(shù), 其定義域為[0,1]×…×[0,1]. 根據(jù)Sklar定理[15], 它可以連接多個隨機變量的邊際分布構造聯(lián)合分布. 對于二元變量的情形, 其可表示為

H(z1,z2)=C(F1(z1),F2(z2)),

(1)

二元生存Copula函數(shù)和二元失效Copula函數(shù)的關系為

1.2 相依性度量

通常的相依性度量如Spearman’sρ和Kendall’sτ都可采用Copula函數(shù)表示[15]為

1.3 Archimedean Copula族

由于Archimedean Copula族可以準確地描述非對稱隨機變量之間的尾相關性, 因此在實際生活中應用廣泛. 二元Archimedean Copula族函數(shù)定義[15]為C(μ,ν)=φ-1(φ(μ)+φ(ν)), 其中生成函數(shù)φ(·)是一個嚴格單調遞減的凸函數(shù), 滿足φ(μ)+φ(ν)≤φ(0),φ(1)=0,μ,ν∈[0,1],φ-1(·)是φ(·)的逆.

Clayton Copula隸屬于Archimedean Copula族, 其Copula函數(shù)為

CC(μ,ν)=[μ-θ+ν-θ-1]-1/θ,θ∈[-1,∞){0}.

Clayton Copula適用于描述下尾相關、 上尾漸近獨立隨機變量之間的相關性, 由Kendall’sτ及相應Copula函數(shù)的關系可得τ=θ/(θ+2), 當θ→0時, 隨機變量趨于相互獨立.

2 模型描述和基本假設

2.1 模型描述

由于某些系統(tǒng)的失效原因并不確切, 因此用Ki(i=1,2,…,m)表示系統(tǒng)i可能失效原因的集合,ki是Ki的觀測值.ki={1}表示系統(tǒng)i失效是由部件1導致的；ki={2}表示系統(tǒng)i失效是由部件2導致的；ki={1,2}表示系統(tǒng)i的失效原因被屏蔽.因此, 觀測到的系統(tǒng)壽命數(shù)據(jù)為(t1,k1,r1),(t2,k2,r2),…,(tm,km,rm).

2.2 基本假設

1) 設Z1,Z2分別表示部件1和2的失效時間, 各部件的失效時間Zv(v=1,2)服從參數(shù)為γv,α的Burr Ⅻ分布, 簡記為Zv～BX(γv,α), 其中γv>0,α>0.系統(tǒng)的失效時間為各部件失效時間的最小值, 即T=min{Z1,Z2}.

部件v失效時間的累積分布函數(shù)(cumulative distribution function, CDF)、 概率密度函數(shù)(probability density function, PDF)、 生存函數(shù)(survival function, SF)(也稱可靠度函數(shù))分別為

2) 假設系統(tǒng)由兩個相依的部件串聯(lián)而成, 且其失效是由最先失效的部件所致.為考察數(shù)據(jù)間的下尾相依性, 本文選擇Clayton Copula刻畫各部件之間的相依關系, 則系統(tǒng)的SF為

S(t)=SZ1,Z2(z1,z2)|(t,t)=P(Z1>t,Z2>t)=[(1+tα)γ1θ+(1+tα)γ2θ-1]-1/θ.

3) 等概率假設： 屏蔽的發(fā)生與失效原因無關, 即

P(Ki=ki|Ti=ti,Vi=v′)=P(Ki=ki|Ti=ti,Vi=v), ?v′,v∈{1,2},

其中,Ti和Vi分別表示系統(tǒng)i的失效時間和確切失效原因,p?P(Ki=ki|Ti=ti,Vi=v)表示屏蔽概率.為計算簡便, 本文基于等概率假設進行討論, 也有部分研究放寬了該假設[12,17-18].

3 極大似然估計

由部件1導致系統(tǒng)i失效的PDF為

由部件2導致系統(tǒng)i失效的PDF為

則失效系統(tǒng)的似然函數(shù)可表示為

(4)

為簡化似然函數(shù), 定義失效原因示性函數(shù)為

定義輔助示性函數(shù)為

因此基于假設3), 似然函數(shù)可表示為

對數(shù)似然函數(shù)為

對lnL(γ1,γ2,α,θ)關于γ1,γ2,α,θ求一階導數(shù)并令其等于0, 可得似然方程為

由于很難得到模型參數(shù)的顯式表達式, 因此本文考慮迭代算法, 如Newton-Raphson迭代算法或其他迭代算法求解γ1,γ2,α,θ的MLEs.則系統(tǒng)可靠度函數(shù)的估計可表示為

(11)

4 區(qū)間估計

4.1 漸近置信區(qū)間

記參數(shù)γ1,γ2,α,θ的Fisher信息矩陣為

其中,L(·,·)表示lnL(γ1,γ2,α,θ)對參數(shù)γ1,γ2,α,θ的二階偏導數(shù), 其中,

則參數(shù)γ1,γ2,α,θ的100(1-ε)%雙側ACIs為

其中Zε/2是標準正態(tài)分布的ε/2分位數(shù).

4.2 Bootstrap置信區(qū)間

在現(xiàn)代制造業(yè)中, 由于時間和成本的限制, 樣本量通常不大, 基于大樣本數(shù)據(jù)的統(tǒng)計方法可能并不合適, 有時甚至會產生誤導. 因此, 本文采用bootstrap抽樣算法構造參數(shù)的置信區(qū)間.

5) 參數(shù)γ1,γ2,α,θ的100(1-ε)%雙側Boot-BCP CIs為

φ-1(·)是φ(·)的逆函數(shù),I(·)表示示性函數(shù).

5 仿真模擬

下面利用Monte Carlo仿真模擬評估上述統(tǒng)計方法的性能并分析模擬結果. 考慮一個具有兩個相依部件的串聯(lián)系統(tǒng), 假設部件的壽命分布服從Burr Ⅻ分布, 利用二維Clayton Copula函數(shù)刻畫其相依性. 選取參數(shù)γ1=1,γ2=1.2,α=1, 給定樣本量n和失效系統(tǒng)數(shù)m以及截尾方案{r1=r2=…=rm-1=0,rm=n-m}, 在不同的屏蔽概率(p=0.3,0.5)和不同的相依關系(θ=2,4)下, 重復試驗5 000次, 并計算各參數(shù)的MLEs、 均方誤差(mean square errors, MSEs), 以及ACIs和Boot-BCP CIs的平均寬度(average widths, AWs)、 覆蓋率(coverage probabilities, CPs), 結果分別列于表1～表3. 當時間t=1,2,3時系統(tǒng)的可靠度分別為0.347 6,0.213 4,0.152 8.

表1 當p=0.3, θ=2時, 不同樣本量下模型參數(shù)的點估計和區(qū)間估計

表2 當p=0.5, θ=2時, 不同樣本量下模型參數(shù)的點估計和區(qū)間估計

表3 當p=0.3, θ=4時, 不同樣本量下模型參數(shù)的點估計和區(qū)間估計

由表1～表3可見: 隨著樣本量的增大, 點估計中參數(shù)的MLEs逐漸接近真值, MSEs逐漸減小； 區(qū)間估計中, ACI-AWs和BCI-AWs逐漸減小, ACI-CPs和BCI-CPs逐漸增大, 這是因為隨著n或m的增加, 試驗得到的有效樣本量增多, 從而估計結果更準確. 隨著屏蔽概率的增大, MSEs,AWs和CPs的結果表明, 參數(shù)點估計和區(qū)間估計的準確性都會降低, 這是因為屏蔽數(shù)據(jù)越多, 試驗得到的有效信息越少, 因此在分析壽命數(shù)據(jù)時, 屏蔽數(shù)據(jù)的存在不可忽略. 隨著部件壽命變量之間相依性的增強, 參數(shù)的MLEs接近真值, ACI-AWs,BCI-AWs都減小, ACI-CPs,BCI-CPs都增大, 表明在一定程度上考慮部件壽命變量之間的相依性是必要的.

為進一步說明屏蔽概率對系統(tǒng)可靠度估計的影響, 表4列出了當p=0.3和p=0.5時系統(tǒng)可靠度估計值以及其與真實值的誤差.由表4可見, 隨著樣本量的增加, 系統(tǒng)可靠度估計值越來越接近真值.屏蔽概率的增大會導致降低系統(tǒng)可靠度估計的精度.

表4 當p=0.3和p=0.5時系統(tǒng)可靠度的比較

當p=0.5,θ=2時, 刪除失效原因被屏蔽的數(shù)據(jù), 基于失效部件明確的數(shù)據(jù)構造參數(shù)的極大似然估計, 所得結果列于表5. 對比表2和表5可見, 參數(shù)極大似然估計的MSE增大. 表明忽略屏蔽數(shù)據(jù)將導致增大參數(shù)估計的誤差, 因此, 屏蔽數(shù)據(jù)的存在不可忽略.

表5 當p=0.5, θ=2時, 刪除屏蔽數(shù)據(jù)后模型參數(shù)的估計

6 實例分析

考慮用文獻[21]中對36個小型電器進行壽命試驗的實際壽命數(shù)據(jù)集驗證本文方法的可行性. 在該數(shù)據(jù)集中, 失效原因被分為18種不同的模式, 觀察到的大多數(shù)故障都是由模式9引起的, 因此本文將失效模式9視為失效原因1, 將其他所有失效模式視為失效原因2. 給定n=36,m=30以及截尾方案{r1=6,r2=r3=…=r30=0}, 得到含有屏蔽數(shù)據(jù)的逐步Ⅱ型截尾數(shù)據(jù)列于表6.

表6 含有屏蔽數(shù)據(jù)的逐步Ⅱ型截尾數(shù)據(jù)

基于含有屏蔽數(shù)據(jù)的逐步Ⅱ型截尾數(shù)據(jù), 模型參數(shù)γ1,γ2,α,θ相應的MLEs,ACIs和Boot-BCP CIs列于表7. 由表7可見, ACIs和Boot-BCP CIs均包含相應的點估計, Boot-BCP CIs在區(qū)間長度上性能優(yōu)于ACIs.

表7 模型參數(shù)估計

為檢驗本文模型的擬合優(yōu)度, 采用Kolmogorov-Smirnov(K-S)檢驗的K-S距離和p值(表8)以及如圖1所示的Quantile-Quantile(Q-Q)圖說明Burr Ⅻ分布的擬合情況. 由表8可見, 當選取置信度為0.05時, 原因1和原因2的K-S距離均小于相應的臨界值, 且相應的p值均大于0.05. Q-Q圖也表明, 擬合的Burr Ⅻ分布是一個合適的分布. 因此, 本文提出的模型可行且有效.

表8 K-S檢驗

圖1 Q-Q圖Fig.1 Q-Q plots

表9 Copula模型與獨立模型的比較

綜上所述, 本文針對部件壽命變量服從Burr Ⅻ分布的情形, 基于Clayton Copula函數(shù), 討論了含有屏蔽數(shù)據(jù)的相依串聯(lián)系統(tǒng)的可靠性, 得到了模型參數(shù)及系統(tǒng)可靠度函數(shù)的MLEs, 并基于漸近正態(tài)性理論和bootstrap抽樣算法構造了模型參數(shù)的ACIs和Boot-BCP CIs. 仿真模擬和真實數(shù)據(jù)分析驗證了所提出統(tǒng)計方法的可行性與有效性. 在小樣本情形下, 建議采用偏差校正的百分位bootstrap方法進行區(qū)間估計. 在系統(tǒng)壽命試驗的分析中屏蔽數(shù)據(jù)的存在不可忽略, 隨著屏蔽概率的增大模型參數(shù)估計的精度逐漸降低.