何朝兵

(安陽師范學院數(shù)學與統(tǒng)計學院,河南安陽455000)

刪失截斷情形下失效率變點模型的Bayes參數(shù)估計

何朝兵

(安陽師范學院數(shù)學與統(tǒng)計學院,河南安陽455000)

通過添加部分缺失壽命變量數(shù)據(jù),得到了刪失截斷情形下失效率變點模型相對簡單的似然函數(shù).討論了所添加缺失數(shù)據(jù)變量的概率分布和隨機抽樣方法.利用Monte Carlo EM算法對未知參數(shù)進行了迭代.結合Metropolis-Hastings算法對參數(shù)的滿條件分布進行了Gibbs抽樣,基于Gibbs樣本對參數(shù)進行估計,詳細介紹了MCMC方法的實施步驟.隨機模擬試驗的結果表明各參數(shù)Bayes估計的精度較高.

失效率;指數(shù)分布;EM算法;Gibbs抽樣;Metropolis-Hastings算法;截斷正態(tài)分布

§1 引 言

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

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

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

由于刪失截斷數(shù)據(jù)下的失效率變點模型是截尾數(shù)據(jù)模型和變點模型的混合,所以貝葉斯計算方法中的EM算法和MCMC方法很適合研究此變點模型.關于截斷或刪失數(shù)據(jù)下變點問題的研究已有一些成果,可參看文獻[23-26],其中文[23]和[24]分別研究了刪失數(shù)據(jù)下失效率變點問題的非參數(shù)和參數(shù)估計.對刪失并且截斷數(shù)據(jù)下失效率變點模型的研究文只有文[27],但它利用的方法還是傳統(tǒng)的最大似然法,它主要依靠一般的數(shù)值計算方法分區(qū)間逐段尋找變點參數(shù)的最大值點,進而得到參數(shù)的最大似然估計.而利用EM算法和MCMC方法對刪失截斷情形下失效率變點模型的研究還未見到.

下文主要利用EM算法和MCMC方法研究了刪失截斷情形下失效率變點模型的參數(shù)估計問題.EM算法中的E步是通過Monte Carlo方法完成的.結合Metropolis-Hastings算法對各未知參數(shù)的滿條件分布進行了Gibbs抽樣,基于Gibbs樣本對參數(shù)進行估計.隨機模擬試驗的結果表明各參數(shù)Bayes估計的精度較高.

§2 刪失截斷數(shù)據(jù)試驗模型

設(X,Y,T)是一連續(xù)型隨機變量,壽命變量X的分布函數(shù)為F(x|θ)=P(X ≤x),密度函數(shù)為f(x|θ),這里θ是未知參數(shù)向量;Y是一右刪失隨機變量,分布函數(shù)為G(y),密度函數(shù)為g(y);T是一左截斷隨機變量,分布函數(shù)為H(t),密度函數(shù)為h(t),且Y,T的分布與參數(shù)θ無關.假定X,Y,T是相互獨立非負隨機變量.對于n個受試樣品(產品壽命),刪失截斷數(shù)據(jù)試驗模型是,僅在Zi≥Ti時得到觀察數(shù)據(jù),而在Zi＜Ti下無法得到任何觀察值,其中

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

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

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

P(無樣本觀察值)=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ù)比較復雜,參數(shù)估計一般很難處理.而Bayes統(tǒng)計中的MCMC算法處理缺損數(shù)據(jù)非常方便,其方法步驟如下.先引進添加數(shù)據(jù),求出添加數(shù)據(jù)的概率分布,把添加數(shù)據(jù)作為未知參數(shù)處理,然后考慮添加數(shù)據(jù)和未知參數(shù)的后驗分布,未知參數(shù)的滿條件分布可轉化為較簡單的后驗分布,而添加數(shù)據(jù)的滿條件分布則是由添加數(shù)據(jù)的概率分布抽樣,最后分別對各滿條件分布進行Gibbs抽樣即可.

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

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

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

下面利用兩種方法隨機產生z1i.第一種方法稱為舍選法,步驟如下：

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

(1)分別產生來自f(x|θ),g(y),h(t)的隨機數(shù)x,y,t;

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

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

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

§3 失效率變點模型

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

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

(2)對于τ取無信息先驗分布：π(τ)∝1,τ＞0.

假設λ1,λ2,τ相互獨立,則

下面討論截斷變量和刪失變量都服從指數(shù)分布這種特殊情況.

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

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

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

(1)產生U(0,1)隨機數(shù)u1,u2,u3.

(2)計算

§4 EM算法

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

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

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

§5 MCMC方法

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

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

下面介紹一下截斷分布的抽樣.

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

則g(x)是由f(x)定義的在(a,b)上的截斷的分布.設F-1(x)為F(x)的反函數(shù),U為由U(0,1)產生的隨機數(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,∞)上截斷正態(tài)分布的密度函數(shù)為

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

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

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

先考慮左截斷右刪失模型.為了書寫方便,不妨假定后n0個壽命變量沒有觀察值,前n1=n-n0個壽命變量有觀察值z1＜z2＜···＜zn1,則似然函數(shù)為

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

§7 隨機模擬

下面利用R軟件進行隨機模擬試驗.

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

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

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

求解參數(shù)的MLE時,對τ所在的n1+1=119個區(qū)間分別建立對數(shù)似然方程組,得到119個解及其對應的對數(shù)似然函數(shù)值,最大的對數(shù)似然函數(shù)值對應的解即為要求的MLE.這119個解的τ分量及其對應的對數(shù)似然函數(shù)值的對照圖見圖13.

圖1 λ1的EM迭代過程

圖2 λ2的EM迭代過程

圖3 τ的EM迭代過程

圖4 λ1的Gibbs抽樣迭代過程

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

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

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

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

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

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

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

§8 結 論

圖5 λ2的Gibbs抽樣迭代過程

圖6 τ的Gibbs抽樣迭代過程

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

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

圖9 τ的兩條迭代鏈軌跡

圖10 λ1的后驗核密度估計

圖11 λ2的后驗核密度估計

圖12 τ的后驗核密度估計

本文主要利用Bayes方法研究了刪失截斷情形下失效率變點模型的參數(shù)估計問題.添加了部分缺失壽命變量數(shù)據(jù)并且討論了它的概率分布和隨機抽樣方法.利用EM算法和MCMC方法得到了變點位置參數(shù)和其它參數(shù)的Bayes估計,隨機模擬數(shù)據(jù)表明參數(shù)估計的精度較高.

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

圖13 τ與對數(shù)似然函數(shù)區(qū)間最大值的對照圖

表1 參數(shù)λ1,λ2,τ的Bayes估計和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ù)理統(tǒng)計(第二版)[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

河南省科技攻關計劃(162102310384);河南省高等學校重點科研項目(16A110001)