雙邊定時(shí)截尾下BS分布的統(tǒng)計(jì)分析

馬怡舟,徐曉嶺,顧蓓青

(上海對(duì)外經(jīng)貿(mào)大學(xué) a.國(guó)際經(jīng)貿(mào)學(xué)院; b.統(tǒng)計(jì)與信息學(xué)院，上海 201620)

在可靠性壽命試驗(yàn)中，產(chǎn)品失效有各種原因，其中有一種是由于在某種周期應(yīng)力作用下，產(chǎn)品產(chǎn)生裂縫，如果裂縫長(zhǎng)度達(dá)到或者超過(guò)某一值時(shí)，則產(chǎn)品失效。兩參數(shù)Birnbaum-Saunders模型是概率物理方法中的一個(gè)重要失效分布模型，Birnbaum和Sauders 1969年在文獻(xiàn)[1]中研究主因裂紋擴(kuò)展導(dǎo)致的材料失效過(guò)程中推導(dǎo)出來(lái)一個(gè)兩參數(shù)疲勞壽命分布，這一失效分布模式在機(jī)械產(chǎn)品可靠性研究中應(yīng)用廣泛，主要應(yīng)用于疲勞失效研究，對(duì)于電子產(chǎn)品性能退化失效分析也有重要應(yīng)用。該兩參數(shù)疲勞壽命分布的分布函數(shù)、密度函數(shù)分別為

t≥0,α>0,β>0

(1)

t≥0,α>0,β>0

(2)

其中，Φ(x),φ(x)分別為標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù)和密度函數(shù)。參數(shù)α和β分別稱為形狀參數(shù)和刻度參數(shù)，記此分布為BS(α,β)。

值得指出的是在國(guó)內(nèi)，文獻(xiàn)[2]根據(jù)混凝土在疲勞荷載下的損傷機(jī)理，提出以應(yīng)變做為表征其損傷程度的量度，認(rèn)為混凝土的疲勞損傷過(guò)程具有馬爾可夫性，通過(guò)求解FOKKER-PLANCK方程，導(dǎo)出了在指定時(shí)間下疲勞損傷分布服從兩參數(shù)Birnbaum Saunders疲勞壽命分布。利用文獻(xiàn)[3]所給出的在直接拉伸作用下混凝土抗拉疲勞壽命的試驗(yàn)數(shù)據(jù)資料，認(rèn)為用兩參數(shù)BS疲勞壽命分布BS(α,β)擬合比用對(duì)數(shù)正態(tài)分布更合適。

由于Birnbaum-Saunders疲勞壽命分布是從疲勞過(guò)程的基本特征出發(fā)導(dǎo)出的，因此它比常用壽命分布如威布爾分布,對(duì)數(shù)正態(tài)分布更適合描述某些由于疲勞而引起失效的產(chǎn)品壽命規(guī)律。此分布已經(jīng)成為可靠性統(tǒng)計(jì)分析中的常用分布之一。關(guān)于兩參數(shù)BS疲勞壽命分布BS(α,β)的統(tǒng)計(jì)分析已有眾多的文獻(xiàn)研究。文獻(xiàn)[4-5]說(shuō)明了在較文獻(xiàn)[1]更弱的條件下產(chǎn)品的壽命分布仍服從BS(α,β)，因此用BS(α,β)分布擬合這類產(chǎn)品的壽命分布比常用的Weibull、對(duì)數(shù)正態(tài)分布更合適。在完全樣本情形，文獻(xiàn)[6]討論了BS(α,β)分布參數(shù)的點(diǎn)估計(jì)，導(dǎo)出了參數(shù)的極大似然估計(jì)。文獻(xiàn)[7]利用模擬方法和極大似然估計(jì)的漸近正態(tài)性討論了參數(shù)的區(qū)間估計(jì)和假設(shè)檢驗(yàn)。文獻(xiàn)[8]研究了BS(α,β)分布的對(duì)數(shù)線性模型，導(dǎo)出了參數(shù)的極大似然估計(jì)和最小二乘估計(jì)，利用極大似然估計(jì)的漸近正態(tài)性給出了參數(shù)的近似置信區(qū)間。文獻(xiàn)[9]討論了可靠度函數(shù)的區(qū)間估計(jì)方法。文獻(xiàn)[10]討論了BS(α,β)分布在截尾試驗(yàn)情形下的統(tǒng)計(jì)分析，給出了參數(shù)的擬最小二乘估計(jì)和極大似然估計(jì)，另外還給出了刻度參數(shù)β的區(qū)間估計(jì)。文獻(xiàn)[11]研究了缺失數(shù)據(jù)場(chǎng)合下兩參數(shù)疲勞壽命分布BS(α,β)刻度參數(shù)β的區(qū)間估計(jì)。文獻(xiàn)[12]提出了一個(gè)新的樞軸量來(lái)構(gòu)造刻度參數(shù)β的經(jīng)典置信區(qū)間，并且具有較為簡(jiǎn)單的顯式表達(dá)式。文獻(xiàn)[13]在失效機(jī)理保持不變的條件下，討論兩參數(shù)BS疲勞壽命分布環(huán)境因子的區(qū)間估計(jì)問(wèn)題。文獻(xiàn)[14]研究了形狀參數(shù)以及均值、分位數(shù)、可靠度等可靠性指標(biāo)的廣義區(qū)間估計(jì)。文獻(xiàn)[15]通過(guò)Monte Carlo模擬說(shuō)明文獻(xiàn)[11-12]所給出的兩種方法可能無(wú)法得到兩參數(shù)疲勞壽命分布BS(α,β)刻度參數(shù)β區(qū)間估計(jì)，同時(shí)指出文獻(xiàn)[14]在利用廣義樞軸量法給出刻度參數(shù)β以及參數(shù)函數(shù)的置信區(qū)間過(guò)程中存在錯(cuò)誤，并用反例進(jìn)行了說(shuō)明，同時(shí)給出了正確的證明。文獻(xiàn)[16]應(yīng)用回歸分析方法給出了兩參數(shù)疲勞壽命分布BS(α,β)參數(shù)的最小二乘估計(jì)和形狀參數(shù)α的簡(jiǎn)單易算的區(qū)間估計(jì)方法，并用計(jì)算機(jī)隨機(jī)模擬方法研究了區(qū)間估計(jì)的效果，模擬結(jié)果顯示效果非常好。文獻(xiàn)[17]針對(duì)兩參數(shù)疲勞壽命分布BS(α,β)利用形狀參數(shù)α的區(qū)間估計(jì)給出了分布變異系數(shù)的區(qū)間估計(jì)和假設(shè)檢驗(yàn)方法。文獻(xiàn)[18]針對(duì)逐步增加的II型截尾數(shù)據(jù)給出了參數(shù)的近似極大似然估計(jì)。文獻(xiàn)[19]通過(guò)對(duì)數(shù)變換給出了求兩參數(shù)疲勞壽命分布BS(α,β)在全樣本場(chǎng)合下參數(shù)的對(duì)數(shù)矩估計(jì)，并通過(guò)大量Monte Carlo模擬比較了各種點(diǎn)估計(jì)方法的精度，基于對(duì)數(shù)變換通過(guò)一階泰勒展開(kāi)，將兩參數(shù)疲勞壽命分布BS(α,β)近似看作兩參數(shù)對(duì)數(shù)正態(tài)分布，由此得到了刻度參數(shù)β和形狀參數(shù)α的近似區(qū)間估計(jì)，且區(qū)間估計(jì)都存在，通過(guò)Monte Carlo模擬比較發(fā)現(xiàn)所給出的近似方法比原有方法更精確，最后通過(guò)若干實(shí)例說(shuō)明方法的可行性。

本文針對(duì)兩參數(shù)Birnbaum-Saunders疲勞壽命分布產(chǎn)品在雙邊定時(shí)截尾壽命試驗(yàn)數(shù)據(jù)場(chǎng)合下，利用“半個(gè)產(chǎn)品失效”的想法，并利用Balakrishan提出的一階泰勒展開(kāi)，將似然方程作適當(dāng)?shù)慕疲玫搅藘蓚€(gè)參數(shù)的近似極大似然估計(jì)，最后通過(guò)一實(shí)際案例說(shuō)明方法的可行性。

1 雙邊定時(shí)截尾壽命試驗(yàn)場(chǎng)合下參數(shù)的點(diǎn)估計(jì)

設(shè)產(chǎn)品的壽命T服從兩參數(shù)BS(α,β)分布，如果做雙邊定時(shí)截尾壽命試驗(yàn)，定時(shí)截尾時(shí)間設(shè)為τ1,τ2其中τ1<τ2，在τ1前共有r1-1個(gè)產(chǎn)品失效，在τ2前共有r2個(gè)產(chǎn)品失效，次序失效時(shí)間為：τ1≤t(r1)≤t(r1+1)≤…≤t(r2)≤τ2，也就是說(shuō)在時(shí)間區(qū)間[τ1,τ2]內(nèi)共有k=r2-r1+1個(gè)產(chǎn)品失效。此時(shí)，似然函數(shù)為(其中C+為正常數(shù))：

L(α,β)=C+[F(τ1;α,β)]r1-1[1-

lnL(α,β)=lnC++(r1-1)ln[F(τ1;α,β)]+

(3)

(4)

式(3)與式(4)為含參數(shù)α,β的二元超越方程組，要從中求解得參數(shù)α,β的極大似然估計(jì)是一個(gè)非常復(fù)雜的問(wèn)題。為解決這一問(wèn)題，下面通過(guò)近似處理，通過(guò)求解僅含參數(shù)β的一元超越方程得到β的近似極大似然估計(jì)，進(jìn)而得到參數(shù)α的近似極大似然估計(jì)。

似然函數(shù)為(其中C+為正常數(shù))：

(5)

lnL(α,β)=lnC++(r1-1)ln[Φ(Z1)]+

(n-r2)ln[1-Φ(Z2)]+

lnL(α,β)=lnC+-klnα-klnβ+(r1-1)ln[Φ(Z1)]+

(n-r2)ln[1-Φ(Z2)]+

(6)

(7)

(8)

(9)

(10)

觀察雙邊定時(shí)截尾壽命試驗(yàn)的數(shù)據(jù)形式，在時(shí)刻τ1前有r1-1個(gè)產(chǎn)品失效，而在時(shí)間t(r1-1)與τ1之間沒(méi)有產(chǎn)品失效，在此將其視作有“半個(gè)產(chǎn)品失效”。在時(shí)刻τ2前有r2個(gè)產(chǎn)品失效，而在時(shí)間t(r2)與τ2之間沒(méi)有產(chǎn)品失效，在此將其視作有“半個(gè)產(chǎn)品失效”。

利用上述泰勒展開(kāi)對(duì)式(9)化簡(jiǎn)得：

(r1-1)Z1(a(1)-b(1)Z1)-k=0

即

也即

(11)

其中，

利用上述泰勒展開(kāi)對(duì)式(10)化簡(jiǎn)得：

即

也即

(12)

(13)

(14)

特別地：當(dāng)r1=1時(shí)，為一般定時(shí)截尾壽命試驗(yàn)場(chǎng)合，即產(chǎn)品試驗(yàn)做到時(shí)間τ0為止，在τ0前有r個(gè)產(chǎn)品失效，此時(shí)τ1=0,τ2=τ0,r2=r。類似在可以給出一般定時(shí)截尾壽命試驗(yàn)參數(shù)的近似極大似然估計(jì)。

2 實(shí)例

本實(shí)例[21]中的數(shù)據(jù)集為6061-T6鋁合金的疲勞壽命。鋁合金的切取位置應(yīng)平行于軋制方向，振蕩頻率為18赫茲。該數(shù)據(jù)集包括101個(gè)觀測(cè)值，最大應(yīng)力為31 000帕。數(shù)據(jù)如下：70， 90，96， 97， 99，100，103，104，104，105，107，108，108，108，109，109，112，112，113，114，114，114，116，119，120，120，120，121，121，123，124，124，124，124，124，128，128，129，129，130，130，130，131，131，131，131，131，132，132，132，133，134，134，134，134，134，136，136，137，138，138，138，139，139，141，141，142，142，142，142，142，142，144，144，145，146，148，148，149，151，151，152，155，156，157，157，157，157，158，159，162，163，163，164，166，166，168，170，174，196，212

文獻(xiàn)[19]在全樣本場(chǎng)合下給出了14種點(diǎn)估計(jì)，如表1所示。

表1 全樣本場(chǎng)合下參數(shù)的點(diǎn)估計(jì)

從雙邊定時(shí)截尾壽命試驗(yàn)與一般定時(shí)截尾壽命試驗(yàn)的計(jì)算結(jié)果看與表1中的結(jié)果近似。