刪失截?cái)嗲樾蜗率首凕c(diǎn)模型的Bayes參數(shù)估計(jì)

何朝兵

(安陽(yáng)師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,河南安陽(yáng)455000)

刪失截?cái)嗲樾蜗率首凕c(diǎn)模型的Bayes參數(shù)估計(jì)

何朝兵

(安陽(yáng)師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,河南安陽(yáng)455000)

通過(guò)添加部分缺失壽命變量數(shù)據(jù),得到了刪失截?cái)嗲樾蜗率首凕c(diǎn)模型相對(duì)簡(jiǎn)單的似然函數(shù).討論了所添加缺失數(shù)據(jù)變量的概率分布和隨機(jī)抽樣方法.利用Monte Carlo EM算法對(duì)未知參數(shù)進(jìn)行了迭代.結(jié)合Metropolis-Hastings算法對(duì)參數(shù)的滿條件分布進(jìn)行了Gibbs抽樣,基于Gibbs樣本對(duì)參數(shù)進(jìn)行估計(jì),詳細(xì)介紹了MCMC方法的實(shí)施步驟.隨機(jī)模擬試驗(yàn)的結(jié)果表明各參數(shù)Bayes估計(jì)的精度較高.

失效率;指數(shù)分布;EM算法;Gibbs抽樣;Metropolis-Hastings算法;截?cái)嗾龖B(tài)分布

§1 引 言

在可靠性統(tǒng)計(jì)中,已工作到某時(shí)刻尚未失效的產(chǎn)品,在這個(gè)時(shí)刻后單位時(shí)間內(nèi)失效的概率稱為該產(chǎn)品(壽命)在這個(gè)時(shí)刻的失效率函數(shù),簡(jiǎn)稱失效率.由于失效率函數(shù)能在給定的時(shí)刻對(duì)瞬間的失效風(fēng)險(xiǎn)進(jìn)行量化,所以失效率函數(shù)在可靠性和生存分析的研究中起著重要的作用.有時(shí)失效率函數(shù)在某個(gè)時(shí)刻會(huì)突然變化,這個(gè)時(shí)刻稱為變點(diǎn),例如特別的治療使病人的身體健康狀況在這個(gè)時(shí)刻有顯著改善,疲勞導(dǎo)致機(jī)器設(shè)備的物理?xiàng)l件在這個(gè)時(shí)刻急劇下降等.對(duì)具有變點(diǎn)的失效率函數(shù)的估計(jì)是一項(xiàng)有趣和富有挑戰(zhàn)性的任務(wù).文獻(xiàn)[1]最早研究了失效率變點(diǎn)模型,得到了參數(shù)的極大似然估計(jì),接著文獻(xiàn)[2]研究了失效率變點(diǎn)的近似置信區(qū)間,文獻(xiàn)[3]研究了自助法的漸近有效性.關(guān)于失效率變點(diǎn)模型的更多研究可參看文獻(xiàn)[4-10].

上述的文獻(xiàn)大多數(shù)都是基于完全數(shù)據(jù)或右刪失數(shù)據(jù)進(jìn)行統(tǒng)計(jì)推斷的.進(jìn)行壽命試驗(yàn)時(shí),會(huì)經(jīng)常出現(xiàn)截?cái)鄶?shù)據(jù),截?cái)嘤址Q為左截?cái)?個(gè)體由于某種原因中途撤離了試驗(yàn),或者被收集的個(gè)體信息中途失去或統(tǒng)計(jì)數(shù)據(jù)時(shí)個(gè)體尚未壽終,這就是右刪失.如果個(gè)體在研究開(kāi)始之前就已經(jīng)壽終或個(gè)體本身因條件限制根本無(wú)法觀察,這樣所獲得的個(gè)體數(shù)據(jù)就是左截?cái)鄶?shù)據(jù).有時(shí)截?cái)嗪蛣h失會(huì)出現(xiàn)在同一個(gè)試驗(yàn)中,例如,對(duì)前期已經(jīng)確診為艾滋病的患者進(jìn)行研究,感興趣的變量是患者從確診到死亡的時(shí)間,如果個(gè)體在研究的最后提前退出研究或仍然存活,這就是右刪失,如果個(gè)體在研究還未開(kāi)始之前就已經(jīng)死亡,結(jié)果導(dǎo)致無(wú)法觀察,這就是左截?cái)?刪失截?cái)鄶?shù)據(jù)模型現(xiàn)已廣泛應(yīng)用于生物學(xué)、醫(yī)學(xué)、人口學(xué)、經(jīng)濟(jì)學(xué)等方面,對(duì)此模型的研究可參看文獻(xiàn)[11-14],其中文[11]研究了此模型的回歸分析和非參數(shù)估計(jì),文[12]研究了此模型下對(duì)數(shù)正態(tài)分布的最大然估計(jì).由于失效率變點(diǎn)模型和刪失截?cái)鄶?shù)據(jù)模型的廣泛應(yīng)用,所以對(duì)刪失截?cái)嗲樾蜗率首凕c(diǎn)模型的研究是一個(gè)有意義的研究方向,具有重要的理論價(jià)值和實(shí)用價(jià)值.

Bayes計(jì)算方法中的數(shù)據(jù)添加算法[15]是近年發(fā)展很快且應(yīng)用很廣的一種算法,它是在觀測(cè)數(shù)據(jù)的基礎(chǔ)上加上一些“潛在數(shù)據(jù)”,從而簡(jiǎn)化計(jì)算并完成一系列簡(jiǎn)單的極大化或隨機(jī)模擬,該“潛在數(shù)據(jù)”可以是“缺損數(shù)據(jù)”或未知參數(shù).EM算法和Markov chain Monte Carlo(MCMC)方法是兩種最常用的數(shù)據(jù)添加算法.EM算法是一種求后驗(yàn)分布的眾數(shù)(即極大似然估計(jì))迭代方法,它的每一步迭代由E步(求期望)和M步(極大化)組成.MCMC方法主要包括Gibbs抽樣和Metropolis-Hastings算法,通常是Gibbs抽樣和Metropolis-Hastings算法相結(jié)合得到Gibbs樣本,然后基于此樣本對(duì)后驗(yàn)分布的各種統(tǒng)計(jì)量進(jìn)行估計(jì).EM算法和MCMC方法使貝葉斯統(tǒng)計(jì)中許多看起來(lái)困難的計(jì)算變得簡(jiǎn)單直觀,其應(yīng)用很廣泛,比如約束參數(shù)模型,變點(diǎn)模型,截尾數(shù)據(jù)和分組數(shù)據(jù),多個(gè)分布的混合等.隨著統(tǒng)計(jì)計(jì)算技術(shù)的發(fā)展,Bayes計(jì)算方法越來(lái)越多地應(yīng)用到變點(diǎn)分析之中.文獻(xiàn)[16-18]利用EM算法對(duì)變點(diǎn)模型進(jìn)行了統(tǒng)計(jì)推斷,文獻(xiàn)[19-22]利用MCMC方法研究了變點(diǎn)模型.

由于刪失截?cái)鄶?shù)據(jù)下的失效率變點(diǎn)模型是截尾數(shù)據(jù)模型和變點(diǎn)模型的混合,所以貝葉斯計(jì)算方法中的EM算法和MCMC方法很適合研究此變點(diǎn)模型.關(guān)于截?cái)嗷騽h失數(shù)據(jù)下變點(diǎn)問(wèn)題的研究已有一些成果,可參看文獻(xiàn)[23-26],其中文[23]和[24]分別研究了刪失數(shù)據(jù)下失效率變點(diǎn)問(wèn)題的非參數(shù)和參數(shù)估計(jì).對(duì)刪失并且截?cái)鄶?shù)據(jù)下失效率變點(diǎn)模型的研究文只有文[27],但它利用的方法還是傳統(tǒng)的最大似然法,它主要依靠一般的數(shù)值計(jì)算方法分區(qū)間逐段尋找變點(diǎn)參數(shù)的最大值點(diǎn),進(jìn)而得到參數(shù)的最大似然估計(jì).而利用EM算法和MCMC方法對(duì)刪失截?cái)嗲樾蜗率首凕c(diǎn)模型的研究還未見(jiàn)到.

下文主要利用EM算法和MCMC方法研究了刪失截?cái)嗲樾蜗率首凕c(diǎn)模型的參數(shù)估計(jì)問(wèn)題.EM算法中的E步是通過(guò)Monte Carlo方法完成的.結(jié)合Metropolis-Hastings算法對(duì)各未知參數(shù)的滿條件分布進(jìn)行了Gibbs抽樣,基于Gibbs樣本對(duì)參數(shù)進(jìn)行估計(jì).隨機(jī)模擬試驗(yàn)的結(jié)果表明各參數(shù)Bayes估計(jì)的精度較高.

§2 刪失截?cái)鄶?shù)據(jù)試驗(yàn)?zāi)Ｐ?/h2>

設(shè)(X,Y,T)是一連續(xù)型隨機(jī)變量,壽命變量X的分布函數(shù)為F(x|θ)=P(X ≤x),密度函數(shù)為f(x|θ),這里θ是未知參數(shù)向量;Y是一右刪失隨機(jī)變量,分布函數(shù)為G(y),密度函數(shù)為g(y);T是一左截?cái)嚯S機(jī)變量,分布函數(shù)為H(t),密度函數(shù)為h(t),且Y,T的分布與參數(shù)θ無(wú)關(guān).假定X,Y,T是相互獨(dú)立非負(fù)隨機(jī)變量.對(duì)于n個(gè)受試樣品(產(chǎn)品壽命),刪失截?cái)鄶?shù)據(jù)試驗(yàn)?zāi)Ｐ褪?僅在Zi≥Ti時(shí)得到觀察數(shù)據(jù),而在Zi＜Ti下無(wú)法得到任何觀察值,其中

當(dāng)有樣本觀察值時(shí),下面求(Zi,Ti,δi)的密度函數(shù).

當(dāng)Zi=zi,Ti=ti,δi=1時(shí),(Zi,Ti,δi)的密度函數(shù)為f(zi|θ)(zi)h(ti),zi≥ti;

當(dāng)Zi=zi,Ti=ti,δi=0時(shí),(Zi,Ti,δi)的密度函數(shù)為g(zi)(zi|θ)h(ti),zi≥ti.

P(無(wú)樣本觀察值)=P(Zi＜Ti)=1-P(Xi≥Ti,Yi≥Ti)

為了研究方便,引入示性變量νi=I(min(Xi,Yi)≥Ti),i=1,2,···,n.

則基于觀察數(shù)據(jù){(zi,ti,δi)：νi=1,1≤i≤n}的似然函數(shù)為

缺損數(shù)據(jù)下的似然函數(shù)比較復(fù)雜,參數(shù)估計(jì)一般很難處理.而B(niǎo)ayes統(tǒng)計(jì)中的MCMC算法處理缺損數(shù)據(jù)非常方便,其方法步驟如下.先引進(jìn)添加數(shù)據(jù),求出添加數(shù)據(jù)的概率分布,把添加數(shù)據(jù)作為未知參數(shù)處理,然后考慮添加數(shù)據(jù)和未知參數(shù)的后驗(yàn)分布,未知參數(shù)的滿條件分布可轉(zhuǎn)化為較簡(jiǎn)單的后驗(yàn)分布,而添加數(shù)據(jù)的滿條件分布則是由添加數(shù)據(jù)的概率分布抽樣,最后分別對(duì)各滿條件分布進(jìn)行Gibbs抽樣即可.

下面添加缺損的部分Xi的值,以獲得較簡(jiǎn)單的似然函數(shù).

當(dāng)νi=0時(shí),添加數(shù)據(jù)Xi=Z1i=z1i.

在min(Xi,Yi)＜Ti的條件下,Z1i的條件密度函數(shù)為

下面利用兩種方法隨機(jī)產(chǎn)生z1i.第一種方法稱為舍選法,步驟如下：

若B(x)計(jì)算煩瑣并且容易產(chǎn)生來(lái)自g(y),h(t)的隨機(jī)數(shù),可以利用第二種方法產(chǎn)生z1i：

(1)分別產(chǎn)生來(lái)自f(x|θ),g(y),h(t)的隨機(jī)數(shù)x,y,t;

(2)如果min(x,y)＜t,令z1i=x,停止;

(3)如果min(x,y)≥t,回到步驟(1).

則添加數(shù)據(jù)后的似然函數(shù)為

§3 失效率變點(diǎn)模型

設(shè)非負(fù)連續(xù)型隨機(jī)變量X的分布函數(shù)和密度函數(shù)分別為F(x|θ)和f(x|θ),稱λ(x|θ)=f(x|θ)/[1-F(x|θ)]為X的失效率函數(shù). 若λ(x)滿足

(1)對(duì)于λi取共軛先驗(yàn)分布Gamma分布Ga(bi,ci),bi,ci已知,即

(2)對(duì)于τ取無(wú)信息先驗(yàn)分布：π(τ)∝1,τ＞0.

假設(shè)λ1,λ2,τ相互獨(dú)立,則

下面討論截?cái)嘧兞亢蛣h失變量都服從指數(shù)分布這種特殊情況.

假設(shè)Y ～ Exp(λy),T ～ Exp(λt),則

下面驗(yàn)證?(x|θ)作為密度函數(shù)的規(guī)范性.

則Z1i的分布是2個(gè)失效率變點(diǎn)模型壽命分布的混合,這個(gè)發(fā)現(xiàn)很有趣.根據(jù)(2)式,還可以利用合成法結(jié)合逆變換法直接產(chǎn)生z1i,步驟如下：

(1)產(chǎn)生U(0,1)隨機(jī)數(shù)u1,u2,u3.

(2)計(jì)算

§4 EM算法

EM算法非常適合處理不完全數(shù)據(jù),下面利用EM算法來(lái)求參數(shù)后驗(yàn)分布的眾數(shù).假設(shè)在第m+1次迭代開(kāi)始時(shí)有估計(jì)值則可通過(guò)E步和M步得到的一個(gè)新的估計(jì).令W表示Z1i組成的向量,為了書(shū)寫方便,簡(jiǎn)記(|θ(m),z,δ,ν)為(|·).E步：

由于?(x|θ(m))為分段函數(shù),積分比較繁瑣,大部分情況下要獲得期望的顯式表示是不可能的,這時(shí)可用Monte Carlo方法完成,這就是所謂的Monte Carlo EM(MCEM)方法.

在βi=1的條件下,Z1i的條件密度函數(shù)為

§5 MCMC方法

下面對(duì)各參數(shù)進(jìn)行利用Gibbs抽樣.首先研究各參數(shù)的滿條件分布.

利用前面介紹的舍選法等方法產(chǎn)生z1i,由于λ1,λ2的滿條件分布是Gamma分布,所以z1i,λ1,λ2都可以直接Gibbs抽樣;但τ的滿條件分布比較復(fù)雜,不能直接進(jìn)行Gibbs抽樣,利用Metropolis-Hastings算法進(jìn)行抽樣,選取區(qū)間(0,+∞)上的截?cái)嗾龖B(tài)分布作為建議分布.

下面介紹一下截?cái)喾植嫉某闃?

設(shè)X是連續(xù)型隨機(jī)變量,其密度函數(shù)及分布函數(shù)分別是f(x),F(x),x∈(-∞,+∞).令

則g(x)是由f(x)定義的在(a,b)上的截?cái)嗟姆植?設(shè)F-1(x)為F(x)的反函數(shù),U為由U(0,1)產(chǎn)生的隨機(jī)數(shù),則F-1{F(a)+U[F(b)-F(a)]}為從g(x)抽取的樣本值.

如果X ～ N(μ,σ2),有F(0)=P(X ≤ 0)= Φ(-μ/σ)=1- Φ(μ/σ),其中Φ為N(0,1)的分布函數(shù).則(0,∞)上截?cái)嗾龖B(tài)分布的密度函數(shù)為

下面介紹MCMC方法的具體步驟.

§6 參數(shù)的最大似然估計(jì)

為了把參數(shù)的Bayes估計(jì)與最大似然估計(jì)(MLE)進(jìn)行比較,下面簡(jiǎn)單介紹MLE的求解步驟.

先考慮左截?cái)嘤覄h失模型.為了書(shū)寫方便,不妨假定后n0個(gè)壽命變量沒(méi)有觀察值,前n1=n-n0個(gè)壽命變量有觀察值z(mì)1＜z2＜···＜zn1,則似然函數(shù)為

再考慮失效率變點(diǎn)模型,假設(shè)τ∈(zk,zk+1],k=0,1,···,n1,其中z0=0,zn1+1=+∞,

§7 隨機(jī)模擬

下面利用R軟件進(jìn)行隨機(jī)模擬試驗(yàn).

隨機(jī)模擬時(shí)取受試樣品的個(gè)數(shù)n=200,θ =(λ1,λ2,τ)=(2,5,0.3),取λ1,λ2的先驗(yàn)分布分別為Ga(0.4,0.15),Ga(2.7,0.5),右刪失變量Yi～Exp(1.5),左截?cái)嘧兞縏i～Exp(5.5).實(shí)際隨機(jī)模擬數(shù)據(jù)中,沒(méi)有觀察值的樣本個(gè)數(shù)n0=82,有觀察值的樣本個(gè)數(shù)n1=118.各參數(shù)的Bayes估計(jì)和MLE見(jiàn)表1.

對(duì)參數(shù)進(jìn)行EM迭代時(shí)選取的初始值為(λ1,λ2,τ)=(5,3,0.8),參數(shù)的EM迭代過(guò)程見(jiàn)圖1至圖3.由于在E步是采用的Monte Carlo方法,所以當(dāng)?shù)祰@一個(gè)定值小幅波動(dòng)時(shí),則可以認(rèn)為算法收斂了,此時(shí)為了增加估計(jì)精度,可增加Monte Carlo抽樣的數(shù)目,再運(yùn)行迭代一段時(shí)間即可停止,停止后,我們把后10次迭代值的均值作為參數(shù)的估計(jì)值.

對(duì)參數(shù)進(jìn)行Gibbs迭代時(shí)選取的初始值為(λ1,λ2,τ)=(8,4,2),取B=10000,M=20000,參數(shù)的Gibbs迭代過(guò)程見(jiàn)圖4至圖6.對(duì)MCMC的收斂性診斷很重要,模擬時(shí)可以對(duì)參數(shù)進(jìn)行多層鏈?zhǔn)降治?即輸入多組初始值,形成多層迭代鏈,當(dāng)抽樣收斂時(shí),迭代圖形重合.在模擬過(guò)程中,輸入兩組初始值分別進(jìn)行10000次迭代.各參數(shù)的兩條迭代鏈軌跡見(jiàn)圖7至圖9.另外,為了和EM估計(jì)(后驗(yàn)眾數(shù)估計(jì))相比較,給出基于后M-B=10000個(gè)Gibbs樣本的核密度眾數(shù)估計(jì),各參數(shù)的核密度估計(jì)見(jiàn)圖10至圖12.

求解參數(shù)的MLE時(shí),對(duì)τ所在的n1+1=119個(gè)區(qū)間分別建立對(duì)數(shù)似然方程組,得到119個(gè)解及其對(duì)應(yīng)的對(duì)數(shù)似然函數(shù)值,最大的對(duì)數(shù)似然函數(shù)值對(duì)應(yīng)的解即為要求的MLE.這119個(gè)解的τ分量及其對(duì)應(yīng)的對(duì)數(shù)似然函數(shù)值的對(duì)照?qǐng)D見(jiàn)圖13.

圖1 λ1的EM迭代過(guò)程

圖2 λ2的EM迭代過(guò)程

圖3 τ的EM迭代過(guò)程

圖4 λ1的Gibbs抽樣迭代過(guò)程

下面進(jìn)行統(tǒng)計(jì)分析：

1)由表1可看出,后驗(yàn)分布的眾數(shù)估計(jì)(即EM估計(jì))、均值估計(jì)、中位數(shù)估計(jì)、核密度眾數(shù)估計(jì)以及MLE的差別很小,λ1,τ估計(jì)值與真值的的誤差較小,相對(duì)誤差小于4%;λ2的誤差稍大,但相對(duì)誤差也沒(méi)超過(guò)11%,所以整體上,各參數(shù)估計(jì)的精度較高.EM估計(jì)與后驗(yàn)均值估計(jì)、MLE相比誤差較小.Gibbs抽樣的MC誤差也較小,估計(jì)效果較好.各參數(shù)可信水平為0.95的可信區(qū)間可近似取為[2.5%分位數(shù),97.5%分位數(shù)],可以看出近似可信區(qū)間的長(zhǎng)度非常短,所以區(qū)間估計(jì)的效果也較好.

2)由圖1至圖3可看出,EM迭代10次左右就收斂了,但由于是利用搜索的方法尋找τ的極大值點(diǎn),計(jì)算機(jī)的迭代速度實(shí)際上并不快.

3)由圖4至圖6可看出,Gibbs抽樣迭代值波動(dòng)較小,估計(jì)效果較好,雖然Gibbs抽樣迭代了20000次,也比EM迭代50次要快得多,這是因?yàn)镚ibbs抽樣迭代完全是利用隨機(jī)抽樣的方法,不必搜索最優(yōu)值,所以從實(shí)際操作角度來(lái)說(shuō),本文介紹的Gibbs抽樣要比EM算法要好.

4)由圖7至圖9可以發(fā)現(xiàn),各參數(shù)的兩條迭代鏈都趨于重合,收斂性較好.

5)由圖10至圖12可看出,參數(shù)核密度的最大值點(diǎn)與對(duì)應(yīng)參數(shù)的真值很接近,核密度眾數(shù)作為后驗(yàn)眾數(shù)的估計(jì)效果較好.

6)由圖13可看出,對(duì)數(shù)似然函數(shù)最大值點(diǎn)與參數(shù)τ的真值很接近,MLE的精度較高,但對(duì)數(shù)似然函數(shù)方程組比較復(fù)雜,所以分區(qū)間求解對(duì)數(shù)似然函數(shù)最大值點(diǎn)時(shí)非常繁瑣,這也是最大似然法的最大缺點(diǎn).

§8 結(jié) 論

圖5 λ2的Gibbs抽樣迭代過(guò)程

圖6 τ的Gibbs抽樣迭代過(guò)程

圖7 λ1的兩條迭代鏈軌跡

圖8 λ2的兩條迭代鏈軌跡

圖9 τ的兩條迭代鏈軌跡

圖10 λ1的后驗(yàn)核密度估計(jì)

圖11 λ2的后驗(yàn)核密度估計(jì)

圖12 τ的后驗(yàn)核密度估計(jì)

本文主要利用Bayes方法研究了刪失截?cái)嗲樾蜗率首凕c(diǎn)模型的參數(shù)估計(jì)問(wèn)題.添加了部分缺失壽命變量數(shù)據(jù)并且討論了它的概率分布和隨機(jī)抽樣方法.利用EM算法和MCMC方法得到了變點(diǎn)位置參數(shù)和其它參數(shù)的Bayes估計(jì),隨機(jī)模擬數(shù)據(jù)表明參數(shù)估計(jì)的精度較高.

本文研究的模型是截尾數(shù)據(jù)模型和變點(diǎn)模型的混合,而EM算法和MCMC方法非常適合處理截尾數(shù)據(jù)和變點(diǎn)問(wèn)題,這也是本文的優(yōu)點(diǎn).但本文利用的MCEM算法也有缺點(diǎn),就是利用優(yōu)化函數(shù)對(duì)τ進(jìn)行極大值搜索時(shí)速度較慢;雖然本文運(yùn)用的MCMC的收斂速度比MCEM算法快,但整體上還是慢的,有時(shí)MCMC是不收斂的,但還沒(méi)有有效的改進(jìn)方法.本文研究變點(diǎn)模型時(shí),假定有1個(gè)變點(diǎn).而如何確定變點(diǎn)的個(gè)數(shù),是尚待解決的問(wèn)題.實(shí)際上,變點(diǎn)個(gè)數(shù)的確定是變點(diǎn)問(wèn)題中的難題,也是變點(diǎn)研究的發(fā)展趨勢(shì),更有研究?jī)r(jià)值和現(xiàn)實(shí)意義.

圖13 τ與對(duì)數(shù)似然函數(shù)區(qū)間最大值的對(duì)照?qǐng)D

表1 參數(shù)λ1,λ2,τ的Bayes估計(jì)和MLE

[1] Matthews D E,Farewell V T.On testing for a constant hazard against a change-point alternative[J].Biometrics,1982,38(2)：463-468.

[2] Loader C R.Inference for a hazard rate change point[J].Biometrika,1991,78(4)：749-757.

[3] Pham T D,Nguyen H T.Bootstrapping the change-point of a hazard rate[J].Annals of the Institute of Statistical Mathematics,1993,45(2)：331-340.

[4] Antoniadis A,Gijbels I,Macgibbon B.Non-parametric estimation for the location of a change-point in an otherwise smooth hazard function under random censoring[J].Scandinavian Journal of Statistics,2000,27(3)：501-519.

[5] Dupuy J F.Detecting change in a hazard regression model with right-censoring[J].Journal of Statistical Planning and Inference,2009,139(5)：1578-1586.

[6] Frobish D,Ebrahimi N.Parametric estimation of change-points for actual event data in recurrent events models[J].Computational Statistics&Data Analysis,2009,53(3)：671-682.

[7]Gijbels I,Gürlerü.Estimation of a change point in a hazard function based on censored data[J].Lifetime Data Analysis,2003,9(4)：395-411.

[8] Karasoy D S,Kadilar C.A new Bayes estimate of the change point in the hazard function[J].Computational statistics&data analysis,2007,51(6)：2993-3001.

[9] Liu Mengling,Lu Wenbin,Shao Yongzhao.A Monte Carlo approach for change-point detection in the Cox proportional hazards model[J].Statistics in Medicine,2008,27(19)：3894-3909.

[10]Luo Xiaolong,Turnbull B W,Clark L C.Likelihood ratio tests for a changepoint with survival data[J].Biometrika,1997,84(3)：555-565.

[11]Gross S T,Lai T L.Nonparametric estimation and regression analysis with left-truncated and right-censored data[J].Journal of the American Statistical Association,1996,91(435)：1166-1180.

[12]Balakrishnan N,Mitra D.Likelihood inference for lognormal data with left truncation and right censoring with an illustration[J].Journal of Statistical Planning and Inference,2011,141(11)：3536-3553.

[13]Cosslett S R.Efficient semiparametric estimation of censored and truncated regressions via a smoothed self-consistency equation[J].Econometrica,2004,72(4)：1277-1293.

[14]Molanes-lopez E M,Cao R,Keilegom I VAN.Smoothed empirical likelihood con fi dence intervals for the relative distribution with left-truncated and right-censored data[J].Canadian Journal of Statistics,2010,38(3)：453-473.

[15] 茆詩(shī)松,王靜龍,濮曉龍.高等數(shù)理統(tǒng)計(jì)(第二版)[M].北京：高等教育出版社,2006.

[16]Yildirim S,Singh S S,Doucet A.An online expectation-maximization algorithm for changepoint models[J].Journal of Computational and Graphical Statistics,2013,22(4)：906-926.

[17]Comert G,Bezuglov A.An online change-point-based model for traffic parameter prediction[J].IEEE Transactions on Intelligent Transportation Systems,2013,14(3)：1360-1369.

[18]Bansal N K,Du H,Hamedani G G.An application of EM algorithm to a change-point problem[J].Communications in Statistics-Theory and Methods,2008,37(13)：2010-2021.

[19]Lavielle M,Lebarbier E.An application of MCMC methods for the multiple change-points problem[J].Signal Processing,2001,81(1)：39-53.

[20]Kim J,Cheon S.Bayesian multiple change-point estimation with annealing stochastic approximation Monte Carlo[J].Computational Statistics,2010,25(2)：215-239.

[21]Yuan Tao,Kuo Yue.Bayesian analysis of hazard rate,change point,and cost-optimal burnin time for electronic devices[J].IEEE Transactions on Reliability 2010,59(1)：132-138.

[22]Liang Faming,Wong W H.Real-parameter evolutionary Monte Carlo with applications to Bayesian mixture models[J].Journal of the American Statistical Association,2001,96(454)：653-666.

[23]Antoniadis A,Grégoire G.Nonparametric estimation in change point hazard rate models for censored data：A counting process approach[J].Journal of Nonparametric Statistics,1993,3(2)：135-154.

[24]Gijbels I,Gürlerü.Estimation of a change point in a hazard function based on censored data[J].Lifetime Data Analysis,2003,9(4)：395-411.

[25]Hu?kov′a M,Neuhaus G.Change point analysis for censored data[J].Journal of Statistical Planning and Inference,2004,126(1)：207-223.

[26]Lim H,Sun Jianguo,Matthews D E.Maximum likelihood estimation of a survival function with a change point for truncated and interval-censored data[J].Statistics in Medicine,2002,21(5)：743-752.

[27]Gürlerü,Yenigün C D.Full and conditional likelihood approaches for hazard change-point estimation with truncated and censored data[J].Computational Statistics&Data Analysis,2011,55(10)：2856-2870.

Bayesian parameter estimation of failure rate model with a change point for truncated and censored data

HE Chao-bing
(School of Mathematics and Statistics,Anyang Normal University,Anyang 455000,China)

By fi lling in some missing data of the life variable,the relatively simple likelihood function of failure rate model with a change point for truncated and censored data is obtained.The probability distribution and random sampling method of the missing data variable fi llled in are discussed.All the unknown parameters are iterated by MCEM algorithm.The parameters are sampled from their full conditional distributions by Gibbs sampler together with Metropolis-Hastings algorithm,and are estimated based on Gibbs sample.The implementation steps of MCMC method are introduced in detail.The random simulation test results show that Bayesian estimations of the parameters are fairly accurate.

failure rate;exponential distribution;EM algorithm;Gibbs sampling;Metropolis-Hastings algorithm;truncated normal distribution

62F15;62N01

O213.2;O212.8

A

：1000-4424(2016)04-0413-15

2016-01-19

2016-07-03

河南省科技攻關(guān)計(jì)劃(162102310384);河南省高等學(xué)校重點(diǎn)科研項(xiàng)目(16A110001)