秦 睿

（廣西師范大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，廣西桂林541004）

0 引言

在產(chǎn)品的可靠性試驗(yàn)，生物醫(yī)學(xué)實(shí)驗(yàn)的壽命評(píng)估，工程學(xué)的疲勞可靠性分析等領(lǐng)域中，受到試驗(yàn)時(shí)間和成本等因素的影響，大量試驗(yàn)產(chǎn)生了不完整數(shù)據(jù)，實(shí)際生活如在臨床醫(yī)學(xué)研究中，由于研究成本、研究對(duì)象失訪或中途退出研究、研究對(duì)象因其他情況死亡等因素，造成研究人員無(wú)法觀察所有研究對(duì)象完整的試驗(yàn)過(guò)程，出現(xiàn)收集到不完整數(shù)據(jù)的情形。例如，定時(shí)截尾、定數(shù)截尾、混合截尾等都是不完整數(shù)據(jù)情形，其中多重Ⅰ型混合截尾情形在實(shí)際生活中廣泛存在。

為了解產(chǎn)品的可靠程度，工程學(xué)的材料耐疲勞程度等各項(xiàng)指標(biāo)，研究人員需要從上述不完整數(shù)據(jù)情形中統(tǒng)計(jì)分析相應(yīng)的結(jié)果。不同于完整數(shù)據(jù)樣本的統(tǒng)計(jì)分析，截尾樣本需要提前對(duì)樣本數(shù)據(jù)進(jìn)行處理后，選擇適當(dāng)?shù)慕y(tǒng)計(jì)方法才能得到比較可靠的結(jié)果，為研究人員的決策提供依據(jù)。

現(xiàn)已有大批專(zhuān)家學(xué)者對(duì)Ⅰ型混合截尾樣本進(jìn)行了大量的研究。在Ⅰ型混合截尾方案測(cè)試中有n個(gè)壽命是獨(dú)立同分布的隨機(jī)變量x，當(dāng)測(cè)試中有m(m≤n)個(gè)單元失效或者測(cè)試時(shí)間達(dá)到預(yù)先設(shè)定的測(cè)試停止時(shí)間T，則整個(gè)測(cè)試停止，此時(shí)分別有如下兩種數(shù)據(jù)情況，第一種當(dāng)測(cè)試中有m(m≤n)個(gè)單元失效，則能夠觀察到m個(gè)完整的壽命數(shù)據(jù)為x1,n≤x2,n≤…≤xm,n，測(cè)試停止時(shí)間為xm,n；第二種當(dāng)測(cè)試時(shí)間達(dá)到預(yù)先設(shè)定時(shí)間T，則能夠觀察到d個(gè)完整的壽命數(shù)據(jù)為x1,n≤x2,n≤…≤xd,n，其中d＜m,xd,n≤T＜xd+1,n，測(cè)試停止時(shí)間為T(mén)。整個(gè)測(cè)試系統(tǒng)的停止時(shí)間定義為min{xm,n,T}。Epstein(1954)在文獻(xiàn)[1]中提出Ⅰ型混合截尾方案次序統(tǒng)計(jì)量的應(yīng)用。在分布已知的情況下，陳家鼎(1989)在文獻(xiàn)[2]證明了基于隨機(jī)截尾數(shù)據(jù)的威布爾分布參數(shù)最大似然估計(jì)具有強(qiáng)相合性；劉煥彬(1994)在文獻(xiàn)[3]中證明了在Ⅰ型混合截尾情況下最大似然估計(jì)的強(qiáng)相合性。Bedbur(2019)等在文獻(xiàn)[4]中研究了多重序貫次序樣本下威布爾分布參數(shù)的統(tǒng)計(jì)推斷，證明了最大似然估計(jì)的一致性和漸近有效性。Górny(2020)等在文獻(xiàn)[5]中討論了在多重Ⅰ型混合刪失序貫樣本情況下指數(shù)分布的最大似然估計(jì)，并用實(shí)際樣本數(shù)據(jù)證明了估計(jì)的有效性。目前針對(duì)多重Ⅰ型混合截尾樣本情形的參數(shù)估計(jì)鮮有人研究。因此，本文利用最大似然估計(jì)方法對(duì)多重Ⅰ型混合截尾樣本情形下的兩參數(shù)威布爾(Weibull)分布參數(shù)進(jìn)行估計(jì)。

1 參數(shù)估計(jì)

假設(shè)壽命隨機(jī)變量X服從一個(gè)比例參數(shù)和形狀參數(shù)分別為λ和c的兩參數(shù)威布爾分布，其關(guān)于隨機(jī)變量X的概率密度函數(shù)(PDF)和累積分布函數(shù)(CDF)分別為

其中，比例參數(shù)(scale parameter)λ≥0 和形狀參數(shù)(shape parameter)c≥0。

多重Ⅰ型混合截尾方案是指，假設(shè)測(cè)試人員能夠觀察到試驗(yàn)樣本的壽命為Xij,1≤i≤k,1≤j≤mi，其中設(shè)測(cè)試系統(tǒng)有k個(gè)測(cè)試組，第i個(gè)測(cè)試組有ni個(gè)試驗(yàn)樣本，設(shè)定第i個(gè)測(cè)試組的固定失效數(shù)為mi(mi≤ni)和測(cè)試停止時(shí)間為T(mén)i＞0，1≤i≤k。則多重Ⅰ型混合截尾的第i組第j個(gè)次序統(tǒng)計(jì)量為

第i組測(cè)試停止時(shí)間為

第i組測(cè)試觀察到失效個(gè)數(shù)為

其中I(-∞,Ti](x)表示x的示性函數(shù)，即x∈(-∞,Ti]示性函數(shù)為1，否則為0。整個(gè)測(cè)試觀察到失效的總個(gè)數(shù)為

參考Kamps U(1995)在文獻(xiàn)[6]中提出的次序統(tǒng)計(jì)的聯(lián)合密度函數(shù)，則關(guān)于多重Ⅰ型混合截尾樣本下威布爾分布的參數(shù)λ和c的似然函數(shù)定義為

將兩參數(shù)威布爾分布的PDF(1)式和CDF(2)式帶入(3)式得到似然函數(shù)

為求解當(dāng)(4)式達(dá)到最大時(shí)的參數(shù)λ和c的估計(jì)值?和?，需要先對(duì)(4)式取對(duì)數(shù)得

其中ai是獨(dú)立于參數(shù)λ和c的常數(shù)。

把對(duì)數(shù)似然函數(shù)(5)式分別關(guān)于參數(shù)λ和c求偏導(dǎo)并令為0 得下式

因此，最大似然估計(jì)?和?只需滿足下式即可

通過(guò)(6)式可以看出，最大似然估計(jì)?和?沒(méi)有具體的顯示表達(dá)式，因此在計(jì)算估計(jì)值上采用簡(jiǎn)單有效的牛頓-拉弗森方法(Newton-Raphson method)迭代求解(6)式方程組。

設(shè)

則牛頓-拉弗森方法迭代方程為

其中

具體參數(shù)λ和c的牛頓-拉弗森方法迭代算法描述如下，假設(shè)預(yù)先給定一個(gè)非常小的常數(shù)γ，本文給定γ=10-8，

(1)給定迭代次數(shù)l=0, 初始值λ0和c0；

(2)根據(jù)(7)式計(jì)算并更新cl+1；

(3)根據(jù)(7)式計(jì)算并更新λl+1；

(4)如果|cl+1-cl| ＜γ成立，則停止迭代，輸出參數(shù)最大似然估計(jì)值

(5)否則重復(fù)步驟(2), (3), (4), (5)。

最大似然估計(jì)的精度可以用偏差和均方誤差來(lái)表示，也能展現(xiàn)出估計(jì)值的穩(wěn)定性，但不夠直觀。置信區(qū)間能簡(jiǎn)單和直觀地給出參數(shù)估計(jì)值的精度和可靠程度。本文基于Thoman D R(1969)等文獻(xiàn)[7]中的定理A 和定理B 構(gòu)造參數(shù)的精確置信區(qū)間。

在參數(shù)λ=1 和c=1 的威布爾分布情形下，設(shè)參數(shù)λ和c的最大似然估計(jì)分別為和。根據(jù)文獻(xiàn)[7]定理Ac獨(dú)立于參數(shù)λ和c的分布，且與具有相同的分布。其中?為多重Ⅰ型混合截尾樣本下威布爾分布參數(shù)c的最大似然估計(jì)。因此，可用的分布來(lái)構(gòu)造未知參數(shù)c的置信區(qū)間。

設(shè)l1,l2為任意兩個(gè)常數(shù)，使得下式在置信水平為1-α的條件下成立

由于與同分布，將=代入(8)式可 得，

最終可得，置信水平為1-α的威布爾分布參數(shù)c的置信區(qū)間為

設(shè)t1,t2為任意兩個(gè)常數(shù)，使得下式在置信水平為1-α的條件下成立

最終可得，置信水平為1-α的威布爾分布參數(shù)λ的置信區(qū)間為

2 模擬研究與實(shí)例分析

假設(shè)兩參數(shù)威布爾分布參數(shù)真值分別為λ=2 和c=2，那么對(duì)應(yīng)參數(shù)估計(jì)的偏差(BIAS)和均方誤差(MSE)見(jiàn)表1。

表1 展示的是樣本總量分別為50, 100, 200, 500 的比例參數(shù)λ和形狀參數(shù)c都為2 時(shí)的兩參數(shù)威布爾分布，在b= 3 5 的多重Ⅰ型混合截尾方案下參數(shù)最大似然估計(jì)的偏差(BIAS)和均方誤差(MSE)。根據(jù)表1 可以看出，最大似然估計(jì)的偏差和均方誤差都隨樣本總量的增大而減小，說(shuō)明參數(shù)的最大似然估計(jì)表現(xiàn)出良好的穩(wěn)定性和有效性。

表1 最大似然估計(jì)的偏差(BIAS)和均方誤差(MSE)

表2 展示的是分組比例分別取1/2, 3/5, 4/5 中不同值時(shí)，比例參數(shù)λ和形狀參數(shù)c都為2 時(shí)的雙參數(shù)威布爾分布，在樣本總量為100 的多重Ⅰ型混合截尾方案下參數(shù)最大似然估計(jì)的偏差(BIAS)和均方誤差(MSE)。根據(jù)表2 可以看出，最大似然估計(jì)的偏差和均方誤差均隨分組比例的增大而減小，說(shuō)明最大似然估計(jì)是符合實(shí)際情形和邏輯，可應(yīng)用性較好。

表2 b 取不同值時(shí)參數(shù)估計(jì)的偏差(BIAS)和均方誤差(MSE)

實(shí)際案例分析數(shù)據(jù)來(lái)源于Proschan(1963)在文獻(xiàn)[8]給出并討論的13 架波音飛機(jī)空調(diào)系統(tǒng)的故障時(shí)間間隔（小時(shí)）數(shù)據(jù)，規(guī)定飛機(jī)完整運(yùn)行2000 小時(shí)后，必須對(duì)飛機(jī)進(jìn)行保養(yǎng)維護(hù)工作。Proschan(1963)對(duì)數(shù)據(jù)進(jìn)行分析發(fā)現(xiàn)每架飛機(jī)空調(diào)系統(tǒng)的故障時(shí)間間隔分布呈形狀參數(shù)c為1 的威布爾分布。

我們考慮使用飛機(jī)編號(hào)為7808 至7814 的7 架飛機(jī)空調(diào)系統(tǒng)故障時(shí)間間隔觀察值數(shù)據(jù)，對(duì)數(shù)據(jù)進(jìn)行處理后應(yīng)用多重Ⅰ型混合截尾方案，即在整個(gè)過(guò)程中，共有k=7 組數(shù)據(jù)，設(shè)定每組觀察都觀察25 次，即ni=25,1≤i≤k，規(guī)定飛機(jī)在保養(yǎng)維護(hù)周期內(nèi)出現(xiàn)14 次及以上空調(diào)故障，就必須進(jìn)行保養(yǎng)維護(hù)工作，即固定觀察失效數(shù)量mi=14,1≤i≤k，固定停止時(shí)間為T(mén)i=2000,1≤i≤k。則通過(guò)實(shí)際觀察值可以得到如下變量的數(shù)值，每架飛機(jī)保養(yǎng)維護(hù)周期的運(yùn)行時(shí)間（停止時(shí)間）為

每架飛機(jī)保養(yǎng)維護(hù)周期內(nèi)空調(diào)故障次數(shù)（失效個(gè)數(shù)）為

利用最大似然估計(jì)方法對(duì)多重Ⅰ型混合截尾樣本情形下的兩參數(shù)威布爾分布參數(shù)進(jìn)行估計(jì)，得到比例參數(shù)的最大似然估計(jì)?=1952.12 和形狀參數(shù)的最大似然估計(jì)?=1.05，其分布的總體均值為1999.07。

從得到的結(jié)果可以看出，將本文的估計(jì)方法運(yùn)用到實(shí)際案例中，得到的結(jié)果是良好，可接受。前人學(xué)者已證明了上述數(shù)據(jù)服從形狀參數(shù)c為1 的威布爾分布，此處的形狀參數(shù)估計(jì)值為?=1.05，兩者相差很小，也符合數(shù)值模擬研究給出的結(jié)果。此外，通過(guò)比例參數(shù)和形狀參數(shù)的最大似然估計(jì)計(jì)算出的總體均值為1999.07，表示在保養(yǎng)維護(hù)周期內(nèi)飛機(jī)運(yùn)行時(shí)間為1999.07 小時(shí)，即飛機(jī)空調(diào)故障時(shí)間滿足飛機(jī)運(yùn)行2000 小時(shí)強(qiáng)制保養(yǎng)維護(hù)規(guī)定，保證了飛機(jī)運(yùn)行時(shí)的安全性和可控性，說(shuō)明方法是可應(yīng)用性良好，符合實(shí)際情況。

3 結(jié)論

本文探究了在多重Ⅰ型混合截尾樣本情形下，兩參數(shù)威布爾分布的參數(shù)估計(jì)問(wèn)題。將多重Ⅰ型混合截尾樣本進(jìn)行次序統(tǒng)計(jì)量處理后，得到次序統(tǒng)計(jì)量的聯(lián)合密度函數(shù)，進(jìn)而使用最大似然估計(jì)方法計(jì)算得到參數(shù)估計(jì)，利用牛頓-拉弗森迭代法求解出參數(shù)估計(jì)值，并給出參數(shù)的精確置信區(qū)間。數(shù)值模擬研究與實(shí)例分析結(jié)果表明，基于最大似然估計(jì)方法能快速得出參數(shù)的估計(jì)和精確置信區(qū)間，同時(shí)參數(shù)估計(jì)值表現(xiàn)出良好的穩(wěn)定性和可應(yīng)用性。此外，在求解參數(shù)估計(jì)值時(shí)，使用了牛頓-拉弗森迭代法，可能會(huì)使參數(shù)估計(jì)值陷入局部最優(yōu)解中，不能得到更加準(zhǔn)確的參數(shù)估計(jì)值，后期需要在參數(shù)估計(jì)值求解的方法上進(jìn)一步加強(qiáng)。另外，后續(xù)可以考慮使用貝葉斯估計(jì)等其他方法估計(jì)參數(shù)。