楊 超

(溫州大學數(shù)理與電子信息工程學院，浙江溫州 325035)

在加速壽命試驗中，通過可靠性分析考慮被測單元、產(chǎn)品和系統(tǒng)的各種失效原因的方法稱為競爭風險模型[1].Wu等[2]通過使用I型漸進混合截尾方法，在競爭風險的壽命服從威布爾分布時，通過恒定應力加速壽命試驗得到最大似然估計和貝葉斯估計[2].Mao等[3]給出了具有指數(shù)競爭失效模型的廣義I型混合截尾數(shù)據(jù)的精確似然推理.基于未知參數(shù)的極大似然估計，利用條件矩母函數(shù)建立參數(shù)的精確條件分布，討論了參數(shù)的矩性質(zhì)和精確約束區(qū)間[3].Balakrishnan等[4]將競爭風險模型引入加速壽命試驗環(huán)境下的單次器件測試分析中，并開發(fā)了EM算法來計算模型參數(shù)的估計量.Zhang等[5]分析了漸進式I型混合審查下具有競爭風險的部分加速壽命試驗中的競爭風險模型，基于被篡改的失效率模型，得到了Weibull未知參數(shù)和加速度因子的最大似然估計、貝葉斯估計、漸近CI和最高后驗密度可信區(qū)間.

但是，競爭風險的獨立性假設(shè)在可靠性分析中具有較大的局限性.作為基本模型，二元競爭風險模型被廣泛用于討論多元競爭風險.在m（＞2）個元件的壽命測試中，m個元件競爭失效的原因可分為兩個部分.一個是導致故障時間的失效原因，另一個原因是其他m-1個元件競爭失效.因此，考慮多元競爭風險之間的相關(guān)性對于統(tǒng)計推斷是非常重要和必要的.

1 模型設(shè)定

1.1 多元Birnbaum-Saunders分布

若T的累積分布函數(shù)可以顯示如（1）式，那么隨機變量T就具有兩個參數(shù)的Birnbaum-Saunders分布（Birnbaum-Saunders分布以下簡稱BS分布，多元Birnbaum-Saunders分布以下簡稱MBS分布）：

其中Φ(·)是標準正態(tài)累積分布函數(shù)，α＞0和β＞0是形狀參數(shù)和尺度參數(shù).也就是說，T～BS(α,β).

T的相應概率密度函數(shù)可以由（2）式給出：

其中參數(shù)α1＞0,α2＞0,…,αp＞0，β1＞0,β2＞0,…,βp＞0，Σ是p×p正定擴散矩陣對于t1＞0,t2＞0,…,tp＞0.Φp(u,v;Σ)是一個標準的多元正常隨機向量(u,v)的累積分布函數(shù)與正定擴散矩陣Σ.

T1,T2,…,Tp的聯(lián)合概率密度函數(shù)由（4）式給出：

聯(lián)合概率分布函數(shù)可以表示為

其中u=(u1,u2,…,up)T.

下面的定理中可以提供多元MBS分布的邊際分布和條件分布.

定理1 若T～MBSp(α,β,Σ)，讓T,α,β,Σ按如下方式劃分

其中T1,α1,β1全是q×1維向量，Σ11是一個q×q維的矩陣，其余的元素與之相應匹配，可得出：

1）T1～MBSq(α1,β1,Σ11)以及T2～MBSp-q(α2,β2,Σ22)；

2）當T2=t2，T1的條件累積分布函數(shù)是

3）當T2=t2，T1的條件概率密度函數(shù)是

證明：

對于1）的證明，可通過令（3）式中tq+1→∞,…,tp→∞簡單得出.

1.2 相依競爭風險模型

考慮p個相依競爭風險因素T1,T2,…,Tp和邊際BS分布，也就是(T1,T2,…,Tp)～MBS(α,β;Σ)，且Tj～BS(αj,βj)，j=1,2.T1,T2,…,Tp之間的相依結(jié)構(gòu)是由多變量隨機向量的MBS分布(T1,T2,…,Tp)描述的，聯(lián)合生存函數(shù)和累積分布函數(shù)(T1,T2,…,Tp)可以分別表示為：

令T=min(T1,T2,…,Tp)作為一個系統(tǒng)的壽命，對于j=1,2,…,p的次密度函數(shù)fj(t)可以被計算為：

T的生存函數(shù)是S(t)=ST1,T2,…,Tp(t,t,…,t)，在相依競爭的風險模型中，潛在的失效時間T=min(T1,T2,…,Tp)和在一組{1,2,…,p}中有一個整數(shù)值的失效索引C，它們形成相依競爭的風險數(shù)據(jù)(T,C).

1.3 基于自適應逐步混合檢測方法的加速壽命試驗的數(shù)據(jù)描述

在k級恒定應力加速壽命試驗中，有k個應力水平s1＜s2＜…＜sk和正常應力水平s0.在每個應力水平si，p依賴失效因子(Ti1,Ti2,…,Tip)，服從MBS(αi1,βi1,αi2,βi2,…,αip,βip;Σ)分布，其中i=1,2,…,k.將n個測試單元置于加速壽命試驗中，并用ni個測試單元將它們分組成k組.在si處對ni個測試單元進行測試，并指定預定時間τi、預期ri故障并在故障時間Ril(l=1,2,…,ri)移除單元.

給定恒應力加速壽命試驗中的預定時間τi，可以在具有相應終止時間min{Ti,ri:ri:ni,τi}和max{Ti,ri:ri:ni,τi}的type-I型漸進混合檢測方法和type-II型PHCS中觀察到應力水平si處的故障時間，其中Ti,ri:ri:ni是應力水平si下的第ri個故障時間.在type-I型PHCS中，當測試單元具有長壽命時，可能發(fā)生零失效時間.在type-II型PHC中，實驗可以在足夠長的時間之后結(jié)束.這些缺點會影響統(tǒng)計推斷的效率.Gasbarra等[6]提出了自適應漸進混合檢測方法，通過重新確定Ril值，自適應PHCS提供了克服PHCS缺陷的有效方法.

如果第ri個逐漸刪失的失效時間發(fā)生在時間τi之前，則si下的實驗在故障時間Ti,ri:ri:ni終止.也就是說則加速測試繼續(xù)，為了盡快到達第ri個逐漸審查的故障時間，在時間τi之后，設(shè)定

如Gasbarra等[6]中提出的這種設(shè)定使我們在時間iτ之后盡快結(jié)束實驗，此時直到在應力水平si處觀察到ri次失效，然后給出加速競爭風險數(shù)據(jù).對于如果失敗指數(shù)cil=1，其中δil=1；否則δil=0.

1）應力水平僅改變BS分布的尺度參數(shù)，因此α1j=α2j=…=αkj=αj,j=1,2,…,p.

2）尺度參數(shù)與應力水平之間的關(guān)系一般是對數(shù)線性的，即

基于上述假設(shè)和數(shù)據(jù)描述，得到未知參數(shù)θ=(α1,α2,…,αp,βi1,βi2,…,βip;Σ)，i=1,2,…,k，其中αj＞0，Σ＞0，β1j＞β2j＞…＞βkj＞0當j=1,2…,p.

2 極大似然估計

在上述加速競爭風險模型的假設(shè)下，似然函數(shù)被表示為

對于j=1,2,…,p在故障時間內(nèi)，在應力水平S1下，f1j(til)是次密度函數(shù)，S1(til)是生存函數(shù).Lil(θ|til)的對數(shù)似然函數(shù)表示為（13）式，因此有關(guān)于單個參數(shù)的似然方程

由（14）給出.

由于（14）式中方程的復雜形式可能在顯式表達式中不能得到未知參數(shù)的顯式極大似然估計量，可通過數(shù)值迭代算法，如Newton-Raphson算法得到極大似然估計的數(shù)值結(jié)果.

3 統(tǒng)計推斷

3.1 初始值設(shè)定

在競爭風險模型中，失敗的原因來自于不同參數(shù)的MBS分布.基于這種相依結(jié)構(gòu)，競爭風險之間的獨立性以及多個競爭風險的形狀和規(guī)模參數(shù)的關(guān)系可能是值得關(guān)注的.討論了下列假設(shè)的似然比檢驗：

在試驗中，H0情況下通過求解似然方程得到參數(shù)的極大似然估計：

3.2 極大似然估計

實際應用中，在壽命試驗中需要對試驗者、制造商和客戶提供單位或產(chǎn)品的界限.具有相關(guān)競爭風險，在故障時間til（l=1,2,…,ri，i=1,2,…,k）有幸存的單元被移除，因此測試單元的故障時間是未觀測到的.在該條件下，用某個測試單位表示某個元件的第s個統(tǒng)計數(shù)據(jù).

已知極大似然估計量

si下的的條件密度函數(shù)可表示為：

其中

4 結(jié)論與討論

考慮了多元MBS分布在三個部件系統(tǒng)的情形，設(shè)定初始值：

在常應力下，我們從MBS分布中分析具有相依競爭風險的自適應漸進混合截尾數(shù)據(jù)，在三個部件的情形下通過設(shè)定抽取樣本數(shù)為5 000，得到相應的極大似然估計和矩估計.在興趣假設(shè)的基礎(chǔ)上，討論了似然比檢驗.進行了仿真研究，得到了估計量和預測值的數(shù)值結(jié)果，見表1.

表1 極大似然估計和矩估計Fig 1 Maximum Likelihood Estimation (MLE) and Moment Estimation (ME)

我們比較了極大似然估計和矩估計下各個參數(shù)的偏差與均方誤差，二者模擬結(jié)果比較接近，偏差與均方誤差數(shù)值均比較小，模擬效果較好.仿真研究表明，所提出的模型和方法可用于基于競爭風險數(shù)據(jù)的統(tǒng)計推斷.在以后的工作中，可以進一步研究加速壽命試驗下自適應漸進混合截尾數(shù)據(jù)下基于相關(guān)競爭風險模型的貝葉斯推理.