杜雪萌，張峰源

（貴州民族大學 數據科學與信息工程學院，貴陽 550025）

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

1 基本假設和模型描述

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

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

基本假設如下：

（1）在應力水平Sj（j=0，1）下，試驗產品失效機理相同。

（2）在應力水平Sj（j=0，1）下，產品的失效時間xji（j=1，2，i=1，2，…，mj）相互獨立。

（3）在應力水平S0下，產品的失效時間x0服從BurrⅢ分布，其分布函數、密度函數、可靠性函數和失效率函數分別為

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

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

2 極大似然估計

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

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

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

式中：M

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

對L（x）取對數后分別對θ，β求偏導，并分別令其偏導數為0可得

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

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

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

3 貝葉斯估計

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

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

因此，參數θ和β的聯合先驗可表示為

結合式（12）和式（18）可得θ和β的聯合后驗密度分布函數為

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

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

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

（2）令i=1；

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

（6）令i=i+1；

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

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

3.1 平方損失函數下的估計

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

3.2 廣義熵損失函數下的估計

3.3 對稱熵損失函數下的估計

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

4 隨機模擬

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

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

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

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

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

5 結論

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