何 亮，徐曉嶺

(上海對外經(jīng)貿(mào)大學(xué) 統(tǒng)計與信息學(xué)院，上海 201620)

隨著科技的不斷發(fā)展，產(chǎn)品更新?lián)Q代速度加快，電子設(shè)備及元器件的可靠性在不斷提高，一般的產(chǎn)品壽命試驗方法已經(jīng)不再適用。所以通常設(shè)計加速壽命試驗對高可靠性、長壽命的產(chǎn)品的壽命進行快速評定。序進應(yīng)力加速壽命試驗(簡稱序加試驗)是加速壽命試驗中的一種重要方法，試驗中受試產(chǎn)品的應(yīng)力水平隨時間連續(xù)上升，是時間的連續(xù)非降函數(shù)，相比于恒定應(yīng)力加速壽命試驗?zāi)芸s短試驗時間。

由Frechet在1927年提出的Frechet分布[1]在各種工程應(yīng)用的統(tǒng)計建模中非常重要，在極端氣候事件中也有相當(dāng)廣泛的應(yīng)用。Mann在文獻[2]中討論了Frechet分布及Gumbel分布間關(guān)系與其參數(shù)估計。Longin在文獻[3]中將Frechet分布運用到擬合股市的極端價格變動下的走勢。Nadarajah在文獻[4]將Frechet分布用于社會學(xué)模型。Zaharim等在文獻[5]中將對數(shù)正態(tài)分布以及Frechet分布用于風(fēng)速數(shù)據(jù)的預(yù)測分析中。Abd-Elfattah在不同的標(biāo)準(zhǔn)下對廣義Frechet分布進行了擬合優(yōu)度檢驗[6]。Mubarak在文獻[7]中給出了逐步增加定數(shù)截尾場合下(移除數(shù)量服從二項分布)的參數(shù)極大似然估計以及區(qū)間估計。Abbas和湯銀才在文獻[8]中對已知形狀參數(shù)的Frechet分布進行Bayes估計，在[9]中研究了在定數(shù)截尾數(shù)據(jù)場合下參數(shù)的極大似然估計以及最小二乘估計。國外的很多學(xué)者對一般化的Frechet分布有較多的研究，Nadarajah與Kotz在文獻[10]中提出指數(shù)化的Frechet分布，與Gupta提出Beta-Frechet分布[11]。Krishna在文獻[12]介紹了Marshall-Olkin Frechet分布。Afify提出Weibull-Frechet 四參數(shù)壽命分布模型[13]，并通過極大似然估計其模型參數(shù)。Mansoor等引入擴展的四參數(shù)Frechet分布，亦稱指數(shù)化的指數(shù)Frechet分布，并與Marshall-Olkin Frechet分布、指數(shù)Frechet分布、Frechet分布基于實際數(shù)據(jù)作比較，說明了新分布更加貼合這些數(shù)據(jù)[14]。在加速壽命試驗方面，Shahab在文獻[15]中基于指數(shù)過程在常應(yīng)力下給出分布的參數(shù)估計以及區(qū)間估計，在[16]中通過圖像研究Frechet 分布簡單步加試驗的最優(yōu)時間。Ghaly在文獻[17]中在常應(yīng)力下對指數(shù)分布族加速壽命試驗的不同方法進行了對比，并基于實際數(shù)據(jù)給出模型選用的建議。

本文基于逆冪律模型對Frechet分布在與時間為一般線性關(guān)系的序進應(yīng)力下的加速壽命試驗進行統(tǒng)計分析，分別運用極大似然估計、MPS估計以及最小距離估計分別給出參數(shù)的點估計，并通過大量的Monte-Carlo模擬對比3種估計方法的精度。

1 基本假設(shè)

統(tǒng)計分析基于以下幾個假設(shè)：

假設(shè)1 在恒定應(yīng)力V>0下，產(chǎn)品壽命X服從Frechet分布，記為X～Fr(α,β(V))，其分布函數(shù)為：FV(t)=exp{-(β(V)/t)α}，其中α>0為形狀參數(shù)，β(V)>0為尺度參數(shù)。

假設(shè)2 不同應(yīng)力水平下，產(chǎn)品的失效機制相同，即形狀參數(shù)α保持不變。

假設(shè)3 特征壽命β(V)與應(yīng)力V滿足逆冪律模型：β(V)=1/(dVc)，其中d>0,c>0為待估參數(shù)，且獨立于應(yīng)力V,φ(V)為已知函數(shù)。對模型兩邊取對數(shù)，可得到β(V)滿足對數(shù)線性關(guān)系：lnβ(V)=a+bφ(V)，其中a=-lnd,b=-c,φ(V)=lnV

假設(shè)4 產(chǎn)品的剩余壽命僅依賴于當(dāng)時已累積失效部分及當(dāng)時的應(yīng)力水平，與累積方式無關(guān)，此即Nelson假定。數(shù)學(xué)表達為：在應(yīng)力水平Vi下產(chǎn)品連續(xù)工作ti時間的累積失效概率Fi(ti)相當(dāng)于該產(chǎn)品在應(yīng)力水平Vj下工作tj時間的累積失效概率，即Fi(ti)=Fj(tj)。

在序加試驗中，將n個產(chǎn)品置于序進應(yīng)力V(t)=Kt+V0(V0≥0)下進行加速壽命試驗，直至所有受測樣本失效為止。關(guān)于在此應(yīng)力下的產(chǎn)品壽命分布，有如下定理。

定理1 若產(chǎn)品壽命滿足上述基本假設(shè)，當(dāng)應(yīng)力為關(guān)于時間的函數(shù)，即產(chǎn)品壽命在應(yīng)力V=V(t)下的分布函數(shù)為：

(1)

當(dāng)ξ→0時，由Riemann積分的定義并由假設(shè)3有β[V(t)]=exp{a+bφ[V(t)]}，則在應(yīng)力V(t)下產(chǎn)品壽命的分布函數(shù)為：

2 逆冪律模型下序進應(yīng)力加速壽命試驗的統(tǒng)計分析

2.1 極大似然估計(方法1)

根據(jù)定理1，產(chǎn)品在應(yīng)力V(t)=Kt+V0(V0>0)下壽命的分布函數(shù)和密度函數(shù)為

(2)

(3)

試驗得到的產(chǎn)品失效數(shù)據(jù)為t(1)≤t(2)≤…≤t(n)。

其似然函數(shù)為

對數(shù)似然函數(shù)為

lnL=nlnα-nαlnd+n(α+1)lnK+

分別令對數(shù)似然函數(shù)對α,d,c求導(dǎo)：

2.2 MPS估計(方法2)

在極值分布的參數(shù)估計中，小樣本場合下的極大似然估計穩(wěn)定性較差，而相比之下由Cheng和Amin提出的MPS估計具有更穩(wěn)定的估計結(jié)果[20-21],甚至能在某些極大似然估計失效的情況下仍能作出良好的評估[22]。

在序進應(yīng)力V(t)=Kt+V0下，產(chǎn)品壽命次序統(tǒng)計量T(i)的分布函數(shù)記為FT(i)(t)，取值為t(i)，則有

FT(i)(t)=IF*(t)(i,n-i+1),i=1,2,…,n

(4)

其中當(dāng)i=1,2,…,n時，

[1-(1+c)ln(Kt(i)+V0)]}·

2.3 最小距離估計(方法3)

Wolfowitz提出了最小距離估計(或稱擬合優(yōu)度估計)[24-25]，并證明了估計的一致性。其思想是通過使經(jīng)驗分布與某些參數(shù)族分布的距離達到最小。常用的距離準(zhǔn)則有卡方距離、Kolmogorov-Smirnov (K-S) 距離、Anderson-Darling (A-D) 距離以及Cramer-Von Mises (CVM) 距離。本文選取Anderson-Darling (A-D) 距離作為擬合優(yōu)度統(tǒng)計量。

經(jīng)驗分布函數(shù)Fn(t)與分布函數(shù)F(t)的距離記為EDF統(tǒng)計量，即：

(5)

其中ψ(t)為權(quán)重函數(shù)。令ψ(t)={F(t)[1-F(t)]}-1，為方便計算，A-D統(tǒng)計量的代數(shù)形式如下：

(6)

其中z(i)=F(t(i)),i=1,2,…,n,t(i)為第i個次序統(tǒng)計量的值，令

3 數(shù)值模擬

針對真值為α=0.5,1,1.5,d=1,c=1,樣本容量為n=20,30,50,在序進應(yīng)力V(t)=t+1下分別進行 50 000次的Monte-Carlo模擬試驗，試驗結(jié)果如表1所示。

表1 序進應(yīng)力為V(t)=t+1下Monte-Carlo模擬試驗結(jié)果

從表1中的結(jié)果，可以得到如下結(jié)論：

1) 關(guān)于參數(shù)α，就均方誤差意義上而言方法2的效果最好。當(dāng)α較小(<1)時，方法3誤差大于方法1，但估計均值更加接近參數(shù)真值；當(dāng)α≥1時，方法3的誤差介于方法1和方法2之間，且估計均值較方法1、方法2更接近真實情況。

2) 關(guān)于參數(shù)d，方法2的估計效果最好。方法1的誤差在3種估計方法中最大，而當(dāng)α=0.5時其估計均值離真值最遠，當(dāng)α=1,1.5時較方法3而言更接近真值。

3) 關(guān)于參數(shù)c，方法1的估計均值離真值較遠，也最不穩(wěn)定；方法2效果次之，方法3效果最佳，甚至在小樣本情況下仍能得到較好的估計結(jié)果。

在實際應(yīng)用中，應(yīng)通過方法2(MPS估計)得到參數(shù)α,d的估計值，通過方法3(最小距離估計)得到參數(shù)c的估計值，以達到最佳的估計效果。下面給出模擬算例。

模擬算例：取真值為α=1,d=1,c=1，在序進應(yīng)力V(t)=t+1下進行序加試驗，利用Monte-Carlo模擬生成30個產(chǎn)品失效數(shù)據(jù)如下：

0.238 50.293 40.325 40.361 70.385 00.445 70.470 40.486 30.546 70.610 30.668 70.772 20.803 70.917 50.941 80.970 50.993 51.009 61.048 31.074 61.143 91.514 11.519 91.659 52.27162.892 02.958 73.098 73.90 749.294 2