何朝兵，劉華文

(1.安陽(yáng)師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，河南 安陽(yáng) 455000；2.山東大學(xué) 數(shù)學(xué)學(xué)院，山東 濟(jì)南 2501000)

IIRCT下泊松分布參數(shù)單變點(diǎn)的貝葉斯估計(jì)

何朝兵1，劉華文2

(1.安陽(yáng)師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，河南 安陽(yáng) 455000；2.山東大學(xué) 數(shù)學(xué)學(xué)院，山東 濟(jì)南 2501000)

首先通過(guò)添加數(shù)據(jù)得到了帶有不完全信息隨機(jī)截尾試驗(yàn)下泊松分布的完全數(shù)據(jù)似然函數(shù)，然后研究了變點(diǎn)位置和其它參數(shù)的滿條件分布，接著利用Gibbs抽樣與Metropolis-Hastings算法相結(jié)合的MCMC方法對(duì)參數(shù)進(jìn)行了估計(jì)，最后進(jìn)行了隨機(jī)模擬，試驗(yàn)結(jié)果表明參數(shù)貝葉斯估計(jì)的精度較高.

完全數(shù)據(jù)似然函數(shù)；滿條件分布；MCMC方法；Gibbs抽樣；Metropolis-Hastings算法

變點(diǎn)分析研究始于20世紀(jì)50年代，自20世紀(jì)70年代初，這一問(wèn)題引起了一些編譯學(xué)家的重視，發(fā)展了一些變點(diǎn)分析方法[1-3]，如最小二乘法、極大似然法、貝葉斯法和非參數(shù)方法等.隨著統(tǒng)計(jì)計(jì)算技術(shù)的快速發(fā)展，貝葉斯方法的應(yīng)用越來(lái)越廣泛，特別是在水文統(tǒng)計(jì)[4,5]和計(jì)量經(jīng)濟(jì)[6]等領(lǐng)域.而貝葉斯計(jì)算方法中的Markov chain Monte Carlo(MCMC)方法，使變點(diǎn)分析中貝葉斯方法的實(shí)際操作變得非常方便.泊松分布是一類應(yīng)用廣泛的離散型壽命分布，它可以作為多種不同類型試驗(yàn)的模型，泊松分布的另一個(gè)應(yīng)用是研究空間分布，例如湖泊中魚(yú)群數(shù)量的分布.文獻(xiàn)[7-12]對(duì)泊松分布進(jìn)行了深入研究，其中文獻(xiàn)[12]研究了II型截尾情形下泊松分布參數(shù)的極大似然估計(jì).近些年來(lái)，關(guān)于隨機(jī)截尾試驗(yàn)的研究比較多，帶有不完全信息隨機(jī)截尾試驗(yàn)(random censoring test model with incomplete information)，簡(jiǎn)稱IIRCT，由文獻(xiàn)[13]首次涉及，接著文獻(xiàn)[14-20]深入研究了IIRCT下連續(xù)型分布的參數(shù)估計(jì).實(shí)際上，文獻(xiàn)[12]研究的II型截尾試驗(yàn)是特殊的IIRCT.而對(duì)IIRCT下壽命分布參數(shù)的變點(diǎn)問(wèn)題，卻很少有文獻(xiàn)研究.下文主要研究了IIRCT下泊松分布參數(shù)單變點(diǎn)的貝葉斯估計(jì).首先通過(guò)添加數(shù)據(jù)得到了IIRCT下泊松分布的完全數(shù)據(jù)似然函數(shù)，然后研究了變點(diǎn)位置和其它參數(shù)的滿條件分布，接著利用Gibbs抽樣與Metropilis-Hastings算法相結(jié)合的MCMC方法得到了參數(shù)的Gibbs樣本，把Gibbs樣本的均值作為各參數(shù)的貝葉斯估計(jì)，最后進(jìn)行了隨機(jī)模擬，試驗(yàn)結(jié)果表明各參數(shù)貝葉斯估計(jì)的精度都較高.

1 離散型壽命IIRCT

設(shè)產(chǎn)品壽命X1，X2,L是獨(dú)立同分布的離散型隨機(jī)變量序列，其分布函數(shù)為F(x,λ)=P(Xi≤x)，分布為f(x,λ)，λ是未知參數(shù).又設(shè)Y1,Y2,L是獨(dú)立的取非負(fù)整數(shù)的離散型隨機(jī)變量序列，分布函數(shù)分別為G1(y),G2(y),L.分布律分別為g1(y),g2(y),L,且gi(y)與參數(shù)λ無(wú)關(guān).{Xi}與{Yi}獨(dú)立.

為估計(jì)參數(shù)λ,n個(gè)樣品的觀察數(shù)據(jù){Zi,1≤i≤n}如下：

(1)當(dāng)Xi≤Yi時(shí)，有兩種情況：Xi以概率ai立即顯示，此時(shí)取Zi=Xi;Xi以概率1-ai不被顯示，此時(shí)取Zi=Yi,稱ai為失效顯示概率.

(2)當(dāng)Xi>Yi時(shí)，取Zi=Yi.

引入隨機(jī)變量αi,βi,i=1,2,L,n.

若Xi≤Yi,αi=1；若Xi>Yi,αi=0.

若Xi≤Yi，且未被顯示時(shí)，βi=0;其它情況，βi=1.

綜上所述，有

設(shè)Zi的觀察值為zi，則基于數(shù)據(jù){(zi,αi,βi),1≤i≤n}的似然函數(shù)為

2 IIRCT下泊松分布參數(shù)的完全數(shù)據(jù)似然函數(shù)

若X的分布律為P(X=k)=λke-λ/k!,k=0,1,L,則稱X服從參數(shù)為λ的泊松分布，記為X～P(λ).假設(shè)IIRCT下產(chǎn)品壽命X～P(λ).

令α表示αi組成的向量，β表示βi組成的向量，z表示z組成的向量.

泊松分布P(λ)基于數(shù)據(jù){(zi,α,βi),1≤i≤n}的似然函數(shù)為

從上式可以看出，基于截尾數(shù)據(jù)的似然函數(shù)比較復(fù)雜，為了方便進(jìn)行參數(shù)估計(jì)，下面添加缺損的Xi的值，以獲得完全數(shù)據(jù)的似然函數(shù).具體如下：

當(dāng)αi=1,βi=0，即第i個(gè)產(chǎn)品失效未被顯示時(shí)，添加數(shù)據(jù)Z1i=Xi=z1i.

在Xi≤zi的條件下，Z1i的條件分布律為

當(dāng)αi=0，即Xi>Yi時(shí)，添加數(shù)據(jù)Z2i=Xi=z2i.

在Xi>zi的條件下，Z2i的條件分布律為

令u表示z1i組成的向量，v表示z2i組成的向量，則完全數(shù)據(jù)似然函數(shù)為

3 IIRCT下泊松分布參數(shù)單變點(diǎn)的貝葉斯估計(jì)

泊松分布參數(shù)單變點(diǎn)模型如下：

其中λ1,λ2兩兩不相等，1≤k≤n-1.

下面對(duì)變點(diǎn)位置k及參數(shù)λ1,λ2進(jìn)行貝葉斯估計(jì).

記θ=(k,λ1,λ2)，則IIRCT下泊松分布參數(shù)單變點(diǎn)問(wèn)題的似然函數(shù)為

下面確定各參數(shù)的先驗(yàn)分布.

(2)取λi的先驗(yàn)分布為共軛先驗(yàn)分布伽瑪分布Ga(bi,ci),bi,ci已知，即

假設(shè)k,λ,λ2獨(dú)立.則

L(θ|z,u,v,α,β)∝π(k)π(λ1)π(λ2)L(z,u,v,α,β|θ)

當(dāng)αi=1,βi=0時(shí)，

其中z-1i={z1j:j≠i}.

當(dāng)αi=0時(shí)，

其中z-2i={z2j:j≠i}.

為了書(shū)寫(xiě)方便，把滿條件分布中的“條件”用“·”代替，例如

π(λ1|k,λ2,z,u,v,α,β)簡(jiǎn)記為π(λ1|·).

各參數(shù)的滿條件分布為

由于得到了各參數(shù)的滿條件分布，下面利用MCMC方法獲得各參數(shù)后驗(yàn)分布的平穩(wěn)分布.參數(shù)z1i,z2i的滿條件分布是截?cái)嗖此煞植?，可以利用逆變換法隨機(jī)抽樣，λ1,λ2的滿條件分布是伽瑪分布，這4個(gè)分布都可以采用Gibbs抽樣；但是k的滿條件分布不是標(biāo)準(zhǔn)分布，進(jìn)行Gibbs抽樣比較困難，可以利用Metropolis-Hastings算法進(jìn)行抽樣，此時(shí)選取建議分布為取值1，2，L，n-1的離散型均勻分布.

下面寫(xiě)出MCMC方法的具體步驟：

從1，2，L，n-1中任取一個(gè)k′，產(chǎn)生一個(gè)隨機(jī)數(shù)u0，若u0≤r(k(t-1),k′),則k(t)=k′，否則k(t)=k(t-1).

4 隨機(jī)模擬

基于上面的討論，下面進(jìn)行隨機(jī)模擬試驗(yàn).取受試樣品的個(gè)數(shù)n=200，樣品壽命Xi:P(10),i=1,2,L,50;Xi:P(30),i=51,52,L,200.則參數(shù)(k,λ1,λ2)的真實(shí)值為(50，10，30).取λ1,λ2的先驗(yàn)分布分別為Ga(7,0.6),Ga(35,1.2),截尾變量Yi～P(32)，失效顯示概率a=0.8.

利用前面推導(dǎo)出的各個(gè)參數(shù)的滿條件分布使用R軟件進(jìn)行MCMC模擬.在模擬運(yùn)行過(guò)程中，先進(jìn)行10000次Gibbs預(yù)迭代，以確保參數(shù)的收斂性，然后丟棄最初的預(yù)迭代，再進(jìn)行10000次Gibbs迭代.迭代從第10001次開(kāi)始至第20000次的R程序的運(yùn)行結(jié)果見(jiàn)表1.

表1 參數(shù)k,λ1,λ2的貝葉斯估計(jì)

最后，進(jìn)行模擬結(jié)果分析.由表1可以看出各參數(shù)的估計(jì)值與真值的相對(duì)誤差都小于3%，精度較高，效果較好；各參數(shù)的MC誤差較小，Gibbs迭代值的波動(dòng)較小，Gibbs迭代比較穩(wěn)定.綜上分析，可以看出通過(guò)MCMC模擬所產(chǎn)生的效果較好.

注：編寫(xiě)R程序時(shí)用到的函數(shù)主要有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] 陳希孺.變點(diǎn)統(tǒng)計(jì)分析簡(jiǎn)介[J].數(shù)理統(tǒng)計(jì)與管理，1991，10(1)：55-59.

[3] 項(xiàng)靜怡，史久恩.非線性系統(tǒng)中數(shù)據(jù)處理的統(tǒng)計(jì)方法[M].北京：科學(xué)出版社，1997.

[4] 熊立華，周芬.水文時(shí)間序列變點(diǎn)分析的貝葉斯方法[J].水電能源科學(xué)，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, 吳喜之.現(xiàn)代貝葉斯統(tǒng)計(jì)學(xué)[M].北京：中國(guó)統(tǒng)計(jì)出版社，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] 劉銀萍，馬曉悅，趙志文.缺失數(shù)據(jù)場(chǎng)合泊松分布參數(shù)的貝葉斯估計(jì)[J].吉林師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2012，33(3)：13-15.

[12] 劉銀萍，宋立新.II型截尾情形下泊松分布參數(shù)的估計(jì)[J].吉林大學(xué)學(xué)報(bào)：理學(xué)版，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下對(duì)數(shù)正態(tài)和正態(tài)分布參數(shù)的MLE[J].南京航空航天大學(xué)學(xué)報(bào)，1996，28(3)：376-383.

[17] 陳怡南，葉爾驊.帶有不完全信息隨機(jī)截尾試驗(yàn)下Weibull分布參數(shù)的MLE[J].數(shù)理統(tǒng)計(jì)與應(yīng)用概率，1996，11(4)：353-363.

[18] 楊紀(jì)龍，葉爾驊.帶有不完全信息隨機(jī)截尾試驗(yàn)下Weibull分布數(shù)參的MLE的相合性及漸近正態(tài)性[J].應(yīng)用概率統(tǒng)計(jì)，2000，16(1)：9-19.

[19] 張曉琴，張虎明.帶有不完全信息隨機(jī)截尾試驗(yàn)下Weibull分布參數(shù)的MLE的重對(duì)數(shù)律[J].應(yīng)用概率統(tǒng)計(jì)，2002，18(1)：101-107.

[20] 宋毅君，李補(bǔ)喜，李濟(jì)洪.帶有不完全信息隨機(jī)截尾試驗(yàn)下最大似然估計(jì)的重對(duì)數(shù)律[J].應(yīng)用概率統(tǒng)計(jì)，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

國(guó)家自然科學(xué)基金(61174099)；河南省教育廳自然科學(xué)基金資助項(xiàng)目(2011B110001).

何朝兵(1975-)，男，河南周口人，講師，碩士，主要從事概率統(tǒng)計(jì)的研究.

何朝兵，劉華文.IIRCT下泊松分布參數(shù)單變點(diǎn)的貝葉斯估計(jì)[J].安徽師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2014，37(4)：335-338.

O213.2；O212.8

A

1001-2443(2014)04-0335-04