兩參數(shù)Birnbaum-Saunders疲勞壽命分布的統(tǒng)計分析方法研究

王蓉華，徐曉嶺，顧蓓青

（1.上海師范大學(xué) 數(shù)理學(xué)院，上海200234;2.上海對外經(jīng)貿(mào)大學(xué) 商務(wù)信息學(xué)院，上海201620）

兩參數(shù)Birnbaum-Saunders疲勞壽命分布的統(tǒng)計分析方法研究

王蓉華1，徐曉嶺2，顧蓓青2

（1.上海師范大學(xué) 數(shù)理學(xué)院，上海200234;2.上海對外經(jīng)貿(mào)大學(xué) 商務(wù)信息學(xué)院，上海201620）

文章首先比較了全樣本場合下兩參數(shù)BS疲勞壽命分布各種參數(shù)的點估計的優(yōu)劣，認(rèn)為當(dāng)樣本容量較大時，用分位數(shù)估計更為方便。通過Monte Carlo模擬說明文獻[6]關(guān)于尺度參數(shù)的區(qū)間估計有可能是不存在的。

兩參數(shù)Birnbaum-Saunders疲勞壽命分布；點估計；區(qū)間估計；尺度參數(shù)

0　引言

Birnbaum-Saunders模型[1]是概率物理方法中的一個重要失效分布模型，這個模型是Birnbaum和Sauders在1969年研究主因裂紋擴展導(dǎo)致的材料失效過程中推導(dǎo)出來的，這一模型在機械產(chǎn)品可靠性研究中應(yīng)用廣泛，主要應(yīng)用于疲勞失效研究，對于電子產(chǎn)品性能退化失效分析也有重要應(yīng)用。

設(shè)X服從兩參數(shù)Birnbaum-Saunders疲勞壽命分布BS(α，β)，其分布函數(shù)F(x)與密度函數(shù)分別為：

其中，α＞0稱為形狀參數(shù)，β＞0稱為刻度參數(shù)（或者稱為尺度參數(shù)），φ(x)，Φ(x)分別為標(biāo)準(zhǔn)正態(tài)分布的密度函數(shù)與分布函數(shù)，即

由于Birnbaum-Saunders疲勞壽命分布是從疲勞過程的基本特征出發(fā)導(dǎo)出的，因此它比常用壽命分布如威布爾分布和對數(shù)正態(tài)分布更適合描述某些由于疲勞而引起失效的產(chǎn)品壽命規(guī)律。此分布已經(jīng)成為可靠性統(tǒng)計分析中的常用分布之一。

1　參數(shù)的幾種點估計方法的綜合比較

1.1參數(shù)的極大似然估計（方法一）

設(shè)X1，X2，…，Xn為來自Birnbaum-Saunders疲勞壽命分布總體X～BS(α，β)的一個容量為n的簡單隨機樣本，其次序觀察值為x1，x2，…，xn。

似然函數(shù)為：

化簡得僅含參數(shù)β的超越方程：

極大似然估計需要解非線性方程組,且計算較為復(fù)雜。由文獻[2]可知，當(dāng)α＞2時，似然函數(shù)有多個極大值點。

1.2矩估計（方法二）

表1　100000次模擬中滿足的頻率

表1　100000次模擬中滿足的頻率

n 20 20 20 30 30 30 40 40 40 α真值0.25 0.75 1 0.25 0.75 1 0.25 0.75 1 β真值0.5 1 1 0.5 1 1 0.5 1 1頻率1.0000 1.0000 0.9999 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

由此可以看出當(dāng)α，β取值較小時矩估計1是可行的。

1.3矩估計（方法三）

由于X-1～BS(α，β-1)，于是可列矩方程組：

1.4逆矩估計（方法四）

由于獨立同分布，服從于N(0，1)，由逆矩估計思想建立如下方程：

化簡可得參數(shù)β，α的逆矩估計分別為：

1.5分位數(shù)估計(方法五)

1.6回歸估計方法(方法六)

孫祝嶺在文獻[3]中利用回歸分析模型，給出了參數(shù)的最小二乘估計如下：

1.7六種點估計的模擬比較

取樣本容量n=15，20，30，40，給定參數(shù)α，β真值，通過10000次模擬，分別計算其六種點估計的均值和均方誤差，模擬結(jié)果如表2所示：

表2　六種估計方法的模擬比較

續(xù)表2

從模擬結(jié)果來看，α，β真值較小的時候效果較好，并且這幾種方法的效果也相近，方法二的精度稍差些。當(dāng)樣本容量n較大時，分位數(shù)估計的效果要好些，尤其是對參數(shù)α的估計較為精確，同時分位數(shù)估計計算也較為方便，適合工程應(yīng)用。

例1[3]：設(shè)總體X～BS(α，β)，來自X的為10的樣本觀察值為：1401.7,2170.3,3762.9,414.0,650.7,441.0,902.8, 806.6,1292.1,760.3，計算α，β的估計值。

例2[4]：數(shù)據(jù)來自Bjerkedal(1960),也被Guptaetal(1997)分析過。數(shù)據(jù)表示幾內(nèi)亞豬注射不同劑量的結(jié)核桿菌的生存時間。對于文中的養(yǎng)殖方法，有72個觀測值如下：

12，15，22，24，24，32，32，33，34，38，38，43，44，48，52，53，54，54，55，56，57，58，58，59，60，60，60，60，61，62，63，65，65，67，68，70，70，72，73，75，76，76，81，83，84，85，87，91，95，96，98，99，109，110，121，127129，131，143，146，146，175，175，211，233，258，258，263，297，341，341，376.

例3[5]：碳纖維的破壞應(yīng)力（GPa）如下所示：

2　文獻[6]關(guān)于參數(shù)β的區(qū)間估計的討論

設(shè)X1，X2，…，Xn為來自Birnbaum-Saunders疲勞壽命分布總體X～BS(α，β)的一個容量為n的簡單隨機樣本，其次序觀察值為x1，x2，…，xn。

又：

則針對函數(shù) g2(β)，存在 β0＞0，且當(dāng)0＜β＜β0時，g2(β)＜0；當(dāng)β＞β0時，g2(β)＞0

記：

也即，當(dāng)0＜β＜β0時嚴(yán)格單調(diào)下降；當(dāng)β＞β0時嚴(yán)格單調(diào)上升。

表3　1000次Monte Carlo模擬比較

[1]Birnbaum Z W,Saunders S C.A New Family of Life Distribution[J]. Appl.Prob,1969,(6).

[2]王炳興,王玲玲.Birnbaum-Saunders疲勞壽命分布的參數(shù)估計[J].華東師范大學(xué)學(xué)報,1996,(4).

[3]孫祝嶺.Birnbaum-Saunders疲勞壽命分布參數(shù)的回歸估計方法[J].兵工學(xué)報，2010，31(9).

[4]Kundu D,Kannan N,Balakrishnan N.On the Hazard Function of Birnbaum-Saunders Distribution and Associated Inference[J].Compu?tational Statistics&Data Analysis,2008,(52).

[5]Gauss M,Cordeiro A.Lemonte J.TheβBirnbaum-Saunders Distribu?tion:An Improved Distribution for Fatigue Life Modeling[J].Computa?tional Statistics and Data Analysis.2011,(55).

[6]孫祝嶺.Birnbaum-Saunders疲勞壽命分布尺度參數(shù)的區(qū)間估計[J].兵工學(xué)報，2009，30（11）.

（責(zé)任編輯/易永生）

O212

A

1002-6487（2016）22-0027-04

國家自然科學(xué)基金資助項目（11671264）；全國統(tǒng)計科學(xué)研究計劃重點項目（2013LZ08）；上海市教育委員會科研創(chuàng)新一般項目（14YZ080）

王蓉華（1972—），女，上海人，博士，副教授，研究方向：應(yīng)用統(tǒng)計。