何朝兵, 劉華文

(1.安陽師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院, 河南 安陽 455000; 2.山東大學(xué) 數(shù)學(xué)學(xué)院, 濟南 250100)

?

帶有不完全信息隨機截尾試驗下伽瑪分布尺度參數(shù)的點估計

何朝兵1*, 劉華文2

(1.安陽師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院, 河南 安陽 455000; 2.山東大學(xué) 數(shù)學(xué)學(xué)院, 濟南 250100)

首先通過添加數(shù)據(jù)得到了帶有不完全信息隨機截尾試驗下伽瑪分布的完全數(shù)據(jù)似然函數(shù),然后分別利用EM算法和MCMC方法對尺度參數(shù)進行了估計,最后進行了隨機模擬試驗,結(jié)果表明尺度參數(shù)點估計的精度比較高.

完全數(shù)據(jù)似然函數(shù); EM算法; 滿條件分布; MCMC方法; Gibbs抽樣

伽瑪分布是一類應(yīng)用廣泛的連續(xù)型分布,尤其在水文、氣象、海洋和可靠性領(lǐng)域,并且伽瑪分布與泊松分布、威布爾分布、卡方分布等分布都有著密切的關(guān)系.關(guān)于伽瑪分布的研究可參看文獻[1-6].近些年來,關(guān)于隨機截尾試驗的研究比較多,帶有不完全信息隨機截尾試驗(random censoring test model with incomplete information),簡稱IIRCT,由文獻[7]首次研究,接著文獻[8-18]深入研究了IIRCT下壽命分布的參數(shù)估計,但對IIRCT下伽瑪分布的參數(shù)估計,卻很少有文獻進行研究.本文通過添加數(shù)據(jù)得到了IIRCT下伽瑪分布的完全數(shù)據(jù)似然函數(shù),給出了由EM算法得到的參數(shù)的迭代公式,接著利用MCMC方法得到了參數(shù)的Gibbs樣本,把Gibbs樣本的均值作為參數(shù)的貝葉斯估計,進行了隨機模擬試驗,結(jié)果表明參數(shù)的EM估計和MCMC估計的精度都較高.

1帶有不完全信息隨機截尾試驗(IIRCT)

設(shè)X1,X2,…是獨立同分布的非負(fù)隨機變量序列,其分布函數(shù)為F(x;θ)=P(Xi≤x),密度函數(shù)為f(x;θ),θ是未知參數(shù).又設(shè)Y1,Y2,…是獨立的非負(fù)隨機變量序列,分布函數(shù)分別為G1(y),G2(y),…，密度函數(shù)分別為g1(y),g2(y),…且gi(y)與參數(shù)θ無關(guān).{Xi}與{Yi}獨立.

為估計參數(shù)θ,n個樣品的觀察數(shù)據(jù){Zi,1≤i≤n}如下:

1) 當(dāng)Xi≤Yi時,有兩種情況:Xi以概率ai立即顯示,此時取Zi=Xi;Xi以概率1-ai不被顯示,此時取Zi=Yi.稱ai為失效顯示概率.

2) 當(dāng)Xi>Yi時,取Zi=Yi.

引入隨機變量αi,βi,i=1,2,…,n.

若Xi≤Yi,αi=1;若Xi>Yi,αi=0.

Xi≤Yi,且未被顯示時,βi=0;其它情況,βi=1.

綜上所述,有

P(βi=1|αi=1)=ai,

P(βi=0|αi=1)=1-ai.

設(shè)Zi的觀察值為zi,則基于數(shù)據(jù){(zi,αi,βi),1≤i≤n}的似然函數(shù)為：

[g(zi)F(zi;θ)(1-ai)]αi(1-βi)·

(1)

2IIRCT下伽瑪分布參數(shù)的點估計

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

假設(shè)IIRCT下產(chǎn)品壽命Xi～Ga(b,θ),且形狀參數(shù)b已知,下面對尺度參數(shù)θ進行估計.

由(1)式可以看出,伽瑪分布基于截尾數(shù)據(jù)的似然函數(shù)比較復(fù)雜,為了方便進行參數(shù)估計,下面添加缺損的Xi的值,以獲得完全數(shù)據(jù)的似然函數(shù).具體如下:

當(dāng)αi=1,βi=0,即第i個產(chǎn)品失效未被顯示時,添加數(shù)據(jù)Z1i=Xi=z1i.

在Xi≤zi的條件下,Z1i的條件密度函數(shù)為:

它是區(qū)間(0,zi]上的截斷伽瑪分布Ga(b,θ),其樣本可表示為：

z1i=F-1{F(0)+U0[F(zi)-F(0)]}=F-1{U0F(zi)},

其中，F(xiàn)是Ga(b,θ)的分布函數(shù),F-1為其反函數(shù),U0為來自均勻分布U(0,1)的一個隨機樣本.

當(dāng)αi=0,即Xi>Yi時,添加數(shù)據(jù)Z2i=Xi=z2i.

Z2i也服從截斷伽瑪分布,其樣本的獲得與Z1i類似.

令α表示αi組成的向量,β表示βi組成的向量,z表示zi組成的向量,u表示z1i組成的向量,v表示z2i組成的向量,則完全數(shù)據(jù)似然函數(shù)為：

[f(z1i;θ)]αi(1-βi)[f(z2i;θ)]1-αi}=

2.2EM算法

EM算法是一種迭代方法,并主要用來求后驗分布的眾數(shù),即極大似然估計,它的每一次迭代由兩步組成:E步(求期望)和M步(極大化).EM算法非常適合處理不完全數(shù)據(jù),下面利用EM算法來求θ的MLE.

取θ的先驗分布為伽瑪分布Ga(c,d),c,d已知,即

π(θ)∝θc-1e-dθ,θ>0,c>0,d>0.

則θ的添加后驗分布為：

L(θ|z,u,v,α,β)∝π(θ)L(z,u,v,α,β|θ)∝

在第m+1次迭代中,假設(shè)有估計值θ(m),則可通過E步和M步得到θ的一個新的估計.

令U表示Z1i組成的向量,V表示Z2i組成的向量.

為了書寫方便,簡記(|θ(m)，z，α,β)為(|·).

E步:

當(dāng)αi=1,βi=0時,

當(dāng)αi=0時,

而要獲得Z1i,Z2i期望的顯式表示是不可能的,這時可用Monte Carlo方法完成,這就是所謂的Monte Carlo EM(MCEM)方法,它將E步改為:

(E1) 利用逆變換法隨機產(chǎn)生li個隨機數(shù)k1,k2,…,kli.利用逆變換法隨機產(chǎn)生ri個隨機數(shù)s1,s2,…,sri;

(E2) 計算

(2)

M步:

得

θ(m+1)=

(3)

(3)式給出了由EM算法得到的形狀參數(shù)θ的迭代公式.

由上述討論可知,EM算法只適用于nb+c-1>0的情況.

2.3參數(shù)的貝葉斯估計

當(dāng)αi=1,βi=0時,

π(z1i|θ,z,z-1i,v,α,β)∝ψ1(z1i;zi,θ),

其中,z-1i={z1j:j≠i}.

當(dāng)αi=0時,

π(z2i|θ,z,u,z-2i,α,β)∝ψ2(z2i;zi,θ),

其中,z-2i={z2j:j≠i}.

參數(shù)θ的滿條件分布為:

π(θ|b,z,u,v,α,β)∝

(4)

由于得到了各參數(shù)的滿條件分布,下面利用MCMC方法獲得參數(shù)θ后驗分布的平穩(wěn)分布.z1i,z2i的滿條件分布都是截斷伽瑪分布,可以利用逆變換法隨機抽樣;θ的滿條件分布是伽瑪分布,是一個標(biāo)準(zhǔn)分布,所以這3個分布都可以采用Gibbs抽樣.

設(shè)θj,j=1,2,…,M1,…,M2為參數(shù)θ的一個容量為M2的Gibbs樣本,其中M1為由于不穩(wěn)定而舍棄的樣本容量,然后利用達到平穩(wěn)狀態(tài)的M2-M1個獨立樣本的均值作為參數(shù)θ的貝葉斯估計,即

(5)

EM算法得到的估計是后驗分布的眾數(shù)(MLE),MCMC方法得到的估計是后驗分布的期望,由于θ的滿條件分布是伽瑪分布,并且滿條件分布抽樣依分布收斂到后驗分布,所以為了討論EM估計和MCMC估計的關(guān)系,下面先討論一下伽瑪分布的眾數(shù)和期望的關(guān)系.

若X～Ga(λ1,λ2),則E(X)=λ1/λ2.當(dāng)0<λ1≤1時,X的密度函數(shù)嚴(yán)格下降,眾數(shù)不存在;當(dāng)λ1>1時,X的密度函數(shù)是單峰函數(shù),眾數(shù)為(λ1-1)/λ2,此時眾數(shù)比期望小,但當(dāng)λ2較大時兩者的差別較小.

3隨機模擬

基于上面的討論,下面進行隨機模擬試驗.設(shè)Xi服從伽瑪分布Ga(11,θ),θ=9,取參數(shù)θ的先驗分布為伽瑪分布Ga(5.2,0.6);截尾變量Yi～Ga(13,8),失效顯示概率a=0.8,樣本容量分別取n=30,50,100,200,300,500,800,1000.對每一固定樣本量隨機產(chǎn)生一組不完全隨機截尾數(shù)據(jù),然后根據(jù)此組數(shù)據(jù)分別利用EM算法和MCMC方法對參數(shù)θ進行估計.

運用EM算法時從θ=16開始迭代;運用MCMC方法時,從θ=12.5開始進行Gibbs抽樣,MCMC模擬運行過程中,先進行10000次Gibbs預(yù)迭代,以確保參數(shù)的收斂性,然后丟棄最初的預(yù)迭代,再進行10000次Gibbs迭代,把第10001次至第20000次的迭代值的算術(shù)平均值作為λ的估計值. EM估計和MCMC估計的R程序運行結(jié)果參見表1.

表1 IIRCT下伽瑪分布參數(shù)估計的隨機模擬結(jié)果Tab.1 Random simulation results of parameter estimation of gamma distribution for IIRCT

當(dāng)n=300時,EM算法迭代10次,Gibbs迭代20000次,迭代過程參見圖1和圖2.

圖1 n=300時,參數(shù)θ的EM迭代過程Fig.1 EM iterations of θ with n=300

圖2 n=300時,參數(shù)θ的Gibbs抽樣迭代過程Fig.2 Gibbs sampling iterations of θ with n=300

由表1可以看出,θ的EM 估計和MCMC估計的差別不大,與真值9的相對誤差都小于7%,MCMC估計比EM估計更接近真實參數(shù)θ;樣本量對估計值的影響也不大,說明得到的估計值是比較穩(wěn)定的,并且精度也較高.

注1 在編寫R程序時用到的函數(shù)主要有:

gamma(), unif(), binom(), min(), max(), sum(), mean(), plot()等.

[1] Gupta S S. Order statistics from the gamma distribution[J]. Technometrics, 1960, 2(2): 243-262.

[2] Choi S C, Wette R. Maximum likelihood estimation of the parameters of the gamma distribution and their bias[J]. Technometrics, 1969, 11(4): 683-690.

[3] Wilks D S. Maximum likelihood estimation for the gamma distribution using data containing zeros[J]. Journal of Climate, 1990, 3(12): 1495-1501.

[4] Khodabin M, Ahmadabadi A. Some properties of generalized gamma distribution[J]. J Math Sci, 2010, 4: 9-28.

[5] Huang W, Shu L, Jiang W, et al. Evaluation of run-length distribution for CUSUM charts under gamma distributions[J]. IIE Transactions, 2013, 45(9): 981-994.

[6] Alam A, Rahman M S, Saadat A H M, et al. Gamma distribution and its application of spatially monitoring meteorological drought in barind, bangladesh[J]. Journal of Environmental Science and Natural Resources, 2013, 5(2): 287-293.

[7] Elperin T I, Gertsbakh I B. Estimation in a random censoring model with incomplete information: exponential lifetime distribution[J]. IEEE Transactions on Reliability, 1988, 37(2): 223-229.

[8] Elperin T, Gertsbakh I. Bayes credibility estimation of an exponential parameter for random censoring and incomplete information[J]. IEEE Transactions on Reliability, 1990, 39(2): 204-208.

[9] Ye E. Consistency of MLE of the parameter of exponential lifetime distribution for random censoring model with incomplete information[J]. Applied Mathematics, 1995, 10(4): 379-386.

[10] 陳怡南, 葉爾驊. IIRCT下對數(shù)正態(tài)和正態(tài)分布參數(shù)的MLE[J].南京航空航天大學(xué)學(xué)報:自然科學(xué)版, 1996,28(3):376-383.

[11] 陳怡南, 葉爾驊. 帶有不完全信息隨機截尾試驗下Weibull分布參數(shù)的MLE[J].數(shù)理統(tǒng)計與應(yīng)用概率, 1996, 11(4):353-363.

[12] 楊紀(jì)龍, 葉爾驊. 帶有不完全信息隨機截尾試驗下Weibull分布參數(shù)的MLE的相合性及漸近正態(tài)性[J].應(yīng)用概率統(tǒng)計, 2000, 16(1):9-19.

[13] 張曉琴, 張虎明. 帶有不完全信息隨機截尾試驗下Weibull分布參數(shù)的MLE的重對數(shù)律[J]. 應(yīng)用概率統(tǒng)計,2002,18(1):101-107.

[14] 宋毅君, 李補喜. 帶有不完全信息隨機截尾試驗下最大似然估計的相合性及漸近正態(tài)性[J].應(yīng)用概率統(tǒng)計, 2003, 19(2):139-149.

[15] 宋毅君, 李補喜, 李濟洪. 帶有不完全信息隨機截尾試驗下最大似然估計的重對數(shù)律[J]. 應(yīng)用概率統(tǒng)計, 2009, 25(2):113-125.

[16] 宋鳳麗, 沈 思. 帶有不完全信息隨機截尾試驗下總體均值型參數(shù)的經(jīng)驗似然[J].應(yīng)用數(shù)學(xué), 2009, 22(1):191-198.

[17] 劉有新, 戴 揚. 帶有不完全信息隨機截尾數(shù)據(jù)的Fisher信息陣[J].重慶工商大學(xué)學(xué)報:自然科學(xué)版, 2008, 25(6):569-572.

[18] 何朝兵. 帶有不完全信息隨機截尾試驗下幾何分布的參數(shù)估計[J].應(yīng)用數(shù)學(xué), 2013, 26(3):574-579.

Point estimation of scale parameter of Gamma distribution for random censoring test model with incomplete information

HE Chaobing1, LIU Huawen2

(1.School of Mathematics and Statistics, Anyang Normal University, Anyang, Henan 455000;2.School of Mathematics, Shandong University, Jinan 250100)

In this paper we firstly obtain the complete data likelihood function of gamma distribution for IIRCT after adding data， then estimate the scale parameter by EM algorithm and MCMC method， respectively. Finally random simulation tests are conducted， and the results show that the point estimations of the scale parameter are fairly accurate．

complete data likelihood function; EM algorithm; full conditional distribution; MCMC method; Gibbs sampling

2014-01-13.

國家自然科學(xué)基金項目(61174099); 河南省教育廳自然科學(xué)基金項目(2011B110001).

1000-1190(2014)04-0479-04

O213.2; O212.8

A

*E-mail: chaobing5@163.com.