左刪失恒定應(yīng)力部分加速壽命試驗(yàn)下逆Rayleigh分布的參數(shù)估計(jì)

龍兵，張忠占

（1.荊楚理工學(xué)院數(shù)理學(xué)院，湖北荊門448000；2.北京工業(yè)大學(xué)應(yīng)用數(shù)理學(xué)院，北京100124）

0 引 言

在可靠性壽命試驗(yàn)中經(jīng)常會(huì)遇到一些長(zhǎng)壽命的產(chǎn)品，在正常使用條件下進(jìn)行試驗(yàn)，通常需要花費(fèi)很長(zhǎng)的時(shí)間。為了更快地獲得失效數(shù)據(jù)進(jìn)行可靠性分析，可以選擇采用加速壽命試驗(yàn)，即將產(chǎn)品置于更嚴(yán)酷的應(yīng)力條件下試驗(yàn)，這樣可以加速性能的退化，使產(chǎn)品更快地失效。將所有的受試產(chǎn)品都放在加速環(huán)境下試驗(yàn)，稱為加速壽命試驗(yàn)，若只放一部分產(chǎn)品，則稱部分加速壽命試驗(yàn)。由恒定應(yīng)力部分的加速壽命試驗(yàn)數(shù)據(jù)可估算得到加速因子，從而實(shí)現(xiàn)兩種應(yīng)力水平下壽命數(shù)據(jù)的轉(zhuǎn)換，節(jié)省試驗(yàn)時(shí)間和經(jīng)費(fèi)，并能很快得到統(tǒng)計(jì)分析結(jié)果。部分加速壽命試驗(yàn)已經(jīng)引起了一些統(tǒng)計(jì)學(xué)者的重視，形成了很多研究成果。ISMAIL[1-2]基于截尾樣本，在恒定應(yīng)力部分加速壽命試驗(yàn)?zāi)Ｐ拖拢懻摿薟eibull分布的參數(shù)及加速因子的極大似然估計(jì)，確定了最優(yōu)試驗(yàn)方案。ISMAIL等[3-4]在恒定應(yīng)力部分加速壽命試驗(yàn)?zāi)Ｐ拖拢蟮昧薎I型Pareto分布和逆Weibull分布的參數(shù)及加速因子的極大似然估計(jì)，通過(guò)最小化廣義漸近方差給出了兩種應(yīng)力水平下樣本的最優(yōu)分配方案。文獻(xiàn)[5-6]也在恒定應(yīng)力加速壽命試驗(yàn)?zāi)Ｐ拖拢懻摿薖areto分布和逆Weibull分布的可靠性問(wèn)題。

逆Rayleigh分布作為可靠性領(lǐng)域的一個(gè)重要分支，其統(tǒng)計(jì)性質(zhì)值得深入研究。針對(duì)該分布的研究成果已有很多，JEBELY等[7]在完全觀測(cè)樣本下討論了逆Rayleigh分布的密度函數(shù)和分布函數(shù)的估計(jì)。文獻(xiàn)[8-10]在完全觀測(cè)樣本下用貝葉斯方法研究了逆Rayleigh分布的統(tǒng)計(jì)推斷問(wèn)題。在壽命測(cè)試試驗(yàn)中，截尾是大多數(shù)壽命數(shù)據(jù)集不可避免的特征，NAVID等[11]討論了在單截尾和雙截尾數(shù)據(jù)下逆Rayleigh分布的貝葉斯分析。與已有文獻(xiàn)不同，本文將基于左刪失樣本，在部分加速壽命試驗(yàn)?zāi)Ｐ拖拢謩e用經(jīng)典方法和貝葉斯方法討論逆Rayleigh分布的統(tǒng)計(jì)性質(zhì)，并比較這兩種方法所得到的參數(shù)估計(jì)值。

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

在應(yīng)力水平S0，S1(S0＜S1)下進(jìn)行壽命試驗(yàn)，其中S0為正常應(yīng)力水平，S1為加速應(yīng)力水平。基本假設(shè)如下：

A 1 在應(yīng)力水平S0,S1下，受試產(chǎn)品有相同的失效機(jī)理。

A2 在應(yīng)力水平S0下，產(chǎn)品的失效時(shí)間為t；在應(yīng)力水平S1下，產(chǎn)品的失效時(shí)間為x，滿足t=βx，其中，β(β＞1)為加速因子，且t與x相互獨(dú)立。

A3 在應(yīng)力水平S0下，產(chǎn)品的失效時(shí)間t服從逆Rayleigh分布，其分布函數(shù)和概率密度函數(shù)分別為

其中，θ＞0為尺度參數(shù)。

根據(jù)假設(shè)A 1～A 3，在應(yīng)力水平S1下，產(chǎn)品的失效時(shí)間x服從逆Rayleigh分布，其分布函數(shù)和概率密度函數(shù)分別為

即在應(yīng)力水平S1下，產(chǎn)品的失效時(shí)間x服從參數(shù)為λ的逆Rayleigh分布。

現(xiàn)有n個(gè)受試產(chǎn)品被隨機(jī)分成2組，它們的容量分別為n0,n1。其中n0個(gè)樣品被安排到正常應(yīng)力水平S0下進(jìn)行壽命試驗(yàn)，n1個(gè)樣品被安排到加速應(yīng)力水平S1下進(jìn)行壽命試驗(yàn)。在應(yīng)力水平Si(i=0,1)下進(jìn)行左刪失試驗(yàn)：對(duì)ni(i=0,1)個(gè)產(chǎn)品在0時(shí)刻進(jìn)行壽命試驗(yàn)，直到時(shí)刻τi(i=0,1)時(shí)進(jìn)行觀測(cè)，此時(shí)發(fā)現(xiàn)已經(jīng)有di(i=0,1)個(gè)產(chǎn)品失效，但不知其失效時(shí)刻，此后進(jìn)行跟蹤觀測(cè)，可得到剩余的mi=ni-di(i=0,1)個(gè)產(chǎn)品的失效時(shí)刻。設(shè)觀測(cè)到的失效時(shí)刻分別為t=(t1,t2,…,tm0)（在應(yīng)力水平S0下）和x=(x1,x2,…,xm1)（在應(yīng)力水平S1下）。

2 極大似然估計(jì)

由上述試驗(yàn)數(shù)據(jù)，在應(yīng)力水平S0下的似然函數(shù)為

其中，C0＞ 0，且與θ無(wú)關(guān)。

在應(yīng)力水平S1下的似然函數(shù)為

其中，C1＞0，且與λ無(wú)關(guān)。

由于在不同應(yīng)力水平下進(jìn)行的壽命試驗(yàn)是相互獨(dú)立的，因此在左刪失恒定應(yīng)力部分加速壽命試驗(yàn)下的全似然函數(shù)為

用對(duì)數(shù)似然函數(shù)ln L求關(guān)于θ,λ的偏導(dǎo)數(shù)，得到似然方程組：

解得參數(shù)θ,λ的極大似然估計(jì)，分別為

由此，β的極大似然估計(jì)為

在式（9）中，令 d0=d1=0，m0=n0,m1=n1，不難得到在完全觀測(cè)樣本t=(t1,t2,…,tn0)和x=(x1,x2,…,xn1)下，參數(shù) θ,λ,β的極大似然估計(jì)分別為

3 Fisher信息矩陣及近似置信區(qū)間

利用LOUIS[12]提出的缺失信息原則，計(jì)算觀測(cè)到的Fisher信息矩陣。缺失信息原則的概念可表述為觀測(cè)信息=完全信息-遺失信息。

記 η=(θ,β)，W 表示完全數(shù)據(jù)，X 表示觀測(cè)數(shù)據(jù)，IW(η)為完全信息矩陣，IX(η)為觀測(cè)信息矩陣，IW|X(η)為遺失信息矩陣，則觀測(cè)信息矩陣為

在完全數(shù)據(jù)情形下計(jì)算得到的完全信息矩陣為

在應(yīng)力水平S0下，記Y=(Y1,Y2,…,Yd0)，則Y表示受試樣品在區(qū)間(0,τ0)內(nèi)失效時(shí)刻構(gòu)成的向量，由于無(wú)具體的值，因此視為遺失數(shù)據(jù)。同理，在應(yīng)力水平S1下，記Z=(Z1,Z2,…,Zd1)，且Z表示受試樣品在區(qū)間(0,τ1)內(nèi)失效時(shí)刻構(gòu)成的向量。

根據(jù)條件密度公式，得到Y(jié)j(j=1,2,…,d0)的概率密度函數(shù)為

遺失信息矩陣為

在應(yīng)力水平 S0下，觀測(cè)區(qū)間 (0,τ0)的 Fisher信息矩陣為

在應(yīng)力水平 S1下，觀測(cè)區(qū)間 (0,τ1)的 Fisher信息矩陣為

由文獻(xiàn)[13]知，當(dāng)樣本量很大時(shí)，參數(shù)的極大似然估計(jì) (θ?M,β?M)漸近服從二元正態(tài)分布，且期望為(θ,β)，方差、協(xié)方差矩陣為I-1X(η?)（觀測(cè)信息矩陣IX(η?)的逆矩陣），于是有

因此，參數(shù)θ,β的置信水平為100(1-γ)%的近似置信區(qū)間，分別為

其中，U11和U22為方差、協(xié)方差矩陣的主對(duì)角線元素，Zγ/2為標(biāo)準(zhǔn)正態(tài)分布的γ/2上分位數(shù)。

4 參數(shù)的貝葉斯估計(jì)

當(dāng)有先驗(yàn)信息可以利用時(shí)，采用貝葉斯方法通常可以提高估計(jì)的精度，因此，本節(jié)將計(jì)算參數(shù)θ和加速因子β的貝葉斯估計(jì)。

取參數(shù)θ的先驗(yàn)密度為

其中a0＞0為超參數(shù)。由式（7）和式（11）及貝葉斯公式，可得θ的后驗(yàn)密度為

在平方損失函數(shù)下，θ的貝葉斯估計(jì)為

由于在應(yīng)力水平S0下，產(chǎn)品的可靠度函數(shù)為

因此，在平方損失函數(shù)下R0(t)的貝葉斯估計(jì)為

取參數(shù)λ的先驗(yàn)密度為

其中a1＞0為超參數(shù)。由式（8）和式（15），可得λ的后驗(yàn)密度：

由式（17）無(wú)法求出顯式解，可用數(shù)值積分求解。

下面用極大似然法估計(jì)超參數(shù)a0,a1。

在應(yīng)力水平S0下，根據(jù)先驗(yàn)分布π0(θ)，可得到經(jīng)驗(yàn)分布：

以 Ft(τ0),ft(ti)分 別 代 替 上 式 中 的Ft(τ0;θ),[ft(ti;θ )]，則有

記式（19）的最終迭代值為 a?0。

用類似的方法可求得超參數(shù)a1的估計(jì)值，記為a?1。將 a?0,a?1代入式（13）、（14）、（17），就可求得相應(yīng)的貝葉斯估計(jì)值。

5 Monte Carlo模擬

對(duì)于不同的 n0,n1,θ,λ，從均方誤差（MSE）角度考慮分布參數(shù)和加速因子的估計(jì)性質(zhì)。模擬過(guò)程如下：

第1步給定n0,n1及參數(shù)θ,λ的值。

第2步產(chǎn)生一個(gè)容量為n且服從均勻分布U(0,1)的獨(dú)立同分布樣本U1,U2,…,Un，令

ti=(-θ/ln Ui)1/2，i=1,2,…,n0，

xi=(-λ/ln Ui)1/2,i=n0+1,n0+2,…,n0+n1，則(t1,t2,…,tn0)(xn0+1,xn0+2,…,xn0+n1)就是分別來(lái)自于逆Rayleigh分布式（1）和式（5）的樣本。

第3步給定τ0,τ1的值，得到左刪失樣本。

第4步基于上面得到的左刪失樣本，分別計(jì)算θ,β的點(diǎn)估計(jì)。

第5步步驟2～步驟4重復(fù)1 000次，計(jì)算θ,β的均方誤差。

第 6步在不同的 n0,n1,θ,λ,τ0,τ1下，重復(fù)上述過(guò)程。

將以上步驟產(chǎn)生的隨機(jī)模擬結(jié)果列于表1，從表1中可以看到，隨著n0,n1的增加，參數(shù)θ,β的極大似然估計(jì)和貝葉斯估計(jì)的均方誤差逐漸減小，說(shuō)明極大似然估計(jì)具有大樣本性質(zhì)。在小樣本場(chǎng)合，貝葉斯估計(jì)要明顯優(yōu)于極大似然估計(jì)，隨著樣本量的增加，兩者的差距逐漸縮小。為了比較左刪失樣本下和完全觀測(cè)樣本下參數(shù)估計(jì)精度的差異，當(dāng)θ=2，λ=1時(shí)，模擬計(jì)算了不同樣本容量下參數(shù)θ,β估計(jì)的均方誤差（見(jiàn)表2），結(jié)果顯示，完全觀測(cè)樣本下參數(shù)估計(jì)的精度要優(yōu)于左刪失樣本。

為了了解在具體樣本下各種估計(jì)的效果，取θ=3，λ=2，通過(guò)隨機(jī)模擬的方式分別產(chǎn)生容量為40和50的服從逆Rayleigh分布的樣本，在應(yīng)力水平S0下（θ=3），按從小到大的順序排列：

0.914 ，1.017，1.102，1.109，1.148，1.194，1.277，1.309，1.341，1.381，1.416，1.502，1.515，1.673，1.674，1.687，1.713，1.834，1.875,1.922，2.125，2.211，2.245，2.427，2.617，2.859，2.918，3.022，3.062，3.489，3.567，3.571，3.601，3.725，3.787，3.935，4.965，5.098，5.633，6.288。

在應(yīng)力水平S1下（λ=2），按從小到大的順序排列：0.784，0.793，0.887，0.938，1.056，1.093，1.098，1.108，1.109，1.144，1.184，1.192，1.235，1.253，1.255，1.261，1.263，1.382，1.389，1.401，1.410，1.427，1.451，1.485，1.505，1.513，1.515，1.538，1.546，1.608，1.643，1.673，1.767，1.787，1.884，1.949，2.014，2.192，2.368，2.552，2.792，3.012，3.042，3.065，4.501，4.841，5.121，5.758，6.142，8.505。

表1 左刪失樣本下參數(shù)估計(jì)的均方誤差Table 1 Mean square errors of parameter estimation under left censored samples

表2 完全觀測(cè)樣本下參數(shù)估計(jì)的均方誤差Table 2 Mean square errors of parameter estimation under complete observation samples

表3 參數(shù)的點(diǎn)估計(jì)及置信區(qū)間Table 3 Point estimation and confidence interval of the parameters

由τ0，τ1得到左刪失樣本，計(jì)算θ,β的極大似然估計(jì)、貝葉斯估計(jì)及置信水平為95%的置信區(qū)間，并分別記為θ?M，θ?B，(θ?L,θ?U)和β?M,β?B，(β?L,β?U)，相應(yīng)的估計(jì)結(jié)果見(jiàn)表3。

由于參數(shù)真值θ=3,β=1.224 7，從表3的數(shù)據(jù)中可知，參數(shù)的貝葉斯估計(jì)比極大似然估計(jì)更接近參數(shù)真值，參數(shù)真值位于置信區(qū)間內(nèi)。通常，隨著左刪失數(shù)據(jù)的增多，參數(shù)估計(jì)值與真值之間的差值以及置信區(qū)間的長(zhǎng)度呈逐漸變大趨勢(shì)。

6 結(jié) 語(yǔ)

加速壽命試驗(yàn)是利用產(chǎn)品壽命和應(yīng)力水平之間的某種已知函數(shù)關(guān)系，估算產(chǎn)品在正常應(yīng)力水平時(shí)的壽命。當(dāng)產(chǎn)品壽命和應(yīng)力水平之間的關(guān)系未知時(shí)，可采用恒定應(yīng)力部分加速壽命試驗(yàn)，此方法更適合新產(chǎn)品。用兩種應(yīng)力水平下的試驗(yàn)數(shù)據(jù)估計(jì)加速因子，然后將高應(yīng)力水平下的壽命轉(zhuǎn)換為正常應(yīng)力水平時(shí)的壽命。近年來(lái)，部分加速壽命試驗(yàn)已被應(yīng)用于可靠性統(tǒng)計(jì)推斷。在本文中，當(dāng)產(chǎn)品的失效時(shí)間服從逆Rayleigh分布時(shí)，基于恒定應(yīng)力部分加速壽命試驗(yàn)數(shù)據(jù)，得到了分布參數(shù)和加速因子的極大似然估計(jì)及貝葉斯估計(jì)。由遺失信息原則計(jì)算極大似然估計(jì)的漸近方差，得到參數(shù)的近似置信區(qū)間。通過(guò)數(shù)值模擬說(shuō)明樣本量對(duì)估計(jì)精度的影響，并將極大似然估計(jì)與貝葉斯估計(jì)做了比較，結(jié)果表明，在小樣本場(chǎng)合，貝葉斯估計(jì)要優(yōu)于極大似然估計(jì)。