左截?cái)嘤覄h失數(shù)據(jù)下Pareto分布形狀參數(shù)多變點(diǎn)的貝葉斯估計(jì)

何朝兵　杜保建　劉華文

摘 要 通過(guò)添加缺損的壽命變量數(shù)據(jù)得到了左截?cái)嘤覄h失數(shù)據(jù)下Pareto分布相對(duì)簡(jiǎn)單的似然函數(shù)，給出了形狀參數(shù)變點(diǎn)位置和其他參數(shù)的滿條件分布.利用MCMC方法對(duì)參數(shù)的滿條件分布進(jìn)行了抽樣，把Gibbs樣本的均值作為參數(shù)的貝葉斯估計(jì).隨機(jī)模擬的結(jié)果表明各參數(shù)貝葉斯估計(jì)的精度都較高.

關(guān)鍵詞 似然函數(shù)；滿條件分布；MCMC方法；Gibbs抽樣；Metropolis-Hastings算法

中圖分類號(hào) O213.2； O212.8 文獻(xiàn)標(biāo)識(shí)碼 A 文章編號(hào) 1000-2537（2015）03-0080-05

近年來(lái)變點(diǎn)問(wèn)題成為統(tǒng)計(jì)學(xué)中比較熱門的研究方向，它廣泛應(yīng)用于金融、經(jīng)濟(jì)和地震預(yù)測(cè)等領(lǐng)域.目前變點(diǎn)分析方法主要有極大似然法、貝葉斯法和非參數(shù)方法等.關(guān)于變點(diǎn)問(wèn)題的研究可參看文獻(xiàn)[1-4].而貝葉斯計(jì)算方法中的MCMC方法，使變點(diǎn)分析變得非常方便.Pareto分布在生存分析和可靠性理論等方面仍然具有很多應(yīng)用價(jià)值.關(guān)于Pareto分布的研究可參看文獻(xiàn)[5-7].當(dāng)進(jìn)行壽命試驗(yàn)時(shí)，經(jīng)常會(huì)出現(xiàn)左截?cái)嘤覄h失數(shù)據(jù)，對(duì)此模型的研究可參看文獻(xiàn)[8-10].關(guān)于左截?cái)嗷蛴覄h失數(shù)據(jù)下變點(diǎn)問(wèn)題的研究有一些成果，可參看文獻(xiàn)[11-13]，但對(duì)數(shù)據(jù)既左截?cái)嘤钟覄h失情形下變點(diǎn)問(wèn)題的研究還不多見(jiàn).本文主要利用MCMC方法研究左截?cái)嘤覄h失數(shù)據(jù)下Pareto分布形狀參數(shù)多變點(diǎn)的貝葉斯估計(jì).

1 左截?cái)嘤覄h失數(shù)據(jù)試驗(yàn)?zāi)Ｐ?/p>

設(shè)（X，Y，T）是一連續(xù)型隨機(jī)變量，壽命變量X的分布函數(shù)為F（x；θ）=P（X≤x），密度函數(shù)為f（x；θ），這里θ是未知參數(shù)；Y是一右刪失隨機(jī)變量，分布函數(shù)為G（y），密度函數(shù)為g（y）；T是一左截?cái)嚯S機(jī)變量，分布函數(shù)為H（t），密度函數(shù)為h（t），且Y，T的分布與參

數(shù)θ無(wú)關(guān).假定X，Y，T是相互獨(dú)立非負(fù)隨機(jī)變量.對(duì)于n個(gè)受試樣品（產(chǎn)品壽命），左截?cái)嘤覄h失數(shù)據(jù)的試驗(yàn)?zāi)Ｐ褪莾H在Zi≥Ti時(shí)得到觀察數(shù)據(jù)（Zi，Ti，δi），而在Zi

2 左截?cái)嘤覄h失數(shù)據(jù)下Pareto分布的似然函數(shù)

若X的密度函數(shù)為f（x）=θbθx-（θ+1），x≥b>0，θ>0，則稱X服從尺度參數(shù)為b，形狀參數(shù)為θ的Pareto分布，記為X～Pa（b，θ），易知X的分布函數(shù)為F（x）=1-bθx-θ，x≥b.

由（1）式易知，左截?cái)嘤覄h失數(shù)據(jù)下Pareto分布的似然函數(shù)比較復(fù)雜.

下面添加部分缺損的Xi的值，以獲得較簡(jiǎn)單的似然函數(shù).具體如下：

可以利用篩選法隨機(jī)產(chǎn)生Z1i的值z(mì)1i.z1i抽樣的具體步驟如下：

3 左截?cái)嘤覄h失數(shù)據(jù)下Pareto分布形狀參數(shù)多變點(diǎn)的貝葉斯估計(jì)

Pareto分布形狀參數(shù)多變點(diǎn)模型如下：

Xi～Pa（b，θ1），i=1，2，…，k1，Pa（b，θ2），i=k1+1，…，k2，Pa（b，θ3），i=k2+1，…，n， （3）

其中θ1，θ2，θ3兩兩不相等，1≤k1

下面利用貝葉斯方法對(duì)對(duì)變點(diǎn)位置k1，k2及參數(shù)θ1，θ2，θ3進(jìn)行估計(jì).

令D1={1，2，…，k1}，D2={k1+1，…，k2}，D3={k2+1，…，n}.

記β=（k1，k2，θ1，θ2，θ3），由（2）式和（3）式可得此變點(diǎn)問(wèn)題的似然函數(shù)為

其中nm1=n1（Dm），nm2=n2（Dm），sm=s（Dm），m=1，2，3.

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

（1） 對(duì)于（k1，k2）取無(wú)信息先驗(yàn)分布：π（k1，k2）=1/C2n-1=2[（n-1）（n-2）]-1，1≤k1

（2） 取θi的先驗(yàn)分布為伽瑪分布Ga（ci，di），ci，di已知，即

假設(shè)（k1，k2），θ1，θ2，θ3相互獨(dú)立，則

4 隨機(jī)模擬

下面進(jìn)行隨機(jī)模擬試驗(yàn).取受試樣品的個(gè)數(shù)n=200.

假設(shè)Pareto分布的尺度參數(shù)8是已知的，取θ1，θ2，θ3的先驗(yàn)分布分別為Ga（3.5，6），Ga（20，1.6），Ga（17，2.6）；右刪失變量Yi～Pa（10，2），左截?cái)嘧兞縏i～Pa（6，10）.則（k1，k2，θ1，θ2，θ3）的真實(shí)值為（80，150，0.7，13，6）.

使用R軟件進(jìn)行MCMC模擬，主要對(duì)參數(shù)k1，k2進(jìn)行估計(jì)分析.在模擬過(guò)程中，先進(jìn)行10 000次Gibbs預(yù)迭代，以確保抽樣的收斂性，然后丟棄最初的預(yù)迭代，再進(jìn)行10 000次Gibbs迭代.迭代從第10 001次開(kāi)始至第20 000次的R程序的運(yùn)行結(jié)果見(jiàn)表1.

在模型的分析過(guò)程中，MCMC的收斂性診斷很重要，模擬時(shí)可以對(duì)參數(shù)進(jìn)行多層鏈?zhǔn)降治?，即輸入多組初始值，形成多層迭代鏈，當(dāng)抽樣收斂時(shí)，迭代圖形重合.在模擬過(guò)程中，輸入兩組初始值分別進(jìn)行10 000次迭代，變點(diǎn)位置參數(shù)的多層迭代鏈軌跡見(jiàn)圖3和圖4.

最后，進(jìn)行模擬結(jié)果分析.首先，由表1可以看出位置參數(shù)的貝葉斯估計(jì)與真值的相對(duì)誤差不超過(guò)4%，其他參數(shù)估計(jì)的相對(duì)誤差不超過(guò)9%，整體上估計(jì)的精度較高，效果較好；其次，要判斷所產(chǎn)生的馬爾科夫鏈?zhǔn)欠袷諗?，由圖3和圖4可以發(fā)現(xiàn)，變點(diǎn)位置參數(shù)的兩條迭代鏈都趨于重合，收斂性較好.綜上分析，可以看出通過(guò)MCMC方法模擬所產(chǎn)生的效果較好.

參考文獻(xiàn)：

[1] PAGE E S. Continuous inspection schemes[J]. Biometrika， 1954，41（1）：100-115.

[2] CSRG M， HORVTH L. Limit theorems in change-point analysis[M]. New York： Wiley， 1997.

[3] CHERNOFF H， ZACKS S. Estimating the current mean of a normal distribution which is subjected to changes in time[J]. Ann Math Stat， 1964，35（3）：999-1018.

[4] FEARNHEAD P. Exact and efficient Bayesian inference for multiple changepoint problems[J]. Stat Comput， 2006，16（2）：203-213.

[5] ARNOLD B C. Pareto distribution[M]. New York： John Wiley & Sons， Inc， 1985.

[6] CASTILLO J， DAOUDI J. Estimation of the generalized Pareto distribution[J]. Stat Probab Lett， 2009，79（5）：684-688.

[7] KAMAL M， ZARRIN S， ISLAM A U. Accelerated life testing design using geometric process for Pareto distribution[J]. Int J Adv Stat Probab， 2013，1（2）：25-31.

[8] BALAKRISHNAN N， MITRA D. Likelihood inference for lognormal data with left truncation and right censoring with an illustration[J]. J Stat Plann Inference， 2011，141（11）：3536-3553.

[9] COSSLETT S R. Efficient semiparametric estimation of censored and truncated regressions via a smoothed self-consistency equation[J]. Econometrica， 2004，72（4）：1277-1293.

[10] GROSS S T， LAI T L. Nonparametric estimation and regression analysis with left-truncated and right-censored data[J]. J Am Stat Assoc， 1996，91（435）：1166-1180.

[11] ANTONIADIS A， GRGOIRE G. Nonparametric estimation in change point hazard rate models for censored data： A counting process approach[J]. J Nonparametr Stat， 1993，3（2）：135-154.

[12] ANTONIADIS A， GIJBELS I， MACGIBBON B. Nonparametric estimation for the location of a change-point in an otherwise smooth hazard function under random censoring[J]. Scand J Stat， 2000，27（3）：501-519.

[13] GIJBELS I， GRLER . Estimation of a change point in a hazard function based on censored data[J]. Lifetime Data Anal， 2003，9（4）：395-411.

（編輯 胡文杰）