恒定應(yīng)力部分加速壽命試驗的統(tǒng)計分析

龍兵,張忠占

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

1.引言

隨著科學(xué)技術(shù)的進步,出現(xiàn)了越來越多的高質(zhì)量、高可靠性的產(chǎn)品.在正常條件下不容易得到這些產(chǎn)品的使用壽命,為了縮短壽命或加速性能退化,所有的測試單元或它們中的一部分都受到比通常情況更嚴酷的應(yīng)力,該試驗被稱為加速壽命試驗.如果只放一部分,那么這個試驗叫做部分加速壽命試驗.近年來已經(jīng)有很多文獻研究了加速壽命試驗的統(tǒng)計分析問題.文[1]基于步進應(yīng)力加速壽命試驗?zāi)Ｐ?在逐步混合截尾樣本下給出了Weibull分布中參數(shù)及加速因子的極大似然估計和區(qū)間估計.文[2]在逐步I 型區(qū)間刪失樣本下,討論了k個應(yīng)力水平的步進應(yīng)力加速壽命試驗,得到了失效時間服從指數(shù)型分布模型的參數(shù)估計,從四個標準出發(fā)給出了最優(yōu)試驗方案.文[3] 在產(chǎn)品的的壽命服從Burr XII分布時,基于步進應(yīng)力部分加速壽命試驗?zāi)Ｐ?給出了模型參數(shù)和加速因子的極大似然估計,并得到了參數(shù)的方差、協(xié)方差矩陣.文[4]在恒定應(yīng)力部分加速壽命試驗下,推導(dǎo)出了Weibull分布參數(shù)的極大似然估計,并設(shè)計出最優(yōu)試驗方案.文[5]在逐步II型混合截尾恒加壽命試驗下,給出了指數(shù)分布競爭失效產(chǎn)品未知參數(shù)及可靠性指標的極大似然估計和貝葉斯估計.文[6]基于定時截尾樣本,討論了Pareto分布在恒定應(yīng)力部分加速壽命試驗下的最優(yōu)方案.另外,文[7-12] 也研究了多種壽命分布的可靠性統(tǒng)計分析問題.

加速壽命試驗的特點是把產(chǎn)品放在高于正常應(yīng)力水平下進行試驗,在獲得失效數(shù)據(jù)后,利用產(chǎn)品壽命和應(yīng)力水平之間的關(guān)系得到產(chǎn)品在正常應(yīng)力水平下的壽命特征.實際上,產(chǎn)品壽命和應(yīng)力水平之間的關(guān)系并不一定總是已知的.在這種情況下,可考慮部分加速壽命試驗.部分加速壽命試驗的統(tǒng)計分析不需要加速模型,它通過加速因子,將高應(yīng)力水平下的壽命轉(zhuǎn)化為正常應(yīng)力水平之下,進而進行統(tǒng)計分析.近年來部分加速壽命試驗已經(jīng)被應(yīng)用于可靠性統(tǒng)計推斷中,在試驗方案中一般采用定時截尾、定數(shù)截尾、逐步增加定數(shù)截尾和混合截尾等,對于這些傳統(tǒng)的截尾試驗,理論研究體系比較成熟,相關(guān)成果較多.

本文根據(jù)定時截尾和區(qū)間刪失試驗方案,提出一種新的壽命試驗方案,在本文中被稱為定時區(qū)間刪失試驗.基于試驗樣本對Pareto分布的可靠性問題進行統(tǒng)計分析,并進行了隨機模擬.

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

在2個應(yīng)力水平S0

(I)元件失效只由1個失效機理引起.

(II)在應(yīng)力水平S0下,元件的失效時間服從Pareto分布,其分布函數(shù)、密度函數(shù)和失效率函數(shù)分別為

式中α,θ分別被稱為尺度參數(shù)和形狀參數(shù).

(III)在應(yīng)力水平S1下,元件的失效率函數(shù)為

式中λ(>1)是加速因子.

根據(jù)失效率函數(shù)和可靠度函數(shù)的關(guān)系,可得在應(yīng)力水平S1下的可靠度函數(shù)為

進而,可得在應(yīng)力水平S1下元件失效時間的分布函數(shù)和密度函數(shù)分別為

定時區(qū)間刪失恒定應(yīng)力部分加速壽命試驗?zāi)Ｐ涂擅枋鋈缦拢簭膎個具有上述失效時間的元件中,選取n0=n(1?r),(0

1)在應(yīng)力水平S0下,觀測得到m0個受試元件的失效時刻為t01≤t02≤··· ≤t0m0,這里m0為隨機變量,且假設(shè)m0>0,此后每間隔一段時間進行1次觀測,在觀測時刻Tj發(fā)現(xiàn)有dj個元件在區(qū)間[Tj?1,Tj)(j=1,2,··· ,k)內(nèi)失效.

2)在應(yīng)力水平S1下,觀測得到m1個受試元件的失效時刻為t11≤t12≤··· ≤t1m1,假設(shè)m1>0,此后每間隔一段時間進行1次觀測,在觀測時刻Tj發(fā)現(xiàn)有rj個元件在區(qū)間[Tj?1,Tj)(j=1,2,··· ,k)內(nèi)失效.

定時截尾試驗是可靠性分析中的一種重要壽命試驗方案,它是把一批受試元件跟蹤觀測到預(yù)先設(shè)定的時間T0,其余沒有失效的元件全部撤離試驗,這樣可以得到在T0之前失效元件的具體壽命數(shù)據(jù)和在區(qū)間(T0,+∞)內(nèi)失效元件的個數(shù).如果在T0時刻撤離試驗的元件數(shù)所占的百分比較大,對估計的精度會產(chǎn)生較大的影響,而繼續(xù)進行跟蹤觀測可能面臨某種困難,有時甚至無法辦到,如果可以進行間斷觀測這樣就可以采用本文中的試驗方案.與定時截尾試驗相比較,就是把T0之后的一個區(qū)間刪失推廣到多個區(qū)間刪失,所得到的觀測信息更加準確,從而提高統(tǒng)計推斷的精度.下面將基于這類樣本進行統(tǒng)計分析.

3.極大似然估計

根據(jù)前面的模型和壽命試驗方案,在應(yīng)力水平S0下,忽略常量后失效樣本的似然函數(shù)為

這里I(·)為示性函數(shù).把(2.1-2.2)式代入(3.1)式得

在應(yīng)力水平S1下,忽略常量后失效樣本的似然函數(shù)為

把(2.6-2.7)式代入(3.3)式得

因此,同時在應(yīng)力水平S0和S1下進行壽命試驗,得到全似然函數(shù)為

由于L是關(guān)于α的單調(diào)遞增函數(shù),要使似然函數(shù)L的值最大,要取α= min(t01,t11),因此得到參數(shù)α的極大似然估計為

未知參數(shù)λ和θ的極大似然估計一般通過求解非線性方程組

得到,但計算較復(fù)雜,通常利用Newton-Raphson算法,可是在某些情況下該算法并不收斂,EM算法是處理缺損數(shù)據(jù)的一種較好方法,下面用EM算法來求參數(shù)的估計.

在應(yīng)力水平S0下,記Y=(Y1,Y2,··· ,Yk),Yj=(Yj1,Yj2,··· ,Yjdj),且Yj(j=1,2,··· ,k)表示受試元件在區(qū)間[Tj?1,Tj) 內(nèi)失效時刻構(gòu)成的向量,由于沒有具體的值,被視為遺失數(shù)據(jù).同理,在應(yīng)力水平S1下,記Z= (Z1,Z2,··· ,Zk),Zj= (Zj1,Zj2,··· ,Zjrj),且Zj(j=1,2,··· ,k)表示受試元件在區(qū)間[Tj?1,Tj)內(nèi)失效時刻構(gòu)成的向量.

根據(jù)條件密度公式,可以得到Y(jié)jl(l=1,2,··· ,dj) 的概率密度函數(shù)為

Zjl(l=1,2,··· ,rj) 的概率密度函數(shù)為

令V0= (t01,t02,··· ,t0m0),由V0和Y得W0= (V0,Y),則W0稱為在應(yīng)力水平S0下的偽完全數(shù)據(jù),基于偽完全數(shù)據(jù)得到的似然函數(shù)為

把(3.7)式兩邊取對數(shù)得

令V1=(t11,t12,··· ,t1m1),W1=(V1,Z),則W1稱為在應(yīng)力水平S1下的偽完全數(shù)據(jù),基于W1的似然函數(shù)為

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

同時在應(yīng)力水平S0和S1下進行壽命試驗,全似然函數(shù)為

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

E步：對(3.9)式求期望,得

M步：極大化Q(θ,λ),對Q(θ,λ)分別關(guān)于θ和λ求一階偏導(dǎo)數(shù),并令其為零可得

將方程(3.10)(3.11)聯(lián)立,可解得

由于

因此

根據(jù)(3.12)式,則可得迭代公式

其中θ(s)為θ的第s次迭代值.利用(3.14)式,經(jīng)過多次迭代可得到參數(shù)θ的最終估計值,記為.根據(jù)(3.13)式又可得參數(shù)λ的迭代公式

利用(3.15)式,經(jīng)過多次迭代又可得到參數(shù)λ的最終估計值,記為?λ.

4.Fisher信息矩陣及置信區(qū)間

Fisher信息矩陣反映總體分布的一種特性,它反映總體模型本身所提供的信息量.在參數(shù)估計理論研究中,起到相當(dāng)重要的作用.為了解決不完全數(shù)據(jù)樣本下的參數(shù)估計問題,Louis提出了遺失信息原則用于計算Fisher信息矩陣,其基本原則為:

記η= (θ,λ),W=完全數(shù)據(jù),X=觀測數(shù)據(jù),IW(η)=完全信息矩陣,IX(η)=觀測信息矩陣,IW|X(η)=遺失信息矩陣,則有觀測信息矩陣

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

遺失信息矩陣為

在應(yīng)力水平S0下,Fisher信息矩陣在觀測區(qū)間[Tj?1,Tj)計算為

在應(yīng)力水平S1下,Fisher信息矩陣在觀測區(qū)間[Tj?1,Tj)計算為

由于在大樣本下參數(shù)的極大似然估計漸近服從二元正態(tài)分布,期望為(θ,λ),方差、協(xié)方差矩陣為IX?1(),這里IX?1() 是觀測信息矩陣IX()的逆矩陣,且其中的α用極大似然估計量α?代替,于是有

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

醫(yī)學(xué)技術(shù)的發(fā)展讓外科手術(shù)的治療水平有了顯著提升,但外科手術(shù)患者通常會因為手術(shù)治療出現(xiàn)短暫的功能障礙,在缺乏家屬和護理人員陪伴的情況下易出現(xiàn)安全事故,因而對于此類患者的護理工作應(yīng)重點強化風(fēng)險管理,分析識別可能產(chǎn)生的風(fēng)險因素,以此展開針對性的護理工作。

這里U11和U22是方差、協(xié)方差矩陣的主對角線的元素,Zγ/2是標準正態(tài)分布的γ/2上分位數(shù).

5.數(shù)值模擬

在這一部分將進行模擬研究,論證估計問題的理論結(jié)果.利用絕對值相對偏差(RAB)和均方誤差(MSE)討論了形狀參數(shù)和加速因子估計結(jié)果的性質(zhì),這里

另外,得到了形狀參數(shù)和加速因子的置信水平為95%的近似置信區(qū)間.模擬過程如下:

步1 給定n,r的值;

步2 給定參數(shù)(α,θ,λ)=(10,2,1.5);

步3 產(chǎn)生一個容量為n(n=n0+n1)且服從均勻分布U(0,1)的獨立同分布樣本U1,U2,··· ,Un,令

步4 給定觀測點(T0,T1,T2,T3,T4)=(20,25,30,40,+∞),得到定時區(qū)間刪失樣本;

步5 基于上面得到的定時區(qū)間刪失樣本,分別利用(3.14)(3.15)式計算模型參數(shù)θ和λ的極大似然估計;

步6 把步3-5重復(fù)1000次,計算出參數(shù)的平均絕對值相對偏差和均方誤差;

步7 在不同的n和r下,重復(fù)上述過程,相關(guān)數(shù)據(jù)列于表5.1中;

步8 給定α=5,θ=2,λ=1.8,確定觀測點(T0,T1,T2,T3,T4)=(9,11,13,15,+∞),在不同的n和r下產(chǎn)生了12個定時區(qū)間刪失樣本,列于表5.2中;

步9 計算出參數(shù)θ和λ的點估計和置信水平為95%的近似置信區(qū)間,列于表5.3中.

表5.1 參數(shù)的絕對值相對偏差和均方誤差

表5.2 定時區(qū)間刪失樣本

表5.3 參數(shù)的點估計及95%置信區(qū)間

表5.4 不同應(yīng)力水平下可靠度的估計

從表5.1中可以看到,尺度參數(shù)α的極大似然估計非常接近真值,當(dāng)固定n時,隨著r的增加參數(shù)θ極大似然估計的相對偏差和均方誤差逐漸增大.當(dāng)r固定時,隨著n的增加參數(shù)θ和λ的極大似然估計的平均相對偏差和均方誤差都逐漸減小,這說明極大似然估計具有大樣本性質(zhì).從表5.3中可以看到,參數(shù)真值都位于置信區(qū)間內(nèi)部,隨著樣本量n逐漸增大,置信區(qū)間的長度相應(yīng)變小.

在可靠性理論中可靠度函數(shù)的估計是一個重要的問題,在大多數(shù)實際應(yīng)用和壽命試驗中,如Weibull、Burr XII、Pareto和指數(shù)分布是相當(dāng)合適的模型,已經(jīng)有許多文獻討論了這些分布可靠度函數(shù)的估計問題.下面將給出Pareto分布可靠度函數(shù)的估計,根據(jù)(2.1)式和(2.5)式可以分別得到在應(yīng)力水平S0和S1下可靠度函數(shù)的估計為

6.結(jié)論

由于很多元件有很長的使用壽命,在正常條件下通常需要很長的試驗時間.因此加速壽命試驗或者部分加速壽命試驗?zāi)軌蛟谳^短的時間內(nèi)得到失效數(shù)據(jù).在加速壽命試驗時,很多情況下也涉及到截尾,本文根據(jù)定時截尾試驗的弊端,提出了一種新的壽命試驗方案,即定時區(qū)間刪失試驗.根據(jù)試驗樣本對Pareto分布進行可靠性分析,利用Monte Carlo模擬,從平均絕對值相對偏差和均方誤差方面討論了估計的性質(zhì).另外,在定時區(qū)間刪失試驗下,對其它多參數(shù)壽命分布的統(tǒng)計分析也是一個值得研究的問題.