兩參數拉普拉斯BS 疲勞壽命分布的統計分析

徐曉嶺，顧蓓青*，王蓉華

（1.上海對外經貿大學 統計與信息學院，上海 201620；2.上海師范大學 數理學院，上海 200234）

設T服從兩參數Birnbaum-Saunders 疲勞壽命分布BS(α，β)，其分布函數與密度函數分別為

其中，α＞0，為形狀參數，β＞0，為刻度參數（或尺度參數），φ(x)和Φ(x)分別為標準正態分布的密度函數和分布函數，即

BIRNBAUM 等［1］在研究主因裂紋擴展導致材料失效過程中得到了兩參數BS 疲勞壽命分布，此分布較威布爾分布、對數正態分布等常用的壽命分布更適合描述由疲勞引起失效的產品的壽命規律，已成為可靠性工程中的常用分布之一。值得指出的是，國內學者楊德滋等［2］在混凝土重復荷載作用下，通過求解Fokker-Planck 方程得到的疲勞壽命服從兩參數BS 疲勞壽命分布，分析了混凝土抗拉疲勞壽命試驗資料，并認為兩參數BS 疲勞壽命分布較對數正態分布更接近實際數據。關于兩參數BS 疲勞壽命分布的研究已有很多［3-19］。

隨著研究的深入，兩參數BS 疲勞壽命分布被進一步推廣為廣義BS 疲勞壽命分布。如文獻［20-23］中均涉及廣義BS 疲勞壽命分布，并將原有的兩參數BS 疲勞壽命分布所涉及的標準正態分布N(0，1)改為標準拉普拉斯分布L(1)，即兩參數拉普拉斯BS 疲勞壽命分布。文獻［23］較為詳細地研究了兩參數拉普拉斯BS 疲勞壽命分布，給出了其數字特征，研究了其密度函數、失效率函數的圖像特征及參數的極大似然估計，證明了極大似然估計是唯一存在的。下文將說明其證明過程中存在錯誤。

1 兩參數拉普拉斯BS 疲勞壽命分布的極大似然估計

設非負隨機變量T的分布函數

則稱T服從兩參數拉普拉斯BS 疲勞壽命分布，記為T～LBS(α，β)，其中α＞0，為形狀參數，β＞0，為刻度參數。

易見，其密度函數

另外，兩參數拉普拉斯BS 疲勞壽命分布LBS(α，β)的密度函數在t=β處的二階導數不存在。事實上，易見

綜上所述，文獻［23］引理5.1 中的β0不是唯一的，也可能不存在。

例 1給定樣本容量n，參數真值取α=1，β=1，通過Monte-Carlo 模擬產生3組容量為n的樣本數據，如表1 所示。

表1 模擬產生的3 組樣本數據Table 1 The simulation sample data of three groups

圖1 第1 組樣本數據的gj()圖像Fig.1 Image of gj()for the first sample data

圖2 第2 組樣本數據的gj()圖像Fig.2 Image of gj()for the second sample data

圖3 第3 組樣本數據的gj()圖像Fig.3 Image of gj()for the third sample data

于是，可認為LBS(α，β)分布參數的極大似然估計的唯一性有待進一步研究。

2 兩參數拉普拉斯BS 疲勞壽命分布的參數點估計

2.1 參數的矩估計（方法1）

設總體T～LBS(α，β)，T1，T2，…，Tn為來自總體T的容量為n的簡單隨機樣本。由文獻［23］知，

E(T)=β(α2+1)，D(T)=β2α2(11α2+2)。由矩估計思想，建立方程：

取樣本容量n=5，10，…，50，參數真 值α=0.25，0.50，…，2.00，β=1，通 過10 000 次Monte-Carlo 模擬，得 到滿足a≥11 的次數，見 表2，從 表2可以看出：（i）絕大多數樣本數據滿足a＜11；（ii）固定α，β，隨樣本容量n的增加，滿足a＜11 的次數呈減少趨勢；（iii）固定n，β，隨α的增加，滿足a＜11 的次數呈減少趨勢。

表2 10 000 次模擬中滿足a ≥11 的次數Table 2 The number that satisfies a ≥11 in 10 000 simulations

2.2 參數的矩估計（方法2）

2.3 參數的逆矩估計（方法3）

2.4 參數的分位數估計（方法4）

2.5 對數矩估計（方法5）

引理1[24]設X1，X2，…，Xn是來自總體X的容量為n的簡單隨機樣本，記E(X)=μ，D(X)=σ2＜+∞，該總體X的四階中心矩v4=E(X-EX)4有限，若函數h(x)的四階導數存在且有界，則有

雖然由引理2 知方程（2）有唯一正實根，但考慮方程（2）是含有積分且運算復雜的超越方程，不適合工程應用。在此仍采用極大似然估計的結果，即取相應α的估計為

2.6 參數點估計方法的模擬比較

取樣本容量n=5，10，…，20，參數真 值α=0.25，0.50，…，1.50，β=1，通過1 000 次Monte-Carlo 模擬，計算5 種點估計方法的均值及均方誤差，結果見表3，從表3 可以看出：（?。? 種方法的點估計精度均隨樣本容量的增加而增高；（ⅱ）綜合比較5種方法的均值與均方誤差，其中對數矩估計方法的優勢較明顯，特別是當參數α較大時。

3 兩參數拉普拉斯BS 疲勞壽命分布的參數的近似區間估計

3.1 參數β 的近似區間估計

3.1.1 基于文獻［25］中方法的參數β的近似區間估計

設總體T～LBS(α，β)，T1，T2，…，Tn為來自總體T的容量為n的簡單隨機樣本，其次序統計量記為T(1)≤T(2)≤…≤T(n)。

即G(β) 為樞軸量，其分布與參數無關。通過Monte-Carlo 模擬，計算得到G(β)的上側分位數，并由文獻［25］知，G(β)對β嚴格單 調遞減。對G(β)做恒等變換：

給定置信水平1-γ，記G(β) 分布的上側分位數 分別為ωγ/2，ω1-γ/2，若滿足，則參數β的置信水平為1-γ的區間估計為[β1，β2]，其中β1是G(β)=ωγ/2的根，β2是G(β)=ω1-γ/2的根。

給定顯著性水平γ=0.10，0.05，取樣本容量n=5，10，…，30，參數α=0.25，0.50，…，1.50，β=1，通過10 000 次Monte-Carlo 模擬，滿足ω1-γ/2＞的次數分別記為k1，k2和k，結果見表4，從表4 可以看出，k值很小，即利用樞軸量G(β)求參數β的區間估計在實際中并不可行，利用文獻［25］中的方法不大可能求得參數β的近似區間估計。值得一提的是，針對兩參數BS(α，β)分布，文獻［26］也指出相同的問題。

3.1.2 基于兩參數拉普拉斯分布L(μ，α)的近似方法的參數β的近似區間估計

表4 在10 000 次模擬中滿足的次數Table 4 The number that satisfiesin 10 000 simulations

表4 在10 000 次模擬中滿足的次數Table 4 The number that satisfiesin 10 000 simulations

即X=lnT，可將其近似地看作兩參數拉普拉斯分布L(μ，α)。

設總體T～LBS(α，β)，T1，T2，…，Tn為來自總體T的容量為n的簡單隨機樣本，其次序統計量記為T(1)≤T(2)≤…≤T(n)。令

則可將X1，X2，…，Xn近似地 看作來 自總體X～L(μ，α)的容量為n的簡單隨機樣本，其次序統計量記為X(1)≤X(2)≤…≤X(n)。令

則可將Y1，Y2，…，Yn近似地看作來自標準拉普拉斯總體Y～L(0，1)的容量為n的簡單隨機樣本，從小到大排序，記為Y(1)≤Y(2)≤…≤Y(n)。

于是，可將H(μ)看作僅含參數μ的樞軸量。通過Monte-Carlo 模擬，計算得到H(μ)的上側分位數。注意到H(μ)為μ的嚴格單調增函數，且

給定顯著性水平γ，樞軸量H(μ)的上側1-γ/2，γ/2 的分位數分別記為H1-γ/2和Hγ/2，易見參數μ的置信水平為1-γ的近似區間估計為

于是，可將H(μ)看作僅含參數μ的樞軸量。通過Monte-Carlo 模擬，計算得到H(μ)的上側分位數。注意到H(μ)為μ的嚴格單調增函數，且

給定顯著性水平γ，樞軸量H(μ)的上側1-γ/2，γ/2 分位數分別記為H1-γ/2和Hγ/2，易見參數μ的置信水平為1-γ的近似區間估計為

進而，參數β的置信水平為1-γ的近似區間估計為

給定置信水平0.90，取樣本容量n=10，15，…，30，參數真 值α=0.25，0.50，…，1.50，β=1，通過1 000 次Monte-Carlo 模擬，計算參數β的近似區間估計的平均下限、平均上限、平均長度及近似區間估計包含α的次數（k1），結果見表5，從表5 可以看出：（?。┕潭é粒瑓^間估計的長度隨n的增加而減小，即n越大，區間估計越精確；（ⅱ）k1均在900 以上。

表5 利用兩參數拉普拉斯分布求得到參數β，α 的近似區間估計Table 5 Approximate interval estimations of parameters β，α by two-parameter Laplace distribution

3.2 參數α 的近似區間估計

3.2.1 基于兩參數拉普拉斯分布L(μ，α)近似方法的參數α的近似區間估計

類似于3.1.2 節中參數β的近似區間估計方法，構造以下僅含參數α的樞軸量：

即可將T(α)看作僅含參數α的樞軸量。通過Monte-Carlo模擬，計算得到樞軸量T(α)的上側分位數。

給定顯著性水平γ，樞軸量T(α) 的上側1-γ/2，γ/2 分位數分別記為T1-γ/2和Tγ/2，易得參數α的置信水平為1-γ的近似區間估計為

給定置信水平0.90，取樣本容量n=10，15，…，30，參數真 值α=0.25，0.50，…，1.50，β=1，通 過1 000 次Monte-Carlo 模擬，計算α的近似區間估計的平均下限、平均上限、平均長度，以及近似區間估計包含α的次數（k2），結果見表5，從表5 可以看出：（?。┕潭é?，區間估計的長度隨n的增加而減小，即n愈大，區間估計愈精確；（ⅱ）當α較小時，k2均在900以上，而當α較大時（從模擬結果看，α＞1），k2未達到900，此時，置信水平未達到0.90。

3.2.2 基于參數的參數α的近似區間估計

表6 基于參數的參數α 的近似區間估計Table 6 Approximate interval estimation of parameter α based on parameter

表6 基于參數的參數α 的近似區間估計Table 6 Approximate interval estimation of parameter α based on parameter

給定置信水平0.90，取樣本容量n=10，15，…，30，參數真值α=0.25，0.50，…，1.50，β=1，通過1 000次Monte-Carlo 模擬，計算α的近似區間估計的平均下限、平均上限、平均長度，以及近似區間估計包含α的次數（k），結果見表6，從表6 可以看出：（?。┕潭é?，區間估計的長度隨n的增加而減小，即n越大，區間估計越精確；（ⅱ）k均在900 以上。

比較表5 和表6，發現基于參數的參數α的近似區間估計方法更優，不僅近似區間估計包含α的次數均達到900 以上，而且所得近似區間估計的平均長度較小。