刪失截?cái)嗲樾蜗耊eibull分布多變點(diǎn)模型的參數(shù)估計(jì)

刪失截?cái)嗲樾蜗耊eibull分布多變點(diǎn)模型的參數(shù)估計(jì)何朝兵【期刊名稱】《《高校應(yīng)用數(shù)學(xué)學(xué)報(bào)A輯》》【年(卷)，期】2015(000)002【總頁(yè)數(shù)】12頁(yè)(P127-138)【關(guān)鍵詞】完全數(shù)據(jù)似然函數(shù);滿條件分布;MCMC方法;Gibbs抽樣;Metropolis-Hastings算法【作者】何朝兵【作者單位】安陽(yáng)師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，河南安陽(yáng)455000【正文語(yǔ)種】中文【中圖分類】O213.2;O212.8近年來(lái)變點(diǎn)問(wèn)題成為統(tǒng)計(jì)學(xué)中比較熱門的研究方向，它廣泛應(yīng)用于工業(yè)質(zhì)量控制、水文統(tǒng)計(jì)、經(jīng)濟(jì)和地震預(yù)測(cè)等領(lǐng)域.在變點(diǎn)模型中，將樣本觀察值按出現(xiàn)時(shí)間的先后順序排列，如果在某個(gè)未知時(shí)刻樣本的分布或數(shù)字特征突然改變，那么這個(gè)未知時(shí)刻被稱為變點(diǎn).變點(diǎn)檢測(cè)方法主要有最大似然法、非參數(shù)方法和Bayes方法等.關(guān)于變點(diǎn)模型的研究可參看文獻(xiàn)［1-7］,其中文［1］首次研究了變點(diǎn)模型，文［2］總結(jié)了變點(diǎn)分析理論，文［3］對(duì)變點(diǎn)模型進(jìn)行了準(zhǔn)確高效的Bayes推斷，文［4］給出了水文氣象時(shí)間序列變點(diǎn)問(wèn)題的Bayes變點(diǎn)分析.最大似然法主要依靠樣本的似然函數(shù)來(lái)進(jìn)行統(tǒng)計(jì)推斷，但由于變點(diǎn)位置參數(shù)和分布參數(shù)都是未知的,所以似然函數(shù)很復(fù)雜，這時(shí)似然方程的求解往往要借助于繁瑣的數(shù)值計(jì)算方法來(lái)完成.非參數(shù)方法一般會(huì)繞過(guò)一個(gè)似然函數(shù)的完全規(guī)定，其后果是推斷缺乏效率，估計(jì)有偏.Bayes方法由于能夠綜合各類先驗(yàn)信息，給出變點(diǎn)可能位置等優(yōu)點(diǎn)在國(guó)內(nèi)外得到廣泛應(yīng)用，而計(jì)算的復(fù)雜性是Bayes方法的難點(diǎn).MarkovchainMonteCarlo(MCMC)方法是一種簡(jiǎn)單且行之有效的Bayes計(jì)算方法,主要包括Gibbs抽樣和Metropolis-Hastings算法.MCMC方法使Bayes統(tǒng)計(jì)中許多看起來(lái)困難的計(jì)算變得簡(jiǎn)單直觀.MCMC算法的應(yīng)用很廣泛，比如約束參數(shù)模型，變點(diǎn)模型,截尾數(shù)據(jù)和分組數(shù)據(jù)，多個(gè)分布的混合等.關(guān)于MCMC方法在變點(diǎn)模型中的應(yīng)用可參看文獻(xiàn)［8-11］,其中文［8］給出了MCMC方法在信號(hào)變點(diǎn)問(wèn)題中的應(yīng)用，文［9］利用模擬退火方法研究了Bayes多變點(diǎn)估計(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è)體在研究開始之前就已經(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è)體在研究還未開始之前就已經(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)［12-15］,其中文［12］研究了此模型的回歸分析和非參數(shù)估計(jì)，文［13］研究了此模型下對(duì)數(shù)正態(tài)分布的最大然估計(jì).Weibull分布是可靠性中常用失效分布之一,許多電子與機(jī)械的元器件與設(shè)備等產(chǎn)品的失效分布都是Weibull分布.大量的實(shí)踐說(shuō)明，凡是因某一局部失效而導(dǎo)致全局停止運(yùn)行的元件、器件、設(shè)備等的壽命都可看做或近似看做服從Weibull分布.譬如金屬材料和機(jī)械零件或部件的疲勞壽命就服從Weibull分布.由于變點(diǎn)模型和刪失截?cái)鄶?shù)據(jù)模型的廣泛應(yīng)用，所以對(duì)刪失截?cái)嗲樾蜗伦凕c(diǎn)問(wèn)題的研究是一個(gè)有意義的研究方向,具有重要的理論價(jià)值和實(shí)用價(jià)值.由于此模型是截尾數(shù)據(jù)模型和變點(diǎn)模型的混合，所以Bayes方法中的MCMC方法很適合研究此變點(diǎn)模型.關(guān)于截?cái)嗷騽h失數(shù)據(jù)下變點(diǎn)問(wèn)題的研究已有一些成果，可參看文獻(xiàn)［16-19］,其中文［16］和［17］分別研究了刪失數(shù)據(jù)下風(fēng)險(xiǎn)率變點(diǎn)問(wèn)題的非參數(shù)和參數(shù)估計(jì).文［20］研究了刪失截?cái)嗲樾蜗嘛L(fēng)險(xiǎn)率變點(diǎn)問(wèn)題,但對(duì)刪失截?cái)嗲樾蜗耊eibull分布多變點(diǎn)模型的研究還不多見(jiàn).下文主要利用基于MCMC方法的Bayes方法研究了刪失截?cái)嗲樾蜗耊eibull分布多變點(diǎn)模型的參數(shù)估計(jì)問(wèn)題.把Gibbs樣本的均值作為各參數(shù)的Bayes估計(jì),最后進(jìn)行了隨機(jī)模擬試驗(yàn)，試驗(yàn)結(jié)果表明變點(diǎn)位置參數(shù)和形狀參數(shù)以及尺度參數(shù)Bayes估計(jì)的精度都較高.設(shè)(X,Y,T)是一連續(xù)型隨機(jī)變量，壽命變量X的分布函數(shù)為F(x;0)=P(X<x),密度函數(shù)為f(x;0),這里0是未知參數(shù)向量;Y是一右刪失隨機(jī)變量，分布函數(shù)為G(y),密度函數(shù)為g(y)；T是一左截?cái)嚯S機(jī)變量，分布函數(shù)為H(t),密度函數(shù)為h(t),且Y,T的分布與參數(shù)0無(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,T侄i)的密度函數(shù).當(dāng)Zi=zi,Ti=ti,&=1時(shí)(Zi,Ti,&)的密度函數(shù)為f(zi;0)(zi)h(ti),zi>ti;當(dāng)Zi=zi,Ti=ti,&=0時(shí)(Zi,Ti,&)的密度函數(shù)為g(zi)(zi;0)h(ti),zi>ti.為了研究方便，引入示性變量vi=I(min(Xi,Yi)>Ti),i=12…,n.則基于觀察數(shù)據(jù)｛(zi,t侄i):vi=1,1<i<n｝的似然函數(shù)為若X的密度函數(shù)為F(x;m,入)=1-e-入xm,x>0,m>0,入〉0,則稱X服從參數(shù)為0=(m,A)的Weibull分布，記為X~W(m,入),m稱為形狀參數(shù)，入稱為尺度參數(shù).Weibull分布的密度函數(shù)為在刪失截?cái)嘣囼?yàn)中，假設(shè)產(chǎn)品壽命X~W(m,入).由(1)式易知，刪失截?cái)嗲樾蜗耊eibull分布的似然函數(shù)比較復(fù)雜，所以下面添加缺失的Xi的值以獲得較簡(jiǎn)單的完全數(shù)據(jù)似然函數(shù).具體如下：當(dāng)vi=1,&=0時(shí)Xi有缺失,僅知道Xi>zi,此時(shí)添加數(shù)據(jù)Xi=Z1i=z1i.在Xi>zi的條件下,Z1i的條件密度函數(shù)為Z1i的樣本可表示為z1i=F-1{F(zi)+U[1-F(zi)]},其中F是W(m,入)的分布函數(shù),F-1為其反函數(shù),U為來(lái)自均勻分布U(0,1)的一個(gè)隨機(jī)樣本.當(dāng)進(jìn)行隨機(jī)模擬編寫R程序時(shí),F與F-1可以分別利用函數(shù)pweibull和qweibull進(jìn)行計(jì)算.注當(dāng)vi=1,&=0,即數(shù)據(jù)右刪失時(shí)，上文采取的添加數(shù)據(jù)的方法是隨機(jī)的,實(shí)際上還有其他添加數(shù)據(jù)的方法，如Buckley和James把右刪失數(shù)據(jù)添加為E(Xi|Xi>zi),這時(shí)添加的數(shù)據(jù)為常數(shù)[21].當(dāng)vi=0時(shí)Xi有缺失,僅知道m(xù)in(Xi,Yi)<Ti,此時(shí)添加數(shù)據(jù)Xi=Z2i=z2i..在min(Xi,Yi)<Ti的條件下,Z2i的條件密度函數(shù)為由于f(x;0)是Xi的密度函數(shù)，且Xi與Z2i有相同的取值，符合篩選抽樣的適用條件，故可以利用篩選法隨機(jī)產(chǎn)生Z2i的值z(mì)2i.z2i抽樣的具體步驟如下：由均勻分布U(0,1)抽取u,由f(x;0)抽取x;如果，則z2i=x,停止;否則回到步驟(1).令&z,u1,u2,v分別表示由筑zi,z1i,z2i,vi組成的向量，則完全數(shù)據(jù)似然函數(shù)為當(dāng)對(duì)Weibull壽命產(chǎn)品進(jìn)行壽命試驗(yàn)時(shí)，由于某些原因?qū)е耊eibull分布的參數(shù)發(fā)生改變，這時(shí)如果仍在原來(lái)模型的基礎(chǔ)上進(jìn)行統(tǒng)計(jì)推斷，那么得到的參數(shù)估計(jì)就會(huì)與實(shí)際情況差別很大，并且往往是錯(cuò)誤的，錯(cuò)誤的統(tǒng)計(jì)推斷就會(huì)給工業(yè)生產(chǎn)帶來(lái)錯(cuò)誤的指導(dǎo).所以，對(duì)刪失截?cái)鄶?shù)據(jù)下Weibull壽命多變點(diǎn)問(wèn)題的研究就顯得尤為重要.Weibull分布多變點(diǎn)模型如下：其中(m1,入1),(m2,入2),(m3,入3)兩兩不相等,1<k1<k2<n-1.下面利用Bayes方法對(duì)變點(diǎn)位置k1,k2及參數(shù)入1,入2,入3,m1,m2,m3進(jìn)行估計(jì).令D1={1,2,…,k1},D2={k1+1,…,k2},D3={k2+1，…,n}.記a=(k1,k2,入1,入2,入3,m1,m2,m3),由(2)式和(3)式可得此變點(diǎn)模型的似然函數(shù)為下面確定各參數(shù)的先驗(yàn)分布，記0的先驗(yàn)分布為n(0).(1)易知1<k1<k2<n-1,如果k1與k2無(wú)任何先驗(yàn)信息時(shí)，可把k1與k2看做從1,2,…,n-1中無(wú)放回等可能抽取的兩個(gè)數(shù)，即取k1與k2的聯(lián)合先驗(yàn)分布為基于上面的討論，下面使用R軟件進(jìn)行隨機(jī)模擬試驗(yàn).取兩組受試樣品,第一組受試樣品的個(gè)數(shù)n=200,第二組受試樣品的個(gè)數(shù)n=400.n=200時(shí)右刪失變量Yi~W(0.7,2.5),左截?cái)嘧兞縏i~W(4,0.9).則參數(shù)(k1,k2,入1,入2,入3,m1,m2,m3)的真實(shí)值為(35,11020.8,3,0.651.8).利用日語(yǔ)言編程時(shí),Xi,Yi,Ti的抽樣可通過(guò)函數(shù)rweibull實(shí)現(xiàn).n=400時(shí)右刪失變量Yi~W(3,6.5),左截?cái)嘧兞縏i~W(7,2).則參數(shù)(k1,k2,入1,入2,入3,m1,m2,m3)的真實(shí)值為(120,350,1.5,8,3.5,6,2.5,7.5).根據(jù)各個(gè)參數(shù)的滿條件分布使用R軟件進(jìn)行MCMC模擬，主要對(duì)變點(diǎn)位置參數(shù)進(jìn)行估計(jì)分析.在模擬運(yùn)行過(guò)程中，先進(jìn)行10000次Gibbs預(yù)迭代，以確保抽樣的收斂性，然后丟棄最初的預(yù)迭代，再進(jìn)行10000次Gibbs迭代.迭代從第10001次開始至第20000次的R程序的運(yùn)行結(jié)果見(jiàn)表1和表2.變點(diǎn)位置參數(shù)的Gibbs抽樣迭代過(guò)程見(jiàn)圖1至圖4.在模型的分析過(guò)程中,MCMC的收斂性診斷很重要，模擬時(shí)可以對(duì)參數(shù)進(jìn)行多層鏈?zhǔn)降治?，即輸入多組初始值，形成多層迭代鏈,當(dāng)抽樣收斂時(shí),迭代圖形重合.在模擬過(guò)程中,輸入兩組初始值分別進(jìn)行10000次迭代.變點(diǎn)位置參數(shù)的多層迭代鏈軌跡見(jiàn)圖5至圖8.最后，進(jìn)行模擬結(jié)果分析.首先，由表1和表2可以看出各參數(shù)的Bayes估計(jì)與真值的相對(duì)誤差不超過(guò)5%,MC誤差也很小，所以各參數(shù)Bayes估計(jì)的精度較高,估計(jì)效果較好.并且Gibbs樣本的中位數(shù)與均值差別很小，所以把樣本中位數(shù)作為參數(shù)的Bayes估計(jì)的精度也較高.各參數(shù)可信水平為0.95的可信區(qū)間可近似取為［2.5%分位數(shù),97.5%分位數(shù)］，可以看出近似可信區(qū)間的長(zhǎng)度非常短，所以區(qū)間估計(jì)的效果也較好.由圖1至圖4可以看出Gibbs抽樣波動(dòng)較小，并且很快就趨于平穩(wěn)，呈現(xiàn)出很強(qiáng)的規(guī)律性;其次，要判斷所產(chǎn)生的Markov鏈?zhǔn)欠袷諗?由圖5至圖8可以發(fā)現(xiàn),變點(diǎn)位置參數(shù)的兩條迭代鏈都趨于重合,收斂性較好.綜上分析,可以看出通過(guò)MCMC模擬所產(chǎn)生的效果較好.綜上所述，通過(guò)模擬試驗(yàn)得到的各參數(shù)的Bayes估計(jì)的精度較高，效果較好,并且操作方便.這說(shuō)明了MCMC方法非常適合處理刪失截?cái)嗲樾蜗耊eibull分布多變點(diǎn)問(wèn)題.由于刪失截?cái)鄶?shù)據(jù)下Weibull分布變點(diǎn)問(wèn)題的似然函數(shù)非常復(fù)雜,甚至寫不出顯式的似然函數(shù)，并且參數(shù)較多,所以最大似然方法不太適合處理該模型;另外非參數(shù)方法一般是基于壽命分布類型未知情況下才使用的,并且非參數(shù)方法的精度往往較低，并且推斷缺乏效率.所以與其他方法相比,Bayes方法是一種簡(jiǎn)單且行之有效的方法,并且通過(guò)模擬試驗(yàn)驗(yàn)證了其有較高的精度和較好的效果.論文主要利用Bayes方法研究了刪失截?cái)嗲樾蜗耊eibull分布多變點(diǎn)模型中的參數(shù)估計(jì)問(wèn)題.利用MCMC方法得到了變點(diǎn)位置參數(shù)和形狀參數(shù)以及尺度參數(shù)的Bayes估計(jì)，隨機(jī)模擬數(shù)據(jù)表明Bayes估計(jì)的精度較高.論文中所介紹的統(tǒng)計(jì)推斷是基于Gibbs樣本而進(jìn)行的，所以論文所介紹的Bayes估計(jì)方法具有一般性，同樣適用其他壽命分布，具有一定的推廣價(jià)值.論文研究的模型是截尾數(shù)據(jù)模型和變點(diǎn)模型的混合，而MCMC方法非常適合處理截尾數(shù)據(jù)和變點(diǎn)問(wèn)題，這也是論文的優(yōu)點(diǎn).但MCMC方法也有缺點(diǎn)，例如,MCMC的收斂速度比較慢，有時(shí)MCMC是不收斂的,這些都會(huì)影響MCMC方法在實(shí)際問(wèn)題中的應(yīng)用.論文研究多變點(diǎn)模型時(shí)，假定變點(diǎn)個(gè)數(shù)已知.而如何根據(jù)觀察數(shù)據(jù)確認(rèn)變點(diǎn)個(gè)數(shù)，是論文尚待解決的問(wèn)題，是論文的不足之處.實(shí)際上，變點(diǎn)個(gè)數(shù)的確認(rèn)是現(xiàn)有變點(diǎn)問(wèn)題研究中的難題，也是變點(diǎn)問(wèn)題研究的發(fā)展趨勢(shì).【相關(guān)文獻(xiàn)】PageES.Continuousinspectionschemes[J].Biometrika,1954,41(1):100-115.CsrgM,HorvthL.Limittheoremsinchange-pointanalysis[M].NewYork:Wiley,1997.FearnheadP.ExactandefficientBayesianinferenceformultiplechangepointproblems[J].Statisticsandcomputing,2006,16(2):203-213.PerreaultL,BernierJ,BobeB,etal.Bayesianchange-pointanalysisinhydrometeorologicaltimeseries.Part1.Thenormalmodelrevisited[J].JournalofHydrology,2000,235(3):221-241.陳希孺.變點(diǎn)統(tǒng)計(jì)分析簡(jiǎn)介[J].數(shù)理統(tǒng)計(jì)與管理,1991,10(1):55-59.熊立華，周芬.水文時(shí)間序列變點(diǎn)分析的Bayes方法[J].水電能源科學(xué),2003,21(4):39-41.KotzS,吳喜之.現(xiàn)代Bayes統(tǒng)計(jì)學(xué)[M].北京:中國(guó)統(tǒng)計(jì)出版社,2000.LavielleM,LebarbierE.AnapplicationofMCMCmethodsforthemultiplechange-pointsproblem[J].SignalProcessing,2001,81(1):39-53.KimJ,CheonS.Bayesianmultiplechange-pointestimationwithannealingstochasticapproximationMonteCarlo[J].ComputationalStatistics,2010,25(2):215-239.YuanTao,KuoYue.Bayesiananalysisofhazardrate,changepoint,andcost-optimalburnintimeforelectronicdevices[J].IEEETransactionsonReliability2010,59(1):132-138.LiangFaming,WongWH.Real-parameterevolutionaryMonteCarlowithapplicationstoBayesianmixturemodels[J].JournaloftheAmericanStatisticalAssociation,2001,96(454):653-666.GrossST,LaiTL.Nonparametricestimationandregressionanalysiswithleft-truncatedandright-censoreddata[J].JournaloftheAmericanStatisticalAssociation,1996,91(435):1166-1180.BalakrishnanN,MitraD.Likelihoodinferenceforlognormaldatawithlefttruncationandrightcensoringwithanillustration[J].JournalofStatisticalPlanningandInference,2011,141(11):3536-3553.CosslettSR.Efficientsemiparametricestimationofcensoredandtruncatedregressionsviaasmoothedself-consistencyequation[J].Econometrica,2004,72(4):1277-1293.Molanes-lopezEM,CaoR,KeilegomIVAN.Smoothedempiricallikelihoodconfidenceintervalsfortherelativedistributionwithleft-truncatedandright-censoreddata[J].CanadianJournalofStatistics,2010,38(3):453-473.AntoniadisA,GrgoireG.Nonparametricestimationinchangepointhazardratemodelsforcensoreddata:Acountingprocessapproach[J].JournalofNonparametricStatistics,1993,3(2):135-154.GijbelsI,Grler.Estimationofachangepo