龍兵,張忠占

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

1.引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.極大似然估計(jì)

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

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

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

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

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

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

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

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

在應(yīng)力水平S0下,記Y=(Y1,Y2,··· ,Yk),Yj=(Yj1,Yj2,··· ,Yjdj),且Yj(j=1,2,··· ,k)表示受試元件在區(qū)間[Tj?1,Tj) 內(nèi)失效時(shí)刻構(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)失效時(shí)刻構(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)式兩邊取對(duì)數(shù)得

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

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

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

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

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

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

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

由于

因此

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

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

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

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

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

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

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

遺失信息矩陣為

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

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

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

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

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

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

5.數(shù)值模擬

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

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

步1 給定n,r的值;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6.結(jié)論

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