苑延華

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

?

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

苑延華

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

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

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

0　引　言

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

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

1　刪失數據的定義

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

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

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

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

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

2　有限混合Weibull分布模型

2.1兩組分混合Weibull分布

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

(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的示性函數,即當x∈A時IA(x)=1,否則IA(x)=0。于是,記Ψ=(p1,…,pN-1,α1,λ1,γ1,…,αN,λN,γN)T,參數θk=(αk,λk,γk)T表示第k(k=1,2,…,N)個子總體(組分)的參數,p=(p1,p2,…,pN)T是混合比,滿足

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

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

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

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

(2)

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

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

2.2完全數據的似然函數

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

并設示性向量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的分類示性矩陣,于是,完全數據Z={Zj}={j=1,2,…,n}的似然函數為

(3)

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

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

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

(logpk+logfk(xj))。

(4)

(5)

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

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

3.1模型基本假設與性質

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

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

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

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

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

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

性質2刪失數據D服從多項式分布,滿足

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

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

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

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

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

(6)

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

(7)

其中

另一方面,

(8)

其中,

w=c1+c2+…+cm,

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

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

(9)

定理1Q-函數關于參數Ψ具有一階連續(xù)偏導數。

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

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

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

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

(10)

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

其中Γ′(·)表示Γ函數關于自變量·的導數。此外還有

成立。

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

(11)

(12)

(Ⅱ)

(13)

其中,

3.3算例1

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

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

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

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

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

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

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

圖4　　　　混合Exponential分布的幸存函數、概率密度函數和風險函數

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

圖5　EM算法下參數估計值與迭代步數間的關系

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

圖6　200步以內不同初值EM算法的收斂效果

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

4　基于加權絕對偏差信息準則的EM算法

4.1改進的EM算法

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

首先定義加權絕對偏差信息,并將其作為EM算法的終止準則。

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

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

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

基于加權絕對偏差信息準則的EM算法:

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

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

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

圖7　迭代步數與各被估參數相對誤差絕對值的關系

圖8　估計分布函數的擬合效果與迭代步數的關系

4.3算例3

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

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

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

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

5　結束語

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

[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數據的統(tǒng)計分析[J]. 蘭州大學學報: 自然科學版, 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]鄭明, 楊藝, 鄭宇. 基于分組數據的Weibull分布的參數估計[J]. 高校應用數學學報A輯:中文版, 2003, 18(3): 303-310.

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

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

[9]張頌, 王德輝. 熵損失下定數Progressive刪失情形Weibull分布尺度參數的估計[J]. 吉林大學學報: 理學版, 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分布參數估計的ECM算法[J]. 系統(tǒng)科學與數學, 2010, 30(1): 79-88.

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

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

[15]田玉柱, 田茂再, 陳平. 數據分組和右刪失下混合廣義指數分布的參數估計[J]. 應用概率統(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)計、運籌學與控制論,E-mail:yuanhua-69@sina.com。

10.3969/j.issn.2095-7262.2014.03.021

O213

2095-7262(2014)03-0323-09

A