Geometric分布與壽命分布復(fù)合族參數(shù)估計(jì)的MM算法

熊朝松

(云南師范大學(xué) 數(shù)學(xué)學(xué)院,云南 昆明 650500)

0 引言

壽命分布在可靠性分析中扮演著重要角色.在處理生活中記錄的壽命數(shù)據(jù)時(shí),最直接的方法就是用壽命分布進(jìn)行擬合分析,一個(gè)好的擬合結(jié)果往往能起到重要的指導(dǎo)作用.然而隨著科技的發(fā)展和時(shí)代的進(jìn)步,現(xiàn)實(shí)生活中觀測(cè)到的壽命數(shù)據(jù)越來(lái)越多,結(jié)構(gòu)也越來(lái)越復(fù)雜.經(jīng)典的壽命分布如Exponential分布和Weibull分布等已經(jīng)無(wú)法完全適應(yīng)多種多樣的壽命數(shù)據(jù),在大多數(shù)情況下并不能得到一個(gè)較好的擬合結(jié)果.因此,為了減小在擬合分析中產(chǎn)生的偏差,開(kāi)發(fā)新的壽命分布成為可靠性統(tǒng)計(jì)領(lǐng)域的重點(diǎn)研究方向.學(xué)者們通過(guò)復(fù)合一個(gè)經(jīng)典的壽命分布和不含零的離散分布,研究了大量的模擬壽命數(shù)據(jù)的多參數(shù)分布.Adamidis和Loukas[1]首先將Exponential分布和Geometric分布復(fù)合,得到了一個(gè)新的具有下降風(fēng)險(xiǎn)率的Exponential-Geometric(EG)分布,并詳細(xì)討論了包括生存函數(shù)、風(fēng)險(xiǎn)函數(shù)、似然推斷等在內(nèi)的統(tǒng)計(jì)性質(zhì).EG分布的提出,為開(kāi)發(fā)新壽命分布提供了一個(gè)新的思路.受到Adamidis和Loukas的啟發(fā),Barreto-Souza等[2]和顧蓓青等[3]開(kāi)發(fā)了Weibull-Geometric(WG)分布;姚惠等[4]開(kāi)發(fā)了Pareto-Geometric(PG)分布;高艷紅和周秀輕[5],Okorie等[6]開(kāi)發(fā)了Rayleigh-Geometric(RG)分布;Alkarni[7]歸納總結(jié)了Geometric分布與常用壽命分布的復(fù)合形式,將這類分布統(tǒng)稱為Geometric分布與壽命分布復(fù)合族,并進(jìn)一步研究了這類分布族的相關(guān)統(tǒng)計(jì)性質(zhì)，概率密度函數(shù)和分布函數(shù)形式,可靠性和失效率函數(shù)的顯示公式等.在似然推斷中,Alkarni[7]建議使用EM算法求解這類復(fù)合族分布的極大似然估計(jì).但由于EM算法涉及引入潛變量問(wèn)題,對(duì)于部分初學(xué)者來(lái)說(shuō)不容易理解.因此開(kāi)發(fā)一個(gè)簡(jiǎn)單易懂的參數(shù)估計(jì)方法是有必要的.

與EM算法類似,Minorization-maximization(MM)算法同樣也是一種迭代優(yōu)化算法,其工作原理是通過(guò)找到一個(gè)替代函數(shù)驅(qū)動(dòng)目標(biāo)函數(shù)最大化.值得注意的是,MM算法包含了EM算法,即每個(gè)EM算法都是MM算法,反之不然.且MM算法在算法操作上與EM算法有所區(qū)別,MM算法不需要像EM算法那樣引入潛變量來(lái)尋找替代函數(shù),而是直接利用目標(biāo)函數(shù)的極小化函數(shù)作為替代函數(shù).由于MM算法具有簡(jiǎn)單的結(jié)構(gòu)原理和優(yōu)秀的估計(jì)性質(zhì),因此被廣泛應(yīng)用于統(tǒng)計(jì)學(xué)的各個(gè)領(lǐng)域[8-10].

本文引入MM算法來(lái)計(jì)算Geometric分布與壽命分布復(fù)合族的極大似然估計(jì),推導(dǎo)了這類分布族的MM算法參數(shù)估計(jì)式,隨后將本文提出的方法應(yīng)用于幾個(gè)具體的分布中,最后通過(guò)隨機(jī)模擬實(shí)驗(yàn)驗(yàn)證MM算法估計(jì)的有效性.

1 MM算法

假設(shè)l(θ)是待最大化的目標(biāo)函數(shù),θ∈Θ表示未知的參數(shù)向量,令l(θ)的最大值點(diǎn)為

MM算法的構(gòu)造方式主要包括兩個(gè)步驟,具體內(nèi)容如下.

第一步:Minorization(M)步,主要利用極小化功能找到目標(biāo)函數(shù)的替代函數(shù)Q(θ|θ(t)),或稱為極小化函數(shù),使其滿足

Q(θ|θ(t))≤l(θ),?θ,θ(t)∈Θ;Q(θ(t)|θ(t))=l(θ(t)),

第二步:Maximization(M)步,通過(guò)最大化Q(θ|θ(t))將θ(t)更新為θ(t+1),可以得到

l(θ(t+1))≥Q(θ(t+1)|θ(t))≥Q(θ(t)|θ(t))=l(θ(t)).

MM算法在每次迭代計(jì)算時(shí)都會(huì)增加目標(biāo)函數(shù)值,并且這種上升性質(zhì)為MM算法提供了顯著的數(shù)值穩(wěn)定性.

MM算法的關(guān)鍵步驟就在于第一步中找到合適的替代函數(shù),在此本文介紹一種常用的極小化方法,該方法由Zhou和Lange[11]提出,Tian等[10]對(duì)其做出了詳細(xì)介紹.令r(x)和s(x)是在相同支撐集上的兩個(gè)密度函數(shù),即

{x|r(x)>0}={x|s(x)>0}.

r(x)和s(x)之間的Kullback-Leibler(KL)距離被定義為

對(duì)應(yīng)的離散形式為

在上述公式中，令n=2,(r1,r2)=(δ0,1-δ0)，(s1,s2)=(δ,1-δ),則有

δ0logδ0+(1-δ0)log(1-δ0)≥
δ0logδ+(1-δ0)log(1-δ),

其中δ,δ0∈(0,1),當(dāng)且僅當(dāng)δ=δ0時(shí)等號(hào)成立.將上述不等式重新排列,即有

(1)

2 Geometric分布與壽命分布復(fù)合族的MM算法估計(jì)

2.1 Geometric分布與壽命分布復(fù)合族

定義1給定一個(gè)隨機(jī)變量Z服從一個(gè)零截?cái)郍eometric分布,其概率質(zhì)量函數(shù)(pmf)滿足

fZ(z|p)=(1-p)pz-1,0

假設(shè)T=(T1,…,TZ)是獨(dú)立同分布,概率密度函數(shù)(pdf)為

fTi(x|α)=fT(x|β)|β=(β1,…,βm),m≥1,x,β∈R+

的隨機(jī)變量,其中對(duì)于任意i=1,2,…,隨機(jī)變量Z和Ti是相互獨(dú)立的.令隨機(jī)變量X=min(T1,…,TZ),則根據(jù)上述條件,可以得到X|Z=z的pdf為

同時(shí)可以獲得X和Z的聯(lián)合概率密度函數(shù)為

其中FT(x|β)表示T的累積分布函數(shù)(cdf).通過(guò)計(jì)算容易得到X的邊際pdf和cdf分別為

(2)

和

(3)

如果一個(gè)隨機(jī)變量X的pdf和cdf分別為公式(2)和(3)的形式,則稱X服從Geometric分布與壽命分布復(fù)合族,記作X～GL(p,β).

2.2 極大似然估計(jì)的MM算法

l(p,β|Xobs)=l0(p,β|Xobs)+l1(p,β|Xobs),

其中

觀察到,l1(p,β|Xobs)含有與-log(1-δ)類似的結(jié)構(gòu),則考慮通過(guò)極小化l1(p,β|Xobs)來(lái)獲取l(p,β|Xobs)的替代函數(shù).

已知0

0

將δ=p[1-FT(xi|β)]和δ0=p(t)[1-FT(xi|β(t))]代入公式(1),可以得到l1(p,β|Xobs)的極小化函數(shù)為

其中c是與參數(shù)(p,β)無(wú)關(guān)的常數(shù)項(xiàng).結(jié)合l0(p,β|Xobs)與Q1(p,β|p(t),β(t))可以得到對(duì)數(shù)似然函數(shù)l(p,β|Xobs)的替代函數(shù)為

其中

(4)

3 應(yīng)用

將本文提出的方法應(yīng)用于四個(gè)常見(jiàn)的Geometric分布與壽命分布復(fù)合族分布,分別為Exponential-Geometric分布,Weibull-Geometric分布,Pareto-Geometric分布和Rayleigh-Geometric分布.為了方便推導(dǎo)出基于極大似然的MM算法估計(jì)公式,首先在此介紹Tian等[12]提出的LOG-BETA函數(shù)族和LOG-GAMMA函數(shù)族.

定義2假設(shè)一個(gè)函數(shù)在定義域[0,1]上的表達(dá)式為

g(θ)=c+alog(θ)+blog(1-θ),θ∈[0,1],

其中,c是與參數(shù)θ不相關(guān)的常數(shù)項(xiàng),且有a,b≥0.則稱此函數(shù)屬于LOG-BETA函數(shù)族,記作g(θ)∈LB(θ).log(θ)和log(1-θ)分別為L(zhǎng)OG-BETA函數(shù)族的兩個(gè)基本組裝函數(shù).

注意到,當(dāng)a,b>0時(shí),g(θ)是嚴(yán)格凹函數(shù),則可以得到

(5)

定義3假設(shè)一個(gè)函數(shù)在定義域正實(shí)數(shù)集R+上的表達(dá)式為

g(θ)=c+alog(θ)+b(-θ),θ∈R+,

其中,c是與參數(shù)θ不相關(guān)的常數(shù)項(xiàng),且有a,b≥0.則稱此函數(shù)屬于LOG-GAMMA函數(shù)族,記作g(θ)∈LG(θ).log(θ)和-θ分別為L(zhǎng)OG-GAMMA函數(shù)族的兩個(gè)基本組裝函數(shù).

同樣地,當(dāng)a,b>0時(shí),g(θ)是嚴(yán)格凹函數(shù),則可以得到

(6)

基于上述的LOG-BETA函數(shù)族和LOG-GAMMA函數(shù)族的定義與性質(zhì),下面詳細(xì)討論MM算法在求解幾種常見(jiàn)的復(fù)合壽命分布的極大似然估計(jì)時(shí)的應(yīng)用.

3.1 Exponential-Geometric分布

Adamidis等[1]提出的Exponential-Geometric分布的概率密度函數(shù)為

f(x|p,β)=β(1-p)e-βx(1-pe-βx)-2,x>0，

其中00.則對(duì)應(yīng)的對(duì)數(shù)似然函數(shù)為

利用本文2.2節(jié)的方法,l(p,β|Xobs)的替代函數(shù)為

Q(p,β|p(t),β(t))=Q(p|p(t),β(t))+Q(β|p(t),β(t))+c,

其中

(7)

觀察到,MM算法將一個(gè)兩參數(shù)的優(yōu)化問(wèn)題轉(zhuǎn)變成了兩個(gè)單參數(shù)優(yōu)化問(wèn)題之和,其中關(guān)于p的函數(shù)屬于LOG-BETA函數(shù)族,關(guān)于β的函數(shù)屬于LOG-GAMMA函數(shù)族.

分別令dQ(p|p(t),β(t))/dp=0和dQ(β|p(t),β(t))/dβ=0,同時(shí)利用公式(5)和(6),即可得到參數(shù)p和β的MM算法估計(jì)式.

3.2 Weibull-Geometric分布

Barreto-Souza等[2]將Weibull分布與Geometric分布復(fù)合,得到了三參數(shù)的Weibull-Geometric分布,其概率質(zhì)量函數(shù)滿足

f(x|p,β,α)=αβα(1-p)xα-1exp[-(βx)α]·{1-pexp[-(βx)α]}-2,x>0，

其中00,α>0.令上述的pmf中的參數(shù)β=1/β*,即可得到顧蓓青等[3]考慮的Weibull-Geometric分布.

根據(jù)Weibull-Geometric分布的概率密度函數(shù),對(duì)應(yīng)的對(duì)數(shù)似然函數(shù)可寫(xiě)為

利用本文2.2節(jié)的方法,l(p,β,α|Xobs)的替代函數(shù)為

Q(p,β,α|p(t),β(t),α(t))=Q(p|p(t),β(t),α(t))+Q(β,α|p(t),β(t),α(t))+c,

其中

(8)

容易看出,Q(p|p(t),β(t),α(t))∈LB(p),利用公式(5)可以得到參數(shù)p的估計(jì)式.參數(shù)(β,α)的MM算法迭代公式為

其中α(t+1)為非線性方程

的解.

3.3 Pareto-Geometric分布

姚惠等[4]研究了定義在x>1上的Pareto-Geometric分布.而在現(xiàn)實(shí)中觀測(cè)到的壽命數(shù)據(jù)大多定義在x>0上,因此根據(jù)定義1,可以得到Pareto-Geometric分布的概率密度函數(shù)為

f(x|p,β)=β(1-p)(1+x)-(β+1)[1-p(1+x)-β]-2,x>0,

其中00.故對(duì)應(yīng)的對(duì)數(shù)似然函數(shù)為

利用本文2.2節(jié)的方法,其替代函數(shù)可以寫(xiě)為

Q(p,β|p(t),β(t))=Q(p|p(t),β(t))+Q(β|p(t),β(t))+c,

其中

令dQ(p|p(t),β(t))/dp=0,dQ(β|p(t),β(t))/dβ=0,同時(shí)利用公式(5)和(6),即可得到Pareto-Geometric分布參數(shù)p和β的MM算法估計(jì)式.

3.4 Rayleigh-Geometric分布

高艷紅等[5]提出的Rayleigh-Geometric分布的概率密度函數(shù)為

f(x|p,β)=2β(1-p)xe-βx2(1-pe-βx2)-2,x>0,

其中00.對(duì)應(yīng)的對(duì)數(shù)似然函數(shù)為

l(p,β|Xobs)替代函數(shù)為

Q(p,β|p(t),β(t))=Q(p|p(t),β(t))+Q(β|p(t),β(t))+c,

其中

(10)

令dQ(p|p(t),β(t))/dp=0和dQ(β|p(t),β(t))/dβ=0,同時(shí)利用公式(5)和(6),即可得到Rayleigh-Geometric分布參數(shù)p和β的MM算法估計(jì)式.

4 隨機(jī)模擬

本節(jié)基于第3節(jié)中提到的例子進(jìn)行了數(shù)值實(shí)驗(yàn),以此評(píng)估所提出的MM算法的實(shí)際性能.四組實(shí)驗(yàn)分別生成樣本量為100和500的EG分布、WG分布、PG分布和RG分布獨(dú)立隨機(jī)樣本,并在不同的參數(shù)值組合下,利用本文提出的MM算法計(jì)算參數(shù)的極大似然估計(jì)值.將每組試驗(yàn)重復(fù)進(jìn)行1000次,計(jì)算出參數(shù)的極大似然估計(jì)平均值和估計(jì)的均方誤差.

第一組實(shí)驗(yàn)中,分別設(shè)置EG分布的參數(shù)(p,β)真值為(0.1,0.1),(0.9,0.1),(0.1,1)和(0.9,1).利用3.1節(jié)定義的MM算法來(lái)計(jì)算(p,β)的極大似然估計(jì),隨機(jī)模擬結(jié)果如表1所列.

表1 EG分布參數(shù)的極大似然估計(jì)

第二組實(shí)驗(yàn)中,分別設(shè)置WG分布的參數(shù)(p,β,α)真值為(0.2,0.2,2)和(0.8,2,3).利用3.2節(jié)所定義的MM算法來(lái)計(jì)算(p,β,α)的極大似然估計(jì),隨機(jī)模擬結(jié)果如表2所列.

表2 WG分布參數(shù)的極大似然估計(jì)

第三組實(shí)驗(yàn)中,分別設(shè)置PG分布的參數(shù)(p,β)真值為(0.3,0.3),(0.7,0.3),(0.3,3)和(0.7,3).利用3.3節(jié)定義的MM算法來(lái)計(jì)算(p,β)的極大似然估計(jì),隨機(jī)模擬結(jié)果如表3所列.

表3 PG分布參數(shù)的極大似然估計(jì)

第四組實(shí)驗(yàn)中,分別設(shè)置RG分布的參數(shù)(p,β)真值為(0.4,0.4),(0.6,0.4),(0.4,4)和(0.6,4).利用3.4節(jié)定義的MM算法來(lái)計(jì)算(p,β)的極大似然估計(jì),隨機(jī)模擬結(jié)果如表4所列.

表4 RG分布參數(shù)的極大似然估計(jì)

從表1至表4的結(jié)果可以看出,針對(duì)于4個(gè)不同的分布,在不同的樣本量和參數(shù)真值的組合下,各個(gè)參數(shù)的極大似然估計(jì)值均能收斂到參數(shù)真值,同時(shí)隨著樣本量的增加,估計(jì)的偏差和均方誤差越來(lái)越小,表明參數(shù)估計(jì)的效果和穩(wěn)定性越來(lái)越好.

5 結(jié)束語(yǔ)

本文首先給出了計(jì)算Geometric分布與壽命分布復(fù)合族的參數(shù)極大似然估計(jì)的MM算法,為這類壽命分布族的參數(shù)估計(jì)提供了一個(gè)新的方案.與現(xiàn)有的EM算法不同的是,MM算法是完全基于觀測(cè)數(shù)據(jù)的,不再需要引入潛變量.此外,MM算法將多維優(yōu)化問(wèn)題轉(zhuǎn)變成了低維優(yōu)化問(wèn)題之和,極大地減少了計(jì)算難度和成本.隨后本文將提出的方法應(yīng)用于EG分布、WG分布、PG分布和RG分布的極大似然估計(jì)中,得到了各個(gè)分布參數(shù)的MM算法估計(jì)式.隨機(jī)模擬結(jié)果表明MM算法估計(jì)結(jié)果準(zhǔn)確有效,并且隨著樣本量的增加,估計(jì)的效果也越來(lái)越好.本文提出的方法為Geometric分布與壽命分布復(fù)合族的參數(shù)估計(jì)提供了一個(gè)統(tǒng)一框架,能對(duì)后續(xù)開(kāi)發(fā)的新壽命分布的參數(shù)估計(jì)有所幫助.