苑延華

(黑龍江科技大學(xué) 理學(xué)院, 哈爾濱 150022)

?

廣義I型逐階區(qū)間刪失混合Weibull數(shù)據(jù)的參數(shù)估計

苑延華

(黑龍江科技大學(xué) 理學(xué)院, 哈爾濱 150022)

源于有限混合總體的廣義I型逐階區(qū)間刪失數(shù)據(jù)參數(shù)估計方法的研究不多，基于有限混合Weibull模型，討論Expectation-Maximization(EM)算法對廣義I型逐階區(qū)間刪失數(shù)據(jù)參數(shù)估計的有效性及其改進(jìn)。首先給出估計參數(shù)的EM算法，通過仿真算例，說明EM算法對廣義I型逐階區(qū)間刪失混合數(shù)據(jù)參數(shù)估計產(chǎn)生了過度迭代現(xiàn)象，進(jìn)而，提出了停止EM算法的加權(quán)絕對偏差信息準(zhǔn)則。改進(jìn)的EM算法改善了EM算法無法確定參數(shù)估計停止迭代時刻的不足，在選擇適當(dāng)初值后，可快速獲得滿意的參數(shù)估計結(jié)果。仿真算例驗證了該方法的有效性。

廣義I型逐階區(qū)間刪失; 加權(quán)絕對偏差信息; 有限混合Weibull分布; EM算法

0　引　言

I型刪失、II型刪失和逐階刪失是工業(yè)壽命檢驗和醫(yī)療生存分析中最常見的刪失方案,文中,考慮與I型逐階區(qū)間刪失[1-2](Progressive Type-I Interval Censoring)不同的廣義I型逐階區(qū)間刪失數(shù)據(jù)并假設(shè)個體壽命服從兩組分五參數(shù)混合Weibull分布。

在可靠性分析或壽命分析領(lǐng)域中,帶有逐階刪失的混合分布模型正吸引越來越多的研究者的關(guān)注與研究,取得了一些卓越的成果。文獻(xiàn)[1-9]主要從單一分布總體的角度考慮了總體分布參數(shù)的估計問題,常用參數(shù)估計的方法為最大似然估計、貝葉斯估計等方法。如E. B. Mokhtari等[3]用最大似然估計和漸近最大似然估計法討論了帶有II型逐階混合刪失的Weibull分布的參數(shù)估計。而文獻(xiàn)[10-15]主要討論了關(guān)于混合分布參數(shù)的最大似然、EM算法等估計方法。如C. T. Lin等[10]用最大似然和貝葉斯估計法研究了I型自適應(yīng)逐階混合刪失Weibull分布的參數(shù)估計。受這些研究的啟發(fā),筆者研究廣義I型逐階區(qū)間刪失混合Weibull分布下,參數(shù)最大似然估計EM算法的有效性問題,這個主題國內(nèi)外研究者鮮有觸及。

1　刪失數(shù)據(jù)的定義

圖1　I型逐階區(qū)間刪失數(shù)據(jù)的觀測方案

圖2　廣義I型逐階區(qū)間刪失數(shù)據(jù)的觀測方案

稱這樣的數(shù)據(jù)是廣義I型逐階區(qū)間刪失的,主要有兩個原因:

(2)補充數(shù)據(jù)Ri=di可以理解為iΔ處隨機刪失的個體數(shù)目。

這樣,文中所定義的廣義I型逐階區(qū)間刪失數(shù)據(jù)就是文獻(xiàn)[1]所定義的I型逐階刪失數(shù)據(jù)的一個特殊類型。事實上,廣義I型逐階區(qū)間刪失數(shù)據(jù)大量地存在于跟蹤調(diào)查等實用場合,對此類數(shù)據(jù)的有效統(tǒng)計分析具有很強的實用價值。

2　有限混合Weibull分布模型

2.1兩組分混合Weibull分布

N組分有限混合Weibull分布模型[16]的概率密度函數(shù)和分布函數(shù)分別是

(1)

I[0，+∞)(x-γk)，

其中，fk(x-γk;αk,λk)=αkλkxλk-1exp(-αk(x-γk)λk)I[0,+∞)(x-γk),αk>0,λk>0,IA(x)是集合A的示性函數(shù),即當(dāng)x∈A時IA(x)=1,否則IA(x)=0。于是,記Ψ=(p1,…,pN-1,α1,λ1,γ1,…,αN,λN,γN)T,參數(shù)θk=(αk,λk,γk)T表示第k(k=1,2,…,N)個子總體(組分)的參數(shù),p=(p1,p2,…,pN)T是混合比,滿足

稱隨機變量X服從N組分3N-1參數(shù)(即Ψ)混合Weibull分布,若其概率密度函數(shù)滿足式(1),記作X～FMWeib(x|Ψ;N)。

則隨機變量X～FMWeib(x|Ψ;N)的幸存函數(shù)為

其中sk(x)=exp(-αk(x-γk)λk)·I[0,+∞)(x-γk),k=1,2,…,N。

文中余下部分規(guī)定參數(shù)N=2,γk=0(k=1,2),即隨機變量X服從兩組分混合Weibull分布,比值p表示第一子總體的占比,將其記作

(2)

從圖3可見混合Weibull分布能夠描述比較復(fù)雜的隨機特性,在圖中所給參數(shù)下,這一混合總體的風(fēng)險率呈現(xiàn)出了先降后升再降的先高后低的變化特征,說明混合Weibull分布具有較好的靈活性,適合描述復(fù)雜場合的壽命規(guī)律,這也是對其研究感興趣的一個重要原因。

圖3　　　　混合Weibull分布0.7,0.08,0.5,0.5,2;2)的幸存函數(shù)、概率密度函數(shù)和風(fēng)險函數(shù)

2.2完全數(shù)據(jù)的似然函數(shù)

設(shè)觀測的n個個體X1,X2,…,Xn壽命獨立同分布(2),即

并設(shè)示性向量Ij=(Ij1,Ij2)T∈2×1,Ij的兩個分量中有且只有一個是1,其余分量均為0,即Ijk=1表示Xj來自混合總體中的第k(k=1,2)個子總體,否則Ijk=0。記I=(I11,I12,I21,I22,…,In1,In2)T∈n×2,則矩陣I表示所有被觀察個體壽命X1,X2,…,Xn的分類示性矩陣,于是,完全數(shù)據(jù)Z={Zj}={j=1,2,…,n}的似然函數(shù)為

(3)

2.3廣義I型逐階區(qū)間刪失數(shù)據(jù)的似然函數(shù)

用X=(Xj)表示所有個體的壽命組成的向量,I=(I1,I2,…,In)T∈n×4m表示示性單元矩陣,其中的第j行給出了Xj的屬性特征,C={C1,C2,…,Cm}表示觀測區(qū)間內(nèi)失效的個體數(shù)目,D={R1,R2,…,Rm}表示從每個觀測區(qū)間右端點隨機刪失的個體數(shù)目。于是Z={Zj}={X,I,C,D}為與Y={C1,C2,…,Cm,R}對應(yīng)的完全數(shù)據(jù),這里R=R1+R2+…+Rm,并且完全數(shù)據(jù)Zj={Xj,Ij,Ci,Di}的概率密度函數(shù)為

再由式(3)可得廣義I型逐階區(qū)間刪失補全后的完全數(shù)據(jù)Z所對應(yīng)的對數(shù)似然函數(shù)為

(logpk+logfk(xj))。

(4)

(5)

在下一部分里,討論已知不完全數(shù)據(jù)Y的情況下,參數(shù)Ψ的最大似然估計的EM算法。

3　基于EM算法參數(shù)Ψ的最大似然估計

3.1模型基本假設(shè)與性質(zhì)

根據(jù)引言中廣義I型逐階區(qū)間刪失數(shù)據(jù)的表述,可知文中所研究模型滿足以下三個基本假設(shè),這些假設(shè)與實際觀察的環(huán)境條件是相統(tǒng)一的。

假設(shè)1觀察時間區(qū)間等間隔。

假設(shè)2每個時間節(jié)點進(jìn)入觀測的個體數(shù)目獨立服從同分布。

假設(shè)3文中討論的廣義I型逐階區(qū)間刪失數(shù)據(jù)依據(jù)圖2所示方案獲得。

基于以上假設(shè),設(shè)與不完全觀測數(shù)據(jù)Y={C1,C2,…,Cm,R}對應(yīng)的完全數(shù)據(jù)為Z={X,I,C,D},C={C1,C2,…,Cm},D={R1,R2,…,Rm},i=1,2,…,m,I=(I1,I2,…,In)T∈n×4m,有如下性質(zhì)成立。

性質(zhì)1區(qū)間失效數(shù)據(jù)C服從多項式分布,滿足

性質(zhì)2刪失數(shù)據(jù)D服從多項式分布,滿足

3.2基于EM算法參數(shù)Ψ的最大似然估計

在兩組分混合Weibull分布及廣義I型逐階區(qū)間刪失假設(shè)下,討論在EM算法的框架下求解參數(shù)Ψ=(p,α1,λ1,α2,λ2)的最大似然估計問題。

設(shè)定初始值Ψ(0),根據(jù)式(5),EM算法的第t次迭代為:

E步假定參數(shù)的第t-1次迭代估計為Ψ(t-1),則第t步的Q函數(shù)為

Q(Ψ;Ψ(t-1))=EΨ(t-1)[logLc(Z;Ψ)|Y]=

(6)

這里k=1,2,記p1=p,p2=1-p,于是

(7)

其中

另一方面,

(8)

其中,

w=c1+c2+…+cm,

于是,根據(jù)式(7)和(8),式(6)的Q-函數(shù)可簡記為

Q(Ψ;Ψ(t-1))=EΨ(t-1)[logLc(Ψ)|Y]=

(9)

定理1Q-函數(shù)關(guān)于參數(shù)Ψ具有一階連續(xù)偏導(dǎo)數(shù)。

Q-函數(shù)關(guān)于參數(shù)pk、αk一階偏導(dǎo)存在且連續(xù)是顯然的,這里不進(jìn)行證明。

定義1若對任意正實數(shù)μ,隨機變量Xμ的數(shù)學(xué)期望EXμ存在,則稱隨機變量X具有連續(xù)階原點矩,或稱隨機分布具有連續(xù)階原點矩。

引理1Weibull分布連續(xù)階原點矩存在,且具有一階連續(xù)可微性。

證明設(shè)隨機變量X服從參數(shù)為α(>0),λ(>0)的Weibull分布,μ∈(0,+∞)為任意指定的正實數(shù),則隨機變量X的μ階原點矩為

(10)

式(10)中第二個等號,是引進(jìn)變換t=αxλ實現(xiàn)的,根據(jù)伽馬函數(shù)Γ(s)(s>0)關(guān)于參數(shù)s具有任意階連續(xù)導(dǎo)數(shù)的性質(zhì),以及復(fù)合函數(shù)的可導(dǎo)性,可知服從Weibull分布的隨機變量具有任意階的原點矩,且對階次μ具有任意階連續(xù)偏導(dǎo)數(shù),則關(guān)于μ的一階偏導(dǎo)為

其中?！?·)表示Γ函數(shù)關(guān)于自變量·的導(dǎo)數(shù)。此外還有

成立。

M步對式(9)中Q-函數(shù)的參數(shù)Ψ求偏導(dǎo),得如下方程組:

(11)

(12)

(Ⅱ)

(13)

其中,

3.3算例1

下面通過仿真算例說明EM算法的有效性。這里仿真數(shù)據(jù)生成方案如下:

第1步,指定分組數(shù)m及觀測時間點0

第2步,獨立生成m組同泊松分布P(μ·Δ)的樣本數(shù)據(jù)ni,i=1,2，…,m。

第3步,對于每個i按概率p生成長度為ni的1、0數(shù)據(jù)串,其中數(shù)字1表示該觀測數(shù)據(jù)來自第一個子總體,數(shù)字0表示該觀測數(shù)據(jù)來自第二個子總體,并記錄該數(shù)據(jù)串中數(shù)字1的個數(shù)為ni1,數(shù)字0的個數(shù)為ni0,p為第一個子總體的混合比。

第4步,獨立生成m組容量為ni1,i=1,2,…,m的服從Weibull分布f1(x;α1,λ1)的數(shù)據(jù)xi1j,i=1,2,…,m;j=1,2,…,ni1,再獨立生成m組容量為ni0,i=1,2,…,m的服從Weibull分布f2(x;α2,λ2)的數(shù)據(jù)xi0k,i=1,2,…,m;k=1,2,…,ni0。

第5步,按照圖2所示方法統(tǒng)計分組觀測數(shù)據(jù)Y。

現(xiàn)在,考慮一個特別的算例。令λ1=λ2=1且已知,此時,仿真數(shù)據(jù)總體分布服從兩組分混合指數(shù)分布,仿真總體的其他參數(shù)為α1=0.2、α2=0.04、p=0.3、m=10、μ=80、Δ=1,按此參數(shù)生成的時間區(qū)間i內(nèi)的失效數(shù)據(jù)分別為60、66、45、24、25、17、10、6、5和2,其中,樣本容量n=759,失效總數(shù)w=260,時間區(qū)間i表示觀測區(qū)間((i-1)Δ,iΔ],這里Δ=1。圖4給出了仿真總體的概率密度函數(shù)、生存函數(shù)和風(fēng)險函數(shù)的圖像,從中可見混合指數(shù)分布的風(fēng)險函數(shù)不再是常數(shù)?，F(xiàn)在基于上文仿真數(shù)據(jù),根據(jù)式(12)、(13)進(jìn)行參數(shù)估計。

圖4　　　　混合Exponential分布的幸存函數(shù)、概率密度函數(shù)和風(fēng)險函數(shù)

圖5顯示了在不同初值情況下,EM算法參數(shù)的

圖5　EM算法下參數(shù)估計值與迭代步數(shù)間的關(guān)系

估計的收斂特性,可以看到,不同初值產(chǎn)生的最終的收斂狀況幾乎沒有差別;圖6顯示了200步之內(nèi)的迭代,可見不同初值對估計結(jié)果在最初的計算中是有較大影響的。從圖5和6還可發(fā)現(xiàn),隨著迭代步數(shù)的增加,混合比p越來越大,參數(shù)α1、α2盡管不同,卻是越來越小,這說明以增加迭代步數(shù)為代價獲得更好估計結(jié)果是不合適的。

圖6　200步以內(nèi)不同初值EM算法的收斂效果

產(chǎn)生這一現(xiàn)象的原因在于文中所提的廣義I型逐階區(qū)間刪失數(shù)據(jù)自身沒能提供分類信息,而從圖6的細(xì)節(jié)上會發(fā)現(xiàn)參數(shù)p的估值都有一個先降后升的過程,即p的估計值有最小值;另一方面,對基于有限混合模型的廣義I型逐階區(qū)間刪失數(shù)據(jù)估計的一個重要興趣在于:希望知道短壽子總體的壽命分布信息和占比,這促使提出第4部分所描述的基于加權(quán)絕對偏差信息準(zhǔn)則的改進(jìn)EM算法,圖5和6也顯示了加權(quán)偏差信息準(zhǔn)則的良好收斂特性。

4　基于加權(quán)絕對偏差信息準(zhǔn)則的EM算法

4.1改進(jìn)的EM算法

實際上,在質(zhì)量監(jiān)控領(lǐng)域中,對于廣義I型逐階區(qū)間刪失統(tǒng)計數(shù)據(jù)的興趣在于早期發(fā)現(xiàn)質(zhì)量問題,以便實施質(zhì)量控制。因此,在參數(shù)估計中,希望知道混合比p的合理的最小估值、短壽子總體壽命的合理最短估值。這種觀點是基于這樣的事實:短壽命子總體的最小占比值的估值代表著人們最樂觀的估計,當(dāng)它都不可接受時,糟糕的情形就更危險了,所以,它適于作為風(fēng)險評估的指標(biāo);另外,在EM算法估計的不斷迭代中壽命參數(shù)不斷變小,意味著子總體平均壽命的延長,所以希望知道子總體壽命的最短估計值。基于這樣的事實,文中提出基于加權(quán)絕對偏差信息準(zhǔn)則的改進(jìn)EM算法對參數(shù)進(jìn)行估計。

首先定義加權(quán)絕對偏差信息,并將其作為EM算法的終止準(zhǔn)則。

定義2在有限混合Weibull分布情形下,對于文中所提的廣義I型逐階區(qū)間刪失數(shù)據(jù),定義加權(quán)絕對偏差信息(WeightedAbsoluteDeviationInformation,WADI)為μ,即

其中區(qū)間失效數(shù)據(jù)Ci及其分布滿足性質(zhì)1。

假設(shè)4在觀測時間內(nèi),來自短壽命子總體的個體幾乎都發(fā)生失效;一般長壽命子總體的平均壽命是短壽命子總體平均壽命的三倍以上。

基于加權(quán)絕對偏差信息準(zhǔn)則的EM算法:

第2步,加權(quán)絕對偏差信息準(zhǔn)則(WADIC)為以加權(quán)絕對偏差信息的相對誤差絕對值小于指定閥值為迭代終止條件。

第3步,實施EM算法的估計,重復(fù)迭代,直到獲得可接受的估計結(jié)果。

4.2算例2(算例1續(xù))

圖7　迭代步數(shù)與各被估參數(shù)相對誤差絕對值的關(guān)系

圖8　估計分布函數(shù)的擬合效果與迭代步數(shù)的關(guān)系

4.3算例3

(iii)比較這十三個估計值,選出μ最小的矩形,作為下一步迭代時搜索參數(shù)估計值的區(qū)域。

重復(fù)(i)～(iii)步,直到這些嵌套區(qū)域內(nèi)估計總體的μ的相對誤差絕對值小于閥值0.000 1停止迭代,從而獲得參數(shù)λ1、λ2、p、α1、α2的估計值。

經(jīng)過迭代計算,得到參數(shù)估計值分別為α1=0.188 4、α2=0.035 0、p=0.112 3、λ1=1.798 25、λ2=0.898 3。圖9給出了仿真總體的概率密度函數(shù)和估計總體概率密度函數(shù)的對比圖,可見所提算法對觀測數(shù)據(jù)的擬合效果較好,這種搜索定位法節(jié)約計算資源,使得估計運算的效率較高。

圖9　仿真總體和估計總體的概率分布函數(shù)

5　結(jié)束語

文中討論了兩組分混合Weibull分布在廣義I型逐階區(qū)間刪失情形下參數(shù)估計的EM算法的估計效率問題,發(fā)現(xiàn)當(dāng)刪失比較大時,EM算法存在過度迭代現(xiàn)象,使得估計結(jié)果偏離數(shù)據(jù)來源的分布,為此提出了結(jié)合加權(quán)絕對偏差信息準(zhǔn)則進(jìn)行參數(shù)估計的方法,仿真算例說明了文中所提方法的有效性。

[1]TONG Ng H K,WANG Z. Statistical estimation for the parameters of Weibull distribution based on progressively type-I interval censored sample[J]. Journal of Statistical Computation and Simulation, 2009, 79(2): 145-159.

[2]蘇錦霞, 張藝贏, 田麗娜. Ⅰ型逐階區(qū)間刪失Weibull數(shù)據(jù)的統(tǒng)計分析[J]. 蘭州大學(xué)學(xué)報: 自然科學(xué)版, 2011, 47(5): 109-114, 119.

[3]MOKHTARI E B, RAD A H, YOUSEFZADEH F. Inference for Weibull distribution based on progressively type-II hybrid censored data[J]. Journal of Statistical Planning and Inference, 2011, 141(8): 2824-2838.

[4]JOARDER A, KRISHNA H, KUNDUC D. Inferences on Weibull parameters with conventional type-I censoring[J]. Computational Statistics and Data Analysis, 2011, 55(1): 1-11.

[5]GHITANY M E, TUAN V K, BALAKRISHNANC N. Likelihood estimation for a general class of inverse exponentiated distributions based on complete and progressively censored data[J]. Journal of Statistical Computation and Simulation, 2014, 84(1): 96-106.

[6]鄭明, 楊藝, 鄭宇. 基于分組數(shù)據(jù)的Weibull分布的參數(shù)估計[J]. 高校應(yīng)用數(shù)學(xué)學(xué)報A輯:中文版, 2003, 18(3): 303-310.

[7]吳耀國, 周杰, 王柱, 等. 隨機刪失數(shù)據(jù)下基于EM算法的Weibull分布參數(shù)估計[J].四川大學(xué)學(xué)報: 自然科學(xué)版, 2005, 42(5): 910-913.

[8]任瑞, 周秀輕. 逐步Ⅰ型區(qū)間刪失數(shù)據(jù)下的參數(shù)估計[J]. 南京師大學(xué)報: 自然科學(xué)版, 2011, 34(3): 7-12.

[9]張頌, 王德輝. 熵?fù)p失下定數(shù)Progressive刪失情形Weibull分布尺度參數(shù)的估計[J]. 吉林大學(xué)學(xué)報: 理學(xué)版, 2012, 50(2): 219-226.

[10]LIN C T, CHOU C C, HUANG Y L. Inference for the Weibull distribution with progressive hybrid censoring[J]. Computational Statistics and Data Analysis, 2012, 56(3): 451-467.

[11]PARK B J, LORD D. Application of finite mixture models for vehicle crash data analysis[J]. Accident Analysis and Prevention, 2009, 41(4): 683-691.

[12]蔣卉, 湯銀才. 混合Weibull分布參數(shù)估計的ECM算法[J]. 系統(tǒng)科學(xué)與數(shù)學(xué), 2010, 30(1): 79-88.

[13]木拉提·吐爾德, 胡錫健. EM算法在刪失數(shù)據(jù)分布和混合分布參數(shù)估計中的應(yīng)用[J]. 統(tǒng)計與決策, 2011, 339(15): 161-164.

[14]張曉勤, 王煜, 盧殿軍. 混合指數(shù)威布爾分布的參數(shù)估計[J]. 河南大學(xué)學(xué)報: 自然科學(xué)版, 2012, 42(3): 230-233.

[15]田玉柱, 田茂再, 陳平. 數(shù)據(jù)分組和右刪失下混合廣義指數(shù)分布的參數(shù)估計[J]. 應(yīng)用概率統(tǒng)計, 2012, 28(6): 561-571.

[16]MC LACHLAN G, PEEL D. Finite mixture models[M]. USA,New York: John Wiley & Sons, Inc, 2000.

(編輯王冬)

Parameters estimation of generalized progressive type-I interval-censored mixture Weibull-distributed data

YUANYanhua

(School of Sciences, Heilongjiang University of Science & Technology, Harbin 150022, China)

This paper is a response to the previously insufficient study of the method designed for the parameter estimation for the generalized progressive type-I interval-censored data based on the finite mixture. The study building on finite mixture Weibull-distribution model looks at the effectiveness of Expectation-Maximization(EM) algorithm used for the parameter estimation for generalized progressive type-I interval-censored data and produces a modified EM algorithm. The modification starts with giving the EM algorithm used for the parameter estimation, followed the simulation examples to explain the excessive iteration produced by EM algorithm when used for the parameter estimation for the generalized progressive type-I interval-censored data, and culminates in the necessity for stopping weighted absolute deviance information criterion involved in EM algorithm. The modified EM algorithm is free from the drawback inherent in EM algorithm incapable of determining the time right for stopping iteration when used to perform parameter estimation, thus allowing for a quick production of satisfactory estimators, following an appropriate selection of the initial value. The simulation verifies the viability of the algorithm.

generalized progressive type-I interval censoring; weighted absolute deviation information; finite mixture Weibull distribution; EM algorithm

2013-12-11;

2014-03-25

苑延華(1969-),女,遼寧省本溪人,教授,博士研究生,研究方向:概率統(tǒng)計、運籌學(xué)與控制論,E-mail:yuanhua-69@sina.com。

10.3969/j.issn.2095-7262.2014.03.021

O213

2095-7262(2014)03-0323-09

A