基于Xgamma分布恒加試驗的統計分析

張 鋒,戴林送

(安慶師范大學 數理學院,安徽 安慶 246133)

0 引言

壽命試驗[1～3]是研究、分析和評價產品可靠性的重要方法。隨著時代的不斷進步,一些產品的使用壽命得到極大提高,時至今日,市場上高壽命的產品與日俱增,對產品壽命的驗證如果采用一般的試驗方法耗費時間會進一步增加,故而大多采用加速壽命試驗的方法。本文探討的是加速壽命試驗中的一種模型,恒定應力加速壽命試驗亦稱恒加試驗。同時出于對成本和時間的考慮,加速壽命試驗通常在所有測試元件完全失效之前終止試驗,這種方法被稱之為截尾試驗,其主要有兩大類,一類是定時截尾壽命試驗,即試驗運行到事先既定的時刻時立刻終止試驗;另一類是定數截尾壽命試驗,也是本文所探究的對象,將試驗運行到事先規定好的失效數量就終止試驗,這就使截尾試驗在成本損耗上相比于完全失效下的壽命試驗有較為明顯的優勢。

Xgamma分布[4]于2016年被Subhradev Sen 等人發現并提出,該分布為指數分布和伽馬分布的特殊有限混合,通過對各種數學性質、結構和生存特性的推導,發現在多數情況下,該分布比指數分布具有更大的靈活性,通過對臨床患者的真實數據研究分析發現,相比于指數模型,Xgamma模型更能擬合數據,在醫療等方面更具有普適性。吳清太[5]對Weibull分布在恒加試驗定數截尾條件下的情況進行了討論,利用Gauss-Markov定理對加速方程下未知系數的值進行估計,并對估計結果進行了模擬試驗,但基于運用Gauss-Markov定理的門檻較高,很多情形下需要采用其他方法。2021年,Elshahhat Ahmed等人[6]對Xgamma分布的性質進行了研究,使用Lindley近似和Metropolis-Hastings算法對存在自適應Ⅱ型漸進混合刪失數據的情況下分布中的未知參數和一些生存特征進行了估計。文獻[7]也進一步對Xgamma分布的可靠性特征于三種損失函數下進行貝葉斯估計,模擬的結果表明貝葉斯方法優于經典估計方法。但是,目前有關Xgamma分布的研究大多只是單獨從分布的性質和可靠性角度進行分析,未曾考慮該分布在結合加速壽命試驗時的情況。

為此,本文考慮在恒定應力加速壽命試驗時,定數截尾樣本下,Xgamma模型的統計分析。整體從五個方面進行:1)給出Xgamma分布在恒加試驗定數截尾情況下的基本假定和實施過程;2)從極大似然估計和Bayes估計兩個角度討論定數截尾條件下Xgamma分布的未知參數值;3)鑒于Bayes估計中存在復雜積分,采用了Lindley近似以及馬爾可夫鏈蒙特卡洛(MCMC)方法獲取未知參數的近似Bayes估計;4)通過MATLAB軟件對理論推導的未知參數估計進行數值模擬試驗;5)對第四步模擬估計的結果進行準確度分析。

1 基本假定與實施過程

假設1 假設在正常應力S0和加速應力水平S1

其分布函數為

其中θi大于0,i=0,1,…,k.

假設2[8～11]元件的平均壽命θi與試驗中的應力水平S之間服從以下加速方程,即

lnθi=a+bφ(S),i=0,1,…,k

(3)

其中a,b為參數,φ(S)是應力S的函數.

在元件壽命服從Xgamma分布情況下,定數截尾下恒加試驗的實施過程如下:

1)在正常應力S0和k個加速應力水平S1,S2,…,Sk下,加速應力水平之間滿足:S1

2)根據科學抽樣原理,從一批元件中隨機抽中n=(n1+n2+…+nk)個元件進入恒加試驗,在應力水平Si下放置ni個元件進行定數截尾試驗,試驗安排:1)當有第1個元件失效時,記錄時刻,并且將剩下的(ni-1)個未失效的元件接續進行試驗;2)當發現第2個元件失效時,記錄時刻,將剩下的(ni-2)個未失效的元件接續進入下一步試驗;3)不斷重復下去,直至第ri個測試元件失效時終止試驗,并將剩下的(ni-ri)個元件移離試驗。假定在應力水平Si下ri個失效,則樣本數據遵從

Xi∶1≤Xi∶2≤…≤Xi∶ri,i=1,2,…,k

2 極大似然估計

對于截尾試驗,應力Si下的似然函數Li為

可得在k個應力下的似然函數L為

(4)

l(θi)=lnL(θi)

(5)

將加速方程lnθi=a+bφ(S)代入(5)式,整理可得到含參數a,b的對數似然函數式,為了獲取參數的極大似然估計(MLE),對其分別關于a,b求一階導:

(6)

可推得(6)式中的各項為

(7)

3 貝葉斯估計

在對未知參數(a,b)進行貝葉斯統計分析之前,需要先選取未知參數的先驗分布,未知參數的先驗分布如何進行選取對整個估計有著舉足輕重的作用。

本文采取的是由未知參數(a,b)已有的相關信息獲取先驗分布,由于此時未知參數a,b皆大于0,因而選取Gamma分布作為未知參數的先驗信息分布,并假定兩個參數之間相互獨立,所以不妨假定本文中的未知參數a,b服從參數為(g,d)的Gamma分布,那么a,b的先驗密度為

(8)

根據貝葉斯定理,由先驗分布式和似然函數式推得后驗分布π(a,b|T)的表達式

取平方誤差損失函數,推出未知參數的Bayes估計,對參數(a,b)取g(a,b)=a與g(a,b)=b,得出后驗期望[8]為

令eQ(a,b)=π(a,b)p(T|a,b),則有

(9)

Q(a,b)=l(a,b)+ρ(a,b),l(a,b)=lnp(T|a,b),ρ(a,b)=lnπ(a,b)

由此可推得Lindley近似和MCMC算法.

3.1 Lindley近似算法

Lindley[10]于1980年提出了一種解決近似貝葉斯估計中出現積分的方法。使用該方法可獲得各種壽命模型參數的貝葉斯估計量。貝葉斯分析中的積分比率由下式得出。

(10)

其中Θ={θ1,θ2,…,θm},u(Θ)是關于Θ的函數,lnL(Θ)是對數似然函數,ρ(Θ)取先驗信息的對數,對于大樣本和某些正則條件,式(10)可以通過以下展開式近似

在本文中,由于m=1,因此式(11)可變化為

對參數a的Lindley估計

其中

對ρ(a,b)中的參數a求偏導,得到

同理可對參數b的Lindley估計:

其中:

此時將b′,b″,Lbbb,ρb,wbb一同代入E(b|T),即可以求出E(b|T).

3.2 MCMC算法

在Bayes算法中,還有一種對后驗期望E(a|T)進行估計的經典算法,即對未知參數采用蒙特卡洛馬爾科夫鏈的方法進行計算。由于文中后驗期望E(a|T)是復雜積分無法直接抽樣,因此這里用到了MCMC方法[11]的MH抽樣法。

令目標分布的概率密度函數為g(a),則有

從q分布中產生一個隨機數a*,有

故接受概率為

若

則a1=a*,否則令a1=a0,一直下去產生一組隨機樣本序列a0,a1,a2,a3,…,ai.

4 數值模擬

由此可得基于恒加試驗下Xgamma分布定數截尾條件下的兩種參數估計,現利用MATLAB軟件對估計進行有效性分析和精度評估。

具體步驟如下:

1)本文首先根據對a,b兩個未知參數各自的取值范圍設定初始值,這里取θ=4,a=0.49,b=0.38.

2)從Xgamma分布中生成n個隨機觀測值,本文分別取樣本容量為20、40、60、80、100、120、140、160、180和200,依據定數截尾恒加試驗的試驗原則,取在2個應力加速水平fi=2∶3下的兩組樣本情況進行試驗模擬,并且分別選取不同數量的截尾數,得到在各自情況下的截尾樣本。

3)根據獲取的截尾樣本,對未知參數進行近似估計。在MLE中,將2 000次反復迭代的結果均值作為未知參數a,b的最終估計結果;在Lindley算法中給定初值g1=2,d1=1.5,g2=5,d2=4.1,重復上述步驟并記錄結果;對于MCMC算法,本文給定初值g1=5,d1=3.3,g2=5,d2=4.1,對每組參數進行迭代運算2 000次,取鏈長為30 001,對2 000次迭代運算的第20 001到30 000個數據鏈求平均值作為迭代運算的結果,結果見表1.

表1 模擬試驗結果

4)記錄上述步驟結果并計算未知參數的相對偏差和相對均方誤差,結果見表2.并對a,b兩個參數的RBias結果繪制折線圖,使結果的呈現更加形象直觀。如圖1所示。

圖1 數據模擬對比分析

表2 模擬試驗結果分析

由試驗結果可以看出:比較Lindley、MLE、MCMC三種算法關于Xgamma分布參數的近似估計結果,通過觀察表1與表2的試驗數據,結合圖1中的相對偏差(RBias)與相對均方根誤差(RMSE)的折線趨勢,在本文的假定模型下Bayes估計的精確程度更好,尤其是Lindley近似推出的Bayes估計結果更加穩定,說明Bayes估計的結果是可靠的。

5 結束語

討論了當元件壽命服從Xgamma分布時,在定數截尾方案下恒加試驗的參數估計。利用最大似然估計和貝葉斯估計方法對未知參數進行了推斷,利用Bayes估計能夠減少繁雜的積分問題,進而提升參數估計的精確程度。結果表明,使用貝葉斯估計中的Lindley近似和MCMC算法都能較極大似然估計獲得更加精確的結果。