何朝兵, 劉華文

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

?

帶有不完全信息隨機(jī)截尾試驗(yàn)下伽瑪分布尺度參數(shù)的點(diǎn)估計(jì)

何朝兵1*, 劉華文2

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

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

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

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

1帶有不完全信息隨機(jī)截尾試驗(yàn)(IIRCT)

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

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

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

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

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

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

Xi≤Yi,且未被顯示時(shí),β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ù)的點(diǎn)估計(jì)

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

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

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

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

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

它是區(qū)間(0,zi]上的截?cái)噘が敺植糋a(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為來(lái)自均勻分布U(0,1)的一個(gè)隨機(jī)樣本.

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

Z2i也服從截?cái)噘が敺植?其樣本的獲得與Z1i類(lèi)似.

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

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

2.2EM算法

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

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

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

則θ的添加后驗(yàn)分布為：

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

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

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

為了書(shū)寫(xiě)方便,簡(jiǎn)記(|θ(m)，z，α,β)為(|·).

E步:

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

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

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

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

(E2) 計(jì)算

(2)

M步:

得

θ(m+1)=

(3)

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

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

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

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

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

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

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

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

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

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

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

(4)

由于得到了各參數(shù)的滿(mǎn)條件分布,下面利用MCMC方法獲得參數(shù)θ后驗(yàn)分布的平穩(wěn)分布.z1i,z2i的滿(mǎn)條件分布都是截?cái)噘が敺植?可以利用逆變換法隨機(jī)抽樣;θ的滿(mǎn)條件分布是伽瑪分布,是一個(gè)標(biāo)準(zhǔn)分布,所以這3個(gè)分布都可以采用Gibbs抽樣.

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

(5)

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

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

3隨機(jī)模擬

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

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

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

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

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

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

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

注1 在編寫(xiě)R程序時(shí)用到的函數(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下對(duì)數(shù)正態(tài)和正態(tài)分布參數(shù)的MLE[J].南京航空航天大學(xué)學(xué)報(bào):自然科學(xué)版, 1996,28(3):376-383.

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

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

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

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

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

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

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

[18] 何朝兵. 帶有不完全信息隨機(jī)截尾試驗(yàn)下幾何分布的參數(shù)估計(jì)[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.

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

1000-1190(2014)04-0479-04

O213.2; O212.8

A

*E-mail: chaobing5@163.com.