Weibull定時截尾情形下的矩估計性能

馮自立，陳曉陽，顧家銘

（1.上海大學 軸承研究室，上海 200072；2.上海天安軸承有限公司，上海 200233）

試驗是掌握產(chǎn)品性能最具說服力的手段，其關(guān)鍵在于建立和確認產(chǎn)品性能的概率分布。普遍認為投入試驗的試件越多，估計結(jié)果越精確，但由于受時間和資金的限制，產(chǎn)品的可靠性試驗不僅不能保證大的樣本量，還必須采取截尾試驗方案。在實踐中定時截尾試驗多于定數(shù)截尾試驗。Weibull分布常被用于描述滾動軸承的疲勞壽命，然而對Weibull分布參數(shù)估計的研究，大多集中于完全樣本試驗［1-2］和定數(shù)截尾試驗。文中通過MonteCarlo法研究矩估計量在小樣本量（樣本量不超過10）定時截尾情形下，對Weibull分布參數(shù)的估計性能。

1 小樣本量Weibull定時截尾數(shù)據(jù)的特性

通過MonteCarlo方法［3］研究Weibull分布下的小樣本量定時截尾數(shù)據(jù)，發(fā)現(xiàn)有以下規(guī)律［4］：（1）在一定的樣本容量和形狀參數(shù)β下，使得某一失效數(shù)出現(xiàn)概率最大的截尾時間，可以是某一區(qū)間上的任意取值；（2）對于同一樣本結(jié)構(gòu)（n，r）（n為樣本容量，r為失效數(shù)），不同Weibull分布（形狀參數(shù)β相同，尺度參數(shù)η不同）下的截尾時間區(qū)間端點與尺度參數(shù)的比值相同。將形狀參數(shù)等于1時，使得n＝10的各樣本結(jié)構(gòu)出現(xiàn)概率最大的截尾時間范圍列于表1。

表1 樣本結(jié)構(gòu)與對應的截尾時間范圍

2 常用估計方法

通?？刹捎肳eibull概率紙圖估計法、最小二乘法和極大似然法對Weibull分布下的定時截尾數(shù)據(jù)進行估計。最小二乘法改善了圖估計法易受人為主觀因素影響，精度不高的缺點，在實踐中應用廣泛。但是最小二乘法有一個基本假設(shè)［5］，即每個樣本點都包含了關(guān)于未知母體的相同信息量。然而從Weibull分布母體得到的順序統(tǒng)計量中，前幾個隨機變量的離散程度大于后續(xù)的隨機變量。通過對樣本點合理加權(quán)可以改善對所有試驗數(shù)據(jù)“一視同仁”的不足。文獻［5］通過引入一個中間隨機變量，成功解決了計算樣本點加權(quán)系數(shù)的問題。

對于截尾情形下壽命服從Weibull分布的試件，可在各個截尾樣本上用1-F（t）＝exp［-（t／η）β］代替概率密度函數(shù)與失效樣本所對應的概率密度函數(shù)值相乘，求出似然函數(shù)。

根據(jù)極大似然原理可以得到關(guān)于βMLE（形狀參數(shù)β的極大似然估計）的超越方程：

式中：ti（i＝1，2，…，r）為失效時間。方程的左端是βMLE的嚴格減連續(xù)函數(shù)，可采用二分法進行迭代求解。在高截尾的小樣本量定時截尾數(shù)據(jù)情形下，（1）式的計算不易收斂。

以從W（1，1）（β＝1，η＝1的Weibull分布母體）中抽取的容量為10的定時截尾樣本為例，形狀參數(shù)的極大似然估計量在各個樣本結(jié)構(gòu)下（均保證10 000次有效模擬計算）的5%，50%和95%分位數(shù)與截尾時間（表1）的關(guān)系如圖1所示。

圖1 形狀參數(shù)極大似然估計量的分位數(shù)與截尾時間（均為無量綱）的關(guān)系

由圖1可知：在高、中截尾場合（r＝2～6），極大似然估計量的分位數(shù)基本不隨截尾時間的變化而改變；在低截尾場合（r＝7～9），極大似然估計量的50%和95%分位數(shù)隨著截尾時間的增加明顯減小。這一規(guī)律同樣出現(xiàn)在其他樣本容量所對應的低截尾場合。最小二乘估計量和加權(quán)最小二乘估計量在截尾時間區(qū)間上同樣有這一規(guī)律。為便于比較估計量的性能，對于同一失效數(shù)，均取使其出現(xiàn)概率最大的截尾時間區(qū)間的左端點所對應的估計量分布進行分析。

以對來自于W（1，1），容量為10的定時截尾樣本估計為例，用箱線圖表示形狀參數(shù)的最小二乘估計量、加權(quán)最小二乘估計量以及極大似然估計量在各失效數(shù)下的90%分布區(qū)間（在各樣本結(jié)構(gòu)下對3種估計量均保證10 000次有效模擬計算），見圖2。

圖2 估計量的90%分布區(qū)間箱線圖

由圖2可知：（1）隨著失效數(shù)的增加，各估計量的精度均在提高，最小二乘估計量與加權(quán)最小二乘估計量逐漸趨于無偏（中位數(shù)意義上），而極大似然估計量逐漸變大；（2）極大似然估計法精度最高，加權(quán)最小二乘法次之，最小二乘法精度最低。

3 矩估計

3.1 矩估計原理

極值分布屬于位置-刻度分布族，其分布函數(shù)為：

若A～W（β，η）（隨機變量A服從形狀參數(shù)為β，尺度參數(shù)為η的兩參數(shù)Weibull分布），則ln A服從極值分布，其中參數(shù)μ＝lnη，σ＝1／β。

考慮位置-刻度分布族的好處就在于對定時截尾情形下的數(shù)據(jù)可以采用矩法進行估計［6］。文獻［7］利用Weibull分布與極值分布之間的轉(zhuǎn)換關(guān)系首先研究了Weibull分布下定時截尾數(shù)據(jù)的矩估計問題，并證明了其構(gòu)造的形狀參數(shù)與尺度參數(shù)的矩估計量具有強相合性。

中間變量h（p）為［6-7］：

式中：p＝P（X≤T）（設(shè)隨機變量X服從Weibull分布）。

文獻［7］給出的形狀參數(shù)與尺度參數(shù)的矩估計量為：

式中：T為截尾時間；采用pE（pE＝r／n）作為p的估計，計算出h（pE）作為h（p）的近似值。

為方便使用，給出樣本容量10以內(nèi)所有樣本結(jié)構(gòu)對應的h（pE），見表2。

表2 小樣本量下的h（p E）值

3.2 矩估計性能

在低截尾場合下，矩估計量的5%，50%和95%分位數(shù)同樣隨著截尾時間的增加而減小。以對母體為W（2，1）的截尾樣本（10，9）估計為例，其矩估計量的分位數(shù)隨截尾時間變化的曲線如圖3所示。

對于不同的失效數(shù)，取其對應的截尾時間區(qū)間左端點上的極大似然估計量和矩估計量的分布進行比較。

對來自于W（2，1），容量為10的定時截尾樣本進行估計（在每個樣本結(jié)構(gòu)下均保證10 000次有效模擬計算），形狀參數(shù)的矩估計量與極大似然估計量的90%分布區(qū)間如圖4所示。

圖3 低截尾場合下形狀參數(shù)矩估計量的分位數(shù)隨截尾時間（無量綱）變化的曲線

圖4 估計量的90%分布區(qū)間箱線圖

由圖4可知：（1）在小樣本量定時截尾情形下，矩估計隨著失效數(shù)的增加逐漸變大；（2）矩估計與極大似然估計的精度相當。

4 實例計算

以某型陀螺電動機轉(zhuǎn)子軸承小樣本精度壽命模擬試驗為例說明上述方法的可靠性估計［8］，試驗采用定時截尾方案（截尾時間為4 000 h）。共投入8組（每組2套）軸承，其中5組失效，壽命時間依次為：1 313，2 288，2 472，2 506和3 382 h，其余3組截尾。分布擬合優(yōu)度檢驗結(jié)果表明，可以接受該型軸承精度壽命服從Weibull分布的假設(shè)。分別采用極大似然估計法（1）式和矩估計法（3）式對壽命數(shù)據(jù)進行估計，結(jié)果見表3。

表3 陀螺轉(zhuǎn)子軸承在3種置信水平下的可靠性壽命估計

此外，陀螺生產(chǎn)廠家用10套該型軸承裝配5臺某型高速陀螺電動機，進行了通電壽命試驗。每臺陀螺電動機累計通電工作時間達1 023 h，均滿足性能指標。

陀螺電動機通電壽命試驗結(jié)果表明：兩種估計方法的計算結(jié)果相近，均比較符合實際；矩估計較極大似然估計略偏保守。

5 結(jié)論

（1）極大似然估計較最小二乘估計以及加權(quán)最小二乘估計更加精確。

（2）矩估計無需迭代，在高截尾情形下不會出現(xiàn)極大似然估計難以收斂的情況。

（3）矩估計與極大似然估計精度相當。