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

何 亮，徐曉嶺

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

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

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

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

1 基本假設(shè)

統(tǒng)計(jì)分析基于以下幾個假設(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)品的失效機(jī)制相同，即形狀參數(shù)α保持不變。

假設(shè)3 特征壽命β(V)與應(yīng)力V滿足逆冪律模型：β(V)=1/(dVc)，其中d>0,c>0為待估參數(shù)，且獨(dú)立于應(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é)表達(dá)為：在應(yīng)力水平Vi下產(chǎn)品連續(xù)工作ti時間的累積失效概率Fi(ti)相當(dāng)于該產(chǎn)品在應(yīng)力水平Vj下工作tj時間的累積失效概率，即Fi(ti)=Fj(tj)。

在序加試驗(yàn)中，將n個產(chǎn)品置于序進(jìn)應(yīng)力V(t)=Kt+V0(V0≥0)下進(jìn)行加速壽命試驗(yàn)，直至所有受測樣本失效為止。關(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 逆冪律模型下序進(jìn)應(yīng)力加速壽命試驗(yàn)的統(tǒng)計(jì)分析

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

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

(2)

(3)

試驗(yàn)得到的產(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估計(jì)(方法2)

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

在序進(jìn)應(yīng)力V(t)=Kt+V0下，產(chǎn)品壽命次序統(tǒng)計(jì)量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 最小距離估計(jì)(方法3)

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

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

(5)

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

(6)

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

3 數(shù)值模擬

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

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

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

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

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

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

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

模擬算例：取真值為α=1,d=1,c=1，在序進(jìn)應(yīng)力V(t)=t+1下進(jìn)行序加試驗(yàn)，利用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