何朝兵
(安陽師范學院 數(shù)學與統(tǒng)計學院 河南 安陽 455000)
右刪失左截斷數(shù)據(jù)下離散威布爾分布的參數(shù)估計
何朝兵
(安陽師范學院 數(shù)學與統(tǒng)計學院 河南 安陽 455000)
研究了右刪失左截斷數(shù)據(jù)模型下離散威布爾分布參數(shù)的極大似然估計和漸近置信區(qū)間.介紹了參數(shù)估計的牛頓迭代方法和EM算法,給出了參數(shù)的漸近置信區(qū)間.隨機模擬的結果表明,牛頓迭代方法和EM算法得到的參數(shù)估計結果差別不大.
極大似然估計; 牛頓迭代方法; EM算法; 缺損信息原則; 漸近置信區(qū)間
離散威布爾分布是一種很重要的離散型壽命分布[1-3],它是威布爾分布的離散化,在排隊論和可靠性數(shù)學等分支中有著很廣泛的應用,并且?guī)缀畏植际翘厥獾碾x散威布爾分布.當觀察壽命數(shù)據(jù)時,右刪失與左截斷情形經(jīng)常同時發(fā)生.例如,對確診為艾滋病的患者進行觀察研究,如果患者仍然存活或提前退出研究,這就是右刪失;如果患者在研究開始之前就已經(jīng)死亡,這就是左截斷.右刪失左截斷數(shù)據(jù)模型廣泛應用于醫(yī)學、生物學、經(jīng)濟學等領域,近年來對此模型的研究比較熱.文獻[4-8]研究了此模型下壽命分布為連續(xù)型的情形,但對離散威布爾分布情形的研究尚未見報道.本文主要利用牛頓迭代(NR)方法和EM算法研究了右刪失左截斷數(shù)據(jù)模型下離散威布爾分布參數(shù)的極大似然估計和區(qū)間估計,隨機模擬的結果表明,由這兩種方法得到的參數(shù)估計結果差別不大.
設(X,Y,T)是離散型隨機變量,X的分布函數(shù)為F(x,θ)=P(Xi≤x),分布律為f(x,θ),其中θ是向量參數(shù);Y是右刪失隨機變量,分布函數(shù)為G(y),分布律為g(y);T是左截斷隨機變量,分布函數(shù)為H(t),分布律為h(t),且Y,T的分布與參數(shù)θ無關.假定X,Y,T是相互獨立且取正整數(shù)的隨機變量,X是所感興趣的隨機變量.右刪失左截斷數(shù)據(jù)的試驗模型是:僅在Zi≥Ti時得到觀察數(shù)據(jù)(Zi,Ti,δi),而在Zi 求以下樣本的似然函數(shù): P(無樣本觀察值)=P(Zi 為了敘述與書寫方便,假定前n1個樣本有觀察值,剩下的n2個樣本沒有觀察值(n1+n2=n),則基于數(shù)據(jù){(Zi,Ti,δi),1≤i≤n1}的似然函數(shù)為 若X的分布函數(shù)為F(x;α,β)=1-e-(x/β)α,x>0,且α>0,β>0,則稱X服從尺度參數(shù)為β,形狀參數(shù)為α的威布爾分布,記為X~Wei(α,β).若X的分布律為P(X=k)=q(k-1)α-qkα,k=1,2,…,且0 假設右刪失左截斷數(shù)據(jù)試驗中,所感興趣的變量X服從形狀參數(shù)為α,尺度參數(shù)為q的離散威布爾分布,此時θ=(α,q),下面介紹參數(shù)α與q的估計方法. 2.1 牛頓迭代(NR)方法 基于數(shù)據(jù){(Zi,Ti,δi),1≤i≤n1}的似然函數(shù)為 其對數(shù)似然函數(shù)為 令其兩個偏導數(shù)為零,得到的兩個似然方程為超越方程,只能用數(shù)值方法求解,本文利用NR方法求解. NR方法是通過使似然函數(shù)最大化而獲得極大似然估計的一種直接方法,它需要計算對數(shù)似然函數(shù)關于參數(shù)的一階和二階偏導數(shù),所以十分煩瑣.這里,通過使用R軟件maxLik程序包中的maxNR函數(shù)來獲得q和α的極大似然估計. 注 如果Y服從幾何分布Geo(p0),T服從取值1,2,…,s的離散均勻分布,則 2.2 EM算法 EM算法是一種迭代方法,最初由文獻[9]提出,主要用來求后驗分布的眾數(shù).它的每一步迭代由兩步組成:E步(求期望)和M步(極大化).EM算法處理不完全數(shù)據(jù)非常方便[10-11],下面用EM算法來求α和q的極大似然估計. 若第i個樣本沒有觀察值,添加其觀察值為(Wi,βi),其中:Wi=Xi∧Yi=min(Xi,Yi),βi=I(Xi≤Yi),i=n1+1,n1+2,…,n,則 可得似然函數(shù) 取(α,q)的先驗分布為無信息先驗分布π(α,q)=c,0 E步: 在給定θ,δ,Z和T下,(Wi,βi)的分布律為 則(Wi,αi)關于Wi的邊緣分布律為 P(Wi=k)=ψ1(k,θ)+ψ2(k,θ)ψ3(k,θ),k=1,2,…. 所以 注 如果Y服從幾何分布Geo(p0),T服從取值1,2,…,s的離散均勻分布,則 2.3 極大似然估計的漸近方差和協(xié)方差 利用NR方法求極大似然估計時,I0可以直接得到.而利用EM算法時,由于似然函數(shù)為基于完全數(shù)據(jù)的似然函數(shù),所以I0不能直接得到,但根據(jù)缺損信息原則可以獲得觀察信息陣I0.缺損信息原則為:觀察信息=完全信息-缺損信息.下面求基于EM算法的觀察信息,完全數(shù)據(jù)的似然函數(shù)為 為了方便書寫,記 φ2=n2{βln[q(w-1)α-wα-1]+wαlnq}, ψ(β,w)βψ1(w)+(1-β)ψ2(w),β=0,1, lnL2n2lnψ(β,w). 則完全信息為 (1) 損失信息為 (2) 根據(jù)缺損信息原則,可得觀察信息為 2.4 置信區(qū)間 式中:bq和σq分別為q的1 000個估計值的bootstrap偏差和方差;zβ/2為標準正態(tài)分布的上β/2分位點.α的參數(shù)bootstrap置信區(qū)間的構造與q的類似. 表1 NR方法和EM算法下參數(shù)估計的均值(M)、偏差(B)、均方誤差(MSE)和置信區(qū)間的覆蓋率(CP) 從表1可以看出,通過NR方法和EM算法得到的參數(shù)估計的均值、偏差、均方誤差以及置信水平為0.95的置信區(qū)間的覆蓋率都比較接近,說明這兩種方法差別不大.同時可以看出,樣本容量對估計值的影響也不大,說明得到的估計值是比較穩(wěn)定的,并且精度也較高. 表2 由3種方法構造的置信區(qū)間 由表2可以看出,通過NR方法和EM算法構造的置信區(qū)間是一樣的,雖然它們與參數(shù)bootstrap置信區(qū)間不一樣,但差別不是很大. [1] NEKOUKHOU V, BIDRAM H. The exponentiated discrete Weibull distribution[J]. SORT, 2015, 39(1): 127-146. [2] ALMALKI S J, NADARAJAH S. Modifications of the Weibull distribution: a review[J]. Reliab Eng Syst Safe, 2014, 124: 32-55. [3] ENGLEHARDT J D, LI R. The discrete Weibull distribution: an alternative for correlated counts with confirmation for microbial counts in water[J]. Risk Anal, 2011, 31(3): 370-381. [4] DEWAN I. Comments: EM-based likelihood inference for some lifetime distributions based on left truncated and right censored data and associated model discrimination[J]. S Afr Stat J, 2014, 48(2): 183-185. [5] SHEN P S. Aalen’s additive risk model for left-truncated and right-censored data[J]. Commun Stat Simul Comput, 2014, 43(5): 1006-1019. [6] BALAKRISHNAN N, MITRA D. Likelihood inference based on left truncated and right censored data from a gamma distribution[J]. IEEE Trans Rel, 2013, 62(3): 679-688. [7] SU Y R, WANG J L. Modeling left-truncated and right-censored survival data with longitudinal covariates[J]. Ann Stat, 2012, 40(3): 1465-1488. [8] AHMADI J, DOOSTPARAST M, PARSIAN A. Estimation with left-truncated and right censored data: a comparison study[J]. Stat Probab Lett, 2012, 82(7): 1391-1400. [9] DEMPSTER A P, LAIRD N M, RUBIN D B. Maximum likelihood from incomplete data via the EM algorithm[J]. J R Stat Soc, 1977, 39(1): 1-38. [10] CHUNG Y, LINDSAY B G. Convergence of the EM algorithm for continuous mixing distributions[J]. Stat Probab Lett, 2015, 96(1): 190-195. [11]GRIGOROVA D, ENCHEVA E, GUEORGUIEVA R. EM algorithm for MLE of a probit model for multiple ordinal outcomes[J]. Serdica J Comput, 2013, 7(3): 227-244. (責任編輯:孔 薇) Parameter Estimations of Discrete Weibull Distribution with Left Truncated and Right Censored Data HE Chaobing (SchoolofMathematicsandStatistics,AnyangNormalUniversity,Anyang455000,China) The maximum likelihood estimation and asymptotic confidence intervals of the parameters of discrete Weibull distribution were mainly studied with left truncated and right censored data. Newton-Raphson method and EM algorithm of parameter estimation were introduced, and the asymptotic confidence intervals of the parameters were given. Random simulation test results showed that there was little difference on parameter estimation between Newton-Raphson method and EM algorithm. maximum likelihood estimation; Newton-Raphson method; EM algorithm; missing information principle; asymptotic confidence interval 2015-10-18 國家自然科學基金資助項目(61174099);河南省高等學校重點科研項目(16A110001). 何朝兵(1975—),男,河南周口人,講師,碩士,主要從事概率統(tǒng)計研究,E-mail:chaobing5@163.com. 何朝兵.右刪失左截斷數(shù)據(jù)下離散威布爾分布的參數(shù)估計[J]. 鄭州大學學報(理學版), 2016,48(2): 18-23. O213.2 A 1671-6841(2016)02-0018-06 10.13705/j.issn.1671-6841.20152202 離散威布爾分布的參數(shù)估計方法
0,則稱X服從尺度參數(shù)為q,形狀參數(shù)為α的離散威布爾分布,當α=1時,歸結為幾何分布Geo(1-q).容易證明:若X~Wei(α,β),且β=(-1/lnq)1/α,則Y=[X]+1服從尺度參數(shù)為q,形狀參數(shù)為α的離散威布爾分布.利用此結論可以產(chǎn)生離散威布爾分布的隨機數(shù).
0,c為常數(shù),則θ=(α,q)的添加后驗分布為
3 隨機模擬結果