何朝兵,劉華文
(1.安陽師范學院 數學與統(tǒng)計學院,河南 安陽 455000;2.山東大學 數學學院,山東 濟南 2501000)
IIRCT下泊松分布參數單變點的貝葉斯估計
何朝兵1,劉華文2
(1.安陽師范學院 數學與統(tǒng)計學院,河南 安陽 455000;2.山東大學 數學學院,山東 濟南 2501000)
首先通過添加數據得到了帶有不完全信息隨機截尾試驗下泊松分布的完全數據似然函數,然后研究了變點位置和其它參數的滿條件分布,接著利用Gibbs抽樣與Metropolis-Hastings算法相結合的MCMC方法對參數進行了估計,最后進行了隨機模擬,試驗結果表明參數貝葉斯估計的精度較高.
完全數據似然函數;滿條件分布;MCMC方法;Gibbs抽樣;Metropolis-Hastings算法
變點分析研究始于20世紀50年代,自20世紀70年代初,這一問題引起了一些編譯學家的重視,發(fā)展了一些變點分析方法[1-3],如最小二乘法、極大似然法、貝葉斯法和非參數方法等.隨著統(tǒng)計計算技術的快速發(fā)展,貝葉斯方法的應用越來越廣泛,特別是在水文統(tǒng)計[4,5]和計量經濟[6]等領域.而貝葉斯計算方法中的Markov chain Monte Carlo(MCMC)方法,使變點分析中貝葉斯方法的實際操作變得非常方便.泊松分布是一類應用廣泛的離散型壽命分布,它可以作為多種不同類型試驗的模型,泊松分布的另一個應用是研究空間分布,例如湖泊中魚群數量的分布.文獻[7-12]對泊松分布進行了深入研究,其中文獻[12]研究了II型截尾情形下泊松分布參數的極大似然估計.近些年來,關于隨機截尾試驗的研究比較多,帶有不完全信息隨機截尾試驗(random censoring test model with incomplete information),簡稱IIRCT,由文獻[13]首次涉及,接著文獻[14-20]深入研究了IIRCT下連續(xù)型分布的參數估計.實際上,文獻[12]研究的II型截尾試驗是特殊的IIRCT.而對IIRCT下壽命分布參數的變點問題,卻很少有文獻研究.下文主要研究了IIRCT下泊松分布參數單變點的貝葉斯估計.首先通過添加數據得到了IIRCT下泊松分布的完全數據似然函數,然后研究了變點位置和其它參數的滿條件分布,接著利用Gibbs抽樣與Metropilis-Hastings算法相結合的MCMC方法得到了參數的Gibbs樣本,把Gibbs樣本的均值作為各參數的貝葉斯估計,最后進行了隨機模擬,試驗結果表明各參數貝葉斯估計的精度都較高.
設產品壽命X1,X2,L是獨立同分布的離散型隨機變量序列,其分布函數為F(x,λ)=P(Xi≤x),分布為f(x,λ),λ是未知參數.又設Y1,Y2,L是獨立的取非負整數的離散型隨機變量序列,分布函數分別為G1(y),G2(y),L.分布律分別為g1(y),g2(y),L,且gi(y)與參數λ無關.{Xi}與{Yi}獨立.
為估計參數λ,n個樣品的觀察數據{Zi,1≤i≤n}如下:
(1)當Xi≤Yi時,有兩種情況:Xi以概率ai立即顯示,此時取Zi=Xi;Xi以概率1-ai不被顯示,此時取Zi=Yi,稱ai為失效顯示概率.
(2)當Xi>Yi時,取Zi=Yi.
引入隨機變量αi,βi,i=1,2,L,n.
若Xi≤Yi,αi=1;若Xi>Yi,αi=0.
若Xi≤Yi,且未被顯示時,βi=0;其它情況,βi=1.
綜上所述,有
設Zi的觀察值為zi,則基于數據{(zi,αi,βi),1≤i≤n}的似然函數為
若X的分布律為P(X=k)=λke-λ/k!,k=0,1,L,則稱X服從參數為λ的泊松分布,記為X~P(λ).假設IIRCT下產品壽命X~P(λ).
令α表示αi組成的向量,β表示βi組成的向量,z表示z組成的向量.
泊松分布P(λ)基于數據{(zi,α,βi),1≤i≤n}的似然函數為
從上式可以看出,基于截尾數據的似然函數比較復雜,為了方便進行參數估計,下面添加缺損的Xi的值,以獲得完全數據的似然函數.具體如下:
當αi=1,βi=0,即第i個產品失效未被顯示時,添加數據Z1i=Xi=z1i.
在Xi≤zi的條件下,Z1i的條件分布律為
當αi=0,即Xi>Yi時,添加數據Z2i=Xi=z2i.
在Xi>zi的條件下,Z2i的條件分布律為
令u表示z1i組成的向量,v表示z2i組成的向量,則完全數據似然函數為
泊松分布參數單變點模型如下:
其中λ1,λ2兩兩不相等,1≤k≤n-1.
下面對變點位置k及參數λ1,λ2進行貝葉斯估計.
記θ=(k,λ1,λ2),則IIRCT下泊松分布參數單變點問題的似然函數為
下面確定各參數的先驗分布.
(2)取λi的先驗分布為共軛先驗分布伽瑪分布Ga(bi,ci),bi,ci已知,即
假設k,λ,λ2獨立.則
L(θ|z,u,v,α,β)∝π(k)π(λ1)π(λ2)L(z,u,v,α,β|θ)
當αi=1,βi=0時,
其中z-1i={z1j:j≠i}.
當αi=0時,
其中z-2i={z2j:j≠i}.
為了書寫方便,把滿條件分布中的“條件”用“·”代替,例如
π(λ1|k,λ2,z,u,v,α,β)簡記為π(λ1|·).
各參數的滿條件分布為
由于得到了各參數的滿條件分布,下面利用MCMC方法獲得各參數后驗分布的平穩(wěn)分布.參數z1i,z2i的滿條件分布是截斷泊松分布,可以利用逆變換法隨機抽樣,λ1,λ2的滿條件分布是伽瑪分布,這4個分布都可以采用Gibbs抽樣;但是k的滿條件分布不是標準分布,進行Gibbs抽樣比較困難,可以利用Metropolis-Hastings算法進行抽樣,此時選取建議分布為取值1,2,L,n-1的離散型均勻分布.
下面寫出MCMC方法的具體步驟:
從1,2,L,n-1中任取一個k′,產生一個隨機數u0,若u0≤r(k(t-1),k′),則k(t)=k′,否則k(t)=k(t-1).
基于上面的討論,下面進行隨機模擬試驗.取受試樣品的個數n=200,樣品壽命Xi:P(10),i=1,2,L,50;Xi:P(30),i=51,52,L,200.則參數(k,λ1,λ2)的真實值為(50,10,30).取λ1,λ2的先驗分布分別為Ga(7,0.6),Ga(35,1.2),截尾變量Yi~P(32),失效顯示概率a=0.8.
利用前面推導出的各個參數的滿條件分布使用R軟件進行MCMC模擬.在模擬運行過程中,先進行10000次Gibbs預迭代,以確保參數的收斂性,然后丟棄最初的預迭代,再進行10000次Gibbs迭代.迭代從第10001次開始至第20000次的R程序的運行結果見表1.
表1 參數k,λ1,λ2的貝葉斯估計
最后,進行模擬結果分析.由表1可以看出各參數的估計值與真值的相對誤差都小于3%,精度較高,效果較好;各參數的MC誤差較小,Gibbs迭代值的波動較小,Gibbs迭代比較穩(wěn)定.綜上分析,可以看出通過MCMC模擬所產生的效果較好.
注:編寫R程序時用到的函數主要有pois(),gamma(),min(),runif(),apply(),sum(),mean(),sd(),quantile(),plot(),points(),legend()等.
[1] CS?RG? M, HORVTH L. Limit theorems in change-point analysis[M]. New York: Wiley, 1997.
[2] 陳希孺.變點統(tǒng)計分析簡介[J].數理統(tǒng)計與管理,1991,10(1):55-59.
[3] 項靜怡,史久恩.非線性系統(tǒng)中數據處理的統(tǒng)計方法[M].北京:科學出版社,1997.
[4] 熊立華,周芬.水文時間序列變點分析的貝葉斯方法[J].水電能源科學,2003,21(4):39-41.
[5] PERREAULT L, BERNIER J, BOBéE B, et al. Bayesian change-point analysis in hydrometeorological time series. Part 1. The normal model revisited[J]. Journal of Hydrology, 2000,235(3):221-241.
[6] KOTZ S, 吳喜之.現代貝葉斯統(tǒng)計學[M].北京:中國統(tǒng)計出版社,2000.
[7] LU W, SHI D. A new compounding life distribution: the Weibull-Poisson distribution[J]. Journal of Applied Statistics, 2012,39(1):21-38.
[8] BARRETO-SOUZA W, CRIBARI-NETO F. A generalization of the exponential-poisson distribution[J]. Statistics & Probability Letters, 2009,79(24):2493-2500.
[9] POISSON S D. Recherches sur la probabilité des jugements en matiére criminelle et en matière civile, précédées des règles générales du calcul des probabilités[M]. Bachelier, 1837.
[10] BARBOUR A D, HOLST L, JANSON S. Poisson approximation[M]. Oxford: Clarendon Press, 1992.
[11] 劉銀萍,馬曉悅,趙志文.缺失數據場合泊松分布參數的貝葉斯估計[J].吉林師范大學學報:自然科學版,2012,33(3):13-15.
[12] 劉銀萍,宋立新.II型截尾情形下泊松分布參數的估計[J].吉林大學學報:理學版,2007,45(6):941-944.
[13] 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.
[14] 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.
[15] 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.
[16] 陳怡南,葉爾驊.IIRCT下對數正態(tài)和正態(tài)分布參數的MLE[J].南京航空航天大學學報,1996,28(3):376-383.
[17] 陳怡南,葉爾驊.帶有不完全信息隨機截尾試驗下Weibull分布參數的MLE[J].數理統(tǒng)計與應用概率,1996,11(4):353-363.
[18] 楊紀龍,葉爾驊.帶有不完全信息隨機截尾試驗下Weibull分布數參的MLE的相合性及漸近正態(tài)性[J].應用概率統(tǒng)計,2000,16(1):9-19.
[19] 張曉琴,張虎明.帶有不完全信息隨機截尾試驗下Weibull分布參數的MLE的重對數律[J].應用概率統(tǒng)計,2002,18(1):101-107.
[20] 宋毅君,李補喜,李濟洪.帶有不完全信息隨機截尾試驗下最大似然估計的重對數律[J].應用概率統(tǒng)計,2009,25(2):113-125.
BayesEstimationofParameterofPoissonDistributionwithSingleChangePointforRandomCensoringTestModelwithIncompleteInformation
HE Chao-bing1, LIU Hua-wen2
(1. School of Mathematics and Statistics, Anyang Normal University, Anyang 455000, China; 2. School of Mathematics, Shandong University, Jinan 250100, China)
This paper firstly obtained the complete-data likelihood function of Poisson distribution for IIRCT after adding data, then studied the full conditional distributions of change-point position and other parameters, and estimated the parameters by MCMC method of Gibbs sampling together with Metropolis-Hastings algorithm. Finally random simulation tests were conducted, and the results showed that Bayes estimations of the parameters were fairly accurate.
complete-data likelihood function; full conditional distribution; MCMC method; Gibbs sampling; Metropolis-Hastings algorithm
2013-09-22
國家自然科學基金(61174099);河南省教育廳自然科學基金資助項目(2011B110001).
何朝兵(1975-),男,河南周口人,講師,碩士,主要從事概率統(tǒng)計的研究.
何朝兵,劉華文.IIRCT下泊松分布參數單變點的貝葉斯估計[J].安徽師范大學學報:自然科學版,2014,37(4):335-338.
O213.2;O212.8
A
1001-2443(2014)04-0335-04