杜雪萌，張峰源

（貴州民族大學(xué) 數(shù)據(jù)科學(xué)與信息工程學(xué)院，貴陽 550025）

由于壽命試驗產(chǎn)品的壽命時間越來越長，在進行生存分析與可靠性試驗時很難在正常情況下獲得這些產(chǎn)品的故障數(shù)據(jù)。因此，為了縮短測試時間，加速性能退化，對所有或部分測試產(chǎn)品施加比平時更嚴(yán)重的壓力條件，其被稱為加速壽命測試。如果只將部分產(chǎn)品置于加速條件下，則該試驗又稱為部分加速壽命試驗。近年來，有不少文獻對加速壽命試驗的統(tǒng)計分析進行了研究。Almarashi[1]基于恒定應(yīng)力部分加速壽命試驗，得到了逐步Ⅱ型刪失數(shù)據(jù)下廣義倒指數(shù)分布的極大似然估計和Bayes估計；Kumar等[2]在恒定應(yīng)力加速壽命試驗下，利用9種常用估計方法對Lindley分布進行了參數(shù)估計；武東等[3]基于定數(shù)刪失數(shù)據(jù)，對瑞利分布恒定應(yīng)力加速壽命試驗進行了貝葉斯統(tǒng)計分析。考慮到逐步Ⅱ型刪失數(shù)據(jù)的優(yōu)良性，本文將基于逐步Ⅱ型刪失數(shù)據(jù)，在恒定應(yīng)力部分加速壽命試驗下，對BurrⅢ型分布的形狀參數(shù)和加速因子進行參數(shù)估計。

1 基本假設(shè)和模型描述

逐步Ⅱ型刪失下BurrⅢ分布恒定應(yīng)力部分加速壽命試驗的具體試驗設(shè)計方案如下。

假設(shè)有n個壽命服從BurrⅢ分布且相互獨立的試驗樣品，將其隨機分成2組，其容量分別為n0，n1。其中，n0個樣品被安排到正常應(yīng)力水平S0下進行壽命試驗，n1個樣品被安排到加速應(yīng)力水平S1（S0＜S1）下進行壽命試驗。在應(yīng)力水平Sj（j=0，1）下進行逐步Ⅱ型刪失試驗。當(dāng)觀測到第一個試驗樣品失效時，從剩余的nj-1個未失效的樣品中隨機挑選出Rj1個樣品退出試驗，當(dāng)觀測到第二個試驗樣品失效時，再從剩余的nj-Rj1-2個未失效的樣品中隨機挑選出Rj2個樣品退出試驗，按照這種方法一直試驗，當(dāng)觀測到第mj（mj＜nj）個失效樣品時，將剩下的個樣品全部撤出試驗。共m（m=m0+m1）個失效樣本。記在應(yīng)力水平Sj（j=0，1）觀測到的失效時刻為xj=（xj1，xj2，…，xjmj）。本文假設(shè)從試驗中隨機剔除的Rji完全隨機，保證Rj1+Rj2+…+Rjm=nj-mj即可。

基本假設(shè)如下：

（1）在應(yīng)力水平Sj（j=0，1）下，試驗產(chǎn)品失效機理相同。

（2）在應(yīng)力水平Sj（j=0，1）下，產(chǎn)品的失效時間xji（j=1，2，i=1，2，…，mj）相互獨立。

（3）在應(yīng)力水平S0下，產(chǎn)品的失效時間x0服從BurrⅢ分布，其分布函數(shù)、密度函數(shù)、可靠性函數(shù)和失效率函數(shù)分別為

式中：α＞0，θ＞0分別稱為尺度參數(shù)和形狀參數(shù)，本篇文章假設(shè)尺度參數(shù)α已知。

（4）在應(yīng)力水平S1下，產(chǎn)品的失效率函數(shù)為h1（x）=βh0（x），β＞1，式中：β為加速因子。則應(yīng)力水平S1下產(chǎn)品失效時間x1的分布函數(shù)、密度函數(shù)、可靠性函數(shù)和失效率函數(shù)分別為

2 極大似然估計

基于上述模型描述，通過極大似然方法[4]對BurrⅢ型分布的形狀參數(shù)θ和加速因子β進行估計，可以得到如下定理。

定理1：按照上述逐步Ⅱ型刪失恒定應(yīng)力部分加速壽命試驗，BurrⅢ型分布形狀參數(shù)θ的極大似然估計由方程（9）確定，并且該方程有唯一解，同時，加速因子β極大似然估計式為式（10）。

證明：假設(shè)xj=（xj1，xj2，…，xjmj）是應(yīng)力水平Sj（j=0，1）下服從BurrⅢ型分布的逐步Ⅱ型刪失數(shù)據(jù)，記x=（x0，x1）為恒定應(yīng)力壽命試驗下服從BurrⅢ型分布的逐步Ⅱ型刪失數(shù)據(jù)。Rj=（Rj1，Rj2，…，Rjmj）為應(yīng)力水平Sj（j=0，1）下對應(yīng)的刪失模式，參數(shù)α已知，則關(guān)于x的似然函數(shù)為

式中：M

將式（1）、式（2）、式（7）和式（8）帶入式（11）可得

對L（x）取對數(shù)后分別對θ，β求偏導(dǎo)，并分別令其偏導(dǎo)數(shù)為0可得

參數(shù)θ顯然無法由式（13）得到顯式解，為此，探求此方程在（0，+∞）上是否具有唯一解。

顯然g（θ）′＜0，同時因此方程（13）在（0，+∞）上有唯一解。

式（13）沒有顯式表達式，但有唯一解，故考慮應(yīng)用Newton-Raphson近似法求出近似解。記此解為參數(shù)θ的極大似然估計近似值θ^M，將θ^M帶入式（14）中可得到參數(shù)β的極大似然估計近似值β^M。

3 貝葉斯估計

在一些實際情況下，參數(shù)值的信息可以以獨立的方式獲得。因此，假設(shè)參數(shù)先驗獨立，設(shè)加速度因子β的先驗分布為無信息先驗，即

取θ的先驗分布為伽馬分布，其形狀參數(shù)與尺度參數(shù)分別為a和b，概率密度函數(shù)為

因此，參數(shù)θ和β的聯(lián)合先驗可表示為

結(jié)合式（12）和式（18）可得θ和β的聯(lián)合后驗密度分布函數(shù)為

參數(shù)θ和β的全條件后驗分布函數(shù)分別表示為

基于上述全條件后驗密度分布函數(shù)式（20）和式（21），本文考慮使用Metropolis-Hastings（MH）方法產(chǎn)生Markov Chain Monte Carlo（MCMC）樣本并估計參數(shù)。選取提議分布為正態(tài)分布和具體算法如下：

（1）給出初始值θ(0)，β(0)；

（2）令i=1；

（5）計算θ（i）和β（i）；

（6）令i=i+1；

（7）重復(fù)步驟（2）—（6）M次，生成樣本θ(1)，θ(2)，…，θ(M)與β(1)，β(2)，…，β(M)。

下面利用基于上述方法所產(chǎn)生的樣本θ(1)，θ(2)，…，θ(M)與β(1)，β(2)，…，β(M)對參數(shù)θ和加速因子β進行估計。為方便書寫，記γ＝（θ，β）。

3.1 平方損失函數(shù)下的估計

設(shè)δ是參數(shù)γ的任一決策函數(shù)，則平方損失下L1=（θ-δ）2，γ的貝葉斯估計為

3.2 廣義熵損失函數(shù)下的估計

3.3 對稱熵損失函數(shù)下的估計

不同損失函數(shù)下，γ估計的后驗均方誤差為

4 隨機模擬

利用Monte-Carlo方法通過R軟件產(chǎn)生一個服從BurrⅢ分布的恒定應(yīng)力部分加速壽命試驗逐步Ⅱ型刪失樣本，具體步驟如下：①產(chǎn)生2個容量分別為n0，n1且服從均勻分布U0（0，1），U0（0，1）的獨立同分布樣本，并升序排列為U01，U02，…，U0n0和U11，U12，…，U1n1；②當(dāng)α=2時，設(shè)θ=1.5，β=1.2

則X01，…，X0n0，X11，…，X1n1是一個恒定應(yīng)力部分加速壽命試驗下服從參數(shù)θ=1.5，加速因子β=1.2的BurrⅢ分布的獨立同分布樣本；③根據(jù)不同樣本量n0，n1和觀測次數(shù)m0，m1，隨機生成相應(yīng)的刪失模式R0=（R01，R02，…，R0m0），R1=（R11，R12，…，R1m1）。按所生成的R1，R2進行隨機抽樣，獲得的失效數(shù)據(jù)X0（1），…，X0（m0），X1（1），…，X1（m1），為服從BurrⅢ分布的恒定應(yīng)力部分加速壽命試驗逐步Ⅱ型刪失樣本。

基于上述所生成的刪失樣本，取參數(shù)初值θ=1.8，β=1.4，設(shè)定精度e=0.000 1，則由式（9）和式（10）可求得參數(shù)θ和加速因子β極大似然估計的近似值重復(fù)以上過程1 000次，可求得的均值與均方誤差。根據(jù)上述MH算法步驟，令M=20 000，生成20 000個θ和β，利用所生成的樣本θ(1)，…，θ(20000)，β(1)，…，β(20000)代入式（22）—（24），可得到在3種損失函數(shù)下的θ和β的估計值，其中，取廣義熵損失函數(shù)中參數(shù)c為0.8和1。利用式（25）可求得上述3種損失函數(shù)下估計的后驗均方誤差（MSE）。

表1和表2的模擬結(jié)果表明，在樣本量相同的情況下，觀測次數(shù)不同時，觀測次數(shù)越多時，參數(shù)估計的均方誤差越小，估計效果越好；在觀測次數(shù)相同，樣本量n不同時，n越大，參數(shù)估計的MSE越小。2種估計方法得到的結(jié)果都比較穩(wěn)健。

表1 參數(shù)θ=1.5，β=1.2，不同樣本量n，不同觀測次數(shù)m時的θ的模擬結(jié)果

5 結(jié)論

本文基于恒定應(yīng)力部分壽命試驗，在逐步Ⅱ型刪失數(shù)據(jù)下，首先運用極大似然法對BurrⅢ分布的形狀參數(shù)和加速因子進行估計得到了形狀參數(shù)的估計方程式，并證明該方程式的解是唯一存在的，由于形狀參數(shù)的估計方程式無顯性表達式，采用Newton-Raphson法求出近似解。其次基于MH算法得到了不同損失函數(shù)下未知參數(shù)的貝葉斯估計式。最后運用Monte-Carlo模擬方法，對在不同情況下的形狀參數(shù)與加速因子的估計值進行了模擬比較。結(jié)果表明，定義的極大似然估計和貝葉斯估計的MSE均隨樣本量的增加而減小，2種估計效果近似。