一種基于無信息先驗(yàn)下泊松分布變點(diǎn)的貝葉斯估計(jì)

許 婷，吳有富，張英雪

（1.貴州民族大學(xué)數(shù)據(jù)科學(xué)與信息工程學(xué)院，貴州 貴陽 550025；2.貴州交通職業(yè)技術(shù)學(xué)院 貴州 貴陽 550025）

變點(diǎn)問題的貝葉斯估計(jì)注重的是先驗(yàn)分布的選取，不同先驗(yàn)分布的選取會(huì)對(duì)估計(jì)精度有不同的影響。對(duì)于泊松分布參數(shù)變點(diǎn)的研究，張夢(mèng)琇[1]選取先驗(yàn)分布為共軛先驗(yàn)，并通過二分法實(shí)現(xiàn)了對(duì)泊松分布多變點(diǎn)模型的參數(shù)估計(jì)。范元靜[2]在選擇先驗(yàn)分布為共軛先驗(yàn)伽馬分布的基礎(chǔ)上，運(yùn)用RJMCMC算法對(duì)泊松分布的變點(diǎn)參數(shù)進(jìn)行數(shù)量及位置的估計(jì)。何朝兵[3]選取共軛先驗(yàn)分布作為泊松分布的先驗(yàn)分布，并利用貝葉斯方法對(duì)IIRCT下單變點(diǎn)泊松分布的參數(shù)進(jìn)行估計(jì)。鐘穎和田茂[4]再選取泊松分布的先驗(yàn)分布為共軛先驗(yàn)分布，并使用貝葉斯方法對(duì)單變點(diǎn)泊松分布進(jìn)行參數(shù)估計(jì)。胡興[5]使用Fisher信息陣確定了泊松分布的先驗(yàn)分布，并通過MCMC方法對(duì)推導(dǎo)的后驗(yàn)分布進(jìn)行模擬研究。對(duì)于泊松分布變點(diǎn)問題的研究，人們幾乎都是選擇泊松分布的共軛先驗(yàn)分布伽馬分布作為其先驗(yàn)分布對(duì)其變點(diǎn)問題進(jìn)行研究，在確定無信息先驗(yàn)分布的方法下對(duì)于泊松分布先驗(yàn)分布的選取也僅僅是利用Fisher信息陣確定了其先驗(yàn)分布，而對(duì)于其他無信息先驗(yàn)下的泊松分布變點(diǎn)問題的估計(jì)目前尚未報(bào)道。文章基于此，在Zhang等[6]討論的確定無信息先驗(yàn)分布方法的基礎(chǔ)上，基于“變換不變性”原則推導(dǎo)出Poisson分布參數(shù)的無信息先驗(yàn)分布，在參數(shù)的全條件后驗(yàn)分布下對(duì)Poisson分布尺度參數(shù)變點(diǎn)模型的參數(shù)估計(jì)問題進(jìn)行研究。

1 Poisson分布單變點(diǎn)模型

2 極大似然估計(jì)

3 貝葉斯估計(jì)

貝葉斯估計(jì)作為常用的對(duì)參數(shù)進(jìn)行統(tǒng)計(jì)推斷的方法，在使用貝葉斯方法對(duì)變點(diǎn)問題進(jìn)行研究時(shí)，是將模型中把變點(diǎn)涵蓋在內(nèi)的所有參數(shù)都視為隨機(jī)變量。通過引入先驗(yàn)分布，先驗(yàn)分布是反映在抽樣前對(duì)參數(shù)的認(rèn)識(shí)，對(duì)參數(shù)設(shè)定不同的先驗(yàn)分布，根據(jù)設(shè)定的先驗(yàn)分布和已知樣本服從的分布計(jì)算出每個(gè)參數(shù)的后驗(yàn)分布，后驗(yàn)分布是反映在抽樣后對(duì)參數(shù)的認(rèn)識(shí)，根據(jù)每個(gè)參數(shù)全條件下的后驗(yàn)分布對(duì)參數(shù)進(jìn)行估計(jì)。由于不同的先驗(yàn)分布會(huì)對(duì)參數(shù)估計(jì)結(jié)果有不同的影響，所以在選取先驗(yàn)分布的時(shí)候很重要。文章選用的是無信息先驗(yàn)分布，在此基礎(chǔ)上對(duì)泊松分布的變點(diǎn)進(jìn)行貝葉斯估計(jì)。

3.1 不變性原則下泊松分布的無信息先驗(yàn)分布

3.2 無信息先驗(yàn)下泊松分布單變點(diǎn)參數(shù)的貝葉斯估計(jì)

首先確各參數(shù)的先驗(yàn)分布：

4 數(shù)值仿真

考慮下面的變點(diǎn)模型：

（1）極大似然估計(jì)法

根據(jù)推導(dǎo)出的Poisson分布的對(duì)數(shù)似然函數(shù)式（4）和各參數(shù)的極大似然估計(jì)式（8）對(duì)參數(shù)，k進(jìn)行估計(jì)，使用R軟件進(jìn)行100次模擬的結(jié)果如表1所示：

表1 參數(shù)k12的極大似然估計(jì)

表1 參數(shù)k12的極大似然估計(jì)

參數(shù) 真值 估計(jì)值 相對(duì)誤差1 10 11.26 0.126 20 21.51 0.076 k 40 42.73 0.068 2

（2）貝葉斯估計(jì)

表2 參數(shù)k12的貝葉斯估計(jì)

表2 參數(shù)k12的貝葉斯估計(jì)

參數(shù) 真值Gibbs樣本1 10 10.01 0.0010.000139.672 9.99510.346 20000 2 20 19.410.0295 0.0002919.03 19.41 19.79 20000 k 40 39.980.0005 0.00057 40 39.98 40.00 20000估計(jì)值相對(duì)誤差MC誤差2.5%分位數(shù)中位數(shù)97.5%分位數(shù)

圖1 1迭代軌跡圖

圖2 2迭代軌跡圖

圖3 k迭代軌跡圖

圖4 1后驗(yàn)分布的核密度估計(jì)圖

圖5 2后驗(yàn)分布的核密度估計(jì)圖

圖6 k后驗(yàn)分布的核密度估計(jì)圖

雖然從圖1到圖3中可以看出各參數(shù)產(chǎn)生的馬爾科夫鏈已經(jīng)是平穩(wěn)狀態(tài)了，但是為了判定生成的Gibbs樣本是不是收斂的，在用MCMC方法進(jìn)行模擬的時(shí)候，會(huì)同時(shí)產(chǎn)生多個(gè)迭代鏈，若生成的多個(gè)迭代鏈逐漸穩(wěn)定并且趨于重合，那么可以說明Gibbs抽樣是收斂的。由于文章主要是是為了判定變點(diǎn)的位置，因此在模擬的過程中，固定了變點(diǎn)k的位置，用不同的的值分別進(jìn)行20000次迭代，其兩組初始值分別為，產(chǎn)生的迭代鏈軌跡圖如圖 7—圖 9所示。其中圖 7是下的迭代鏈，從圖7中可以看出，其抽樣基本上都在變點(diǎn)附近波動(dòng)，且波動(dòng)的幅度較小，說明估計(jì)的效果較好；圖8是下的迭代鏈，同樣可以看出，其抽樣基本在變點(diǎn)附近波動(dòng)，抽樣的波動(dòng)程度同樣較小，估計(jì)的效果也較好；圖9是兩組初始值下同時(shí)產(chǎn)生的迭代鏈，從圖中可以看出，產(chǎn)生的兩個(gè)馬爾科夫鏈都是穩(wěn)定且趨于重合的，表明Gibbs樣本是收斂的。

圖7 k第一組初始值下的迭代鏈

圖8 k第二組初始值下的迭代鏈

圖9 k兩組初始值下的迭代鏈

表3 不同樣本量下參數(shù)k12的貝葉斯估計(jì)

表3 不同樣本量下參數(shù)k12的貝葉斯估計(jì)

總樣本量n 100 200 500 1000 2000 5000 1樣本量 40 110 300 400 1100 3000 2樣本量 60 90 200 600 900 2000 1估計(jì)值 10.01 10.30 10.16 9.94 10.05 2估計(jì)值 19.41 20.25 19.94 20.09 19.81 19.85 k估計(jì)值 39.98 111.35300.07 399.861099.64 2999.05

文章還給出了總樣本量 n=100，200，500，1000，2000，5000時(shí)的三個(gè)參數(shù)的估計(jì)值，如表2所示，其中樣本量的值就是變點(diǎn)的位置。從表中可以看出，在給定的的情況下，不管n是小樣本還是大樣本，每個(gè)參數(shù)的估計(jì)值都與真實(shí)值非常接近。且不管樣本量增大到多少，其估計(jì)精度都還是挺高的。因此從模擬的結(jié)果來看，參數(shù)的估計(jì)值與真實(shí)值都是很接近的，表明MCMC的模擬效果是比較好的。

5 實(shí)證分析

我們將我們的方法應(yīng)用于Jarrett[8]給出的英國(guó)煤災(zāi)事故數(shù)據(jù)集，他修正并擴(kuò)展了Maguire[9]等人給出的數(shù)據(jù)集。

表4 英國(guó)年度煤礦災(zāi)害數(shù)據(jù)，1851-1962年

Carlin等人[10]用共軛先驗(yàn)得出的礦難發(fā)生變化的點(diǎn)是在41，也就是在1890年底至1892年初之間礦難次數(shù)開始減少，用我們的方法得到的結(jié)果與其是差不多的，Raftery和Akman[11]用模糊先驗(yàn)也給出了相似的結(jié)果。

從圖10中可以看出，發(fā)生突變的位置是40附近，也就是位置41的后驗(yàn)概率最大，其次是位置40和位置39，對(duì)應(yīng)下來也就是1890年附近，也即是意味著這種變化最有可能發(fā)生在1889年末至1892年初之間。且從圖11中可以看到，發(fā)生變化之前平均每年要發(fā)生3次左右的礦難，而在變化點(diǎn)之后每年發(fā)生的礦難次數(shù)平均是1次左右。

圖10 變點(diǎn)位置后驗(yàn)分布密度曲線

圖11 變點(diǎn)前年均礦難次數(shù)密度曲線

圖12 變點(diǎn)后年均礦難次數(shù)密度曲線

用MCMC抽樣方法對(duì)礦難災(zāi)害的估計(jì)結(jié)果如下：

表5 參數(shù)估計(jì)

從表5中可以看出，MC誤差較小，說明MCMC方法的收斂性較好。將用文章的方法得出的估計(jì)值與Carlin等人研究的結(jié)果相比，相對(duì)誤差都小于0.1，說明估計(jì)效果較好。從估計(jì)結(jié)果來看，英國(guó)煤礦災(zāi)害變點(diǎn)位置之后發(fā)生災(zāi)難的次數(shù)很大程度降低的原因可能是與19世紀(jì)英國(guó)頒布的一系列有關(guān)采礦安全等法案有關(guān)，實(shí)施之后對(duì)采礦安全有了一定程度上的影響。

6 結(jié)論

文章主要在已有的推導(dǎo)無信息先驗(yàn)分布的方法下，基于不變性理論下推導(dǎo)出Poisson分布的無信息先驗(yàn)分布，利用貝葉斯方法實(shí)現(xiàn)對(duì)參數(shù)的估計(jì)，并和極大似然方法做比較。

根據(jù)表1和表2中的估計(jì)值可以看出，兩種方法都能估計(jì)出變點(diǎn)位置參數(shù) k的位置和尺度參數(shù)的值。但是通過比較兩種方法的相對(duì)誤差可以知道極大似然估計(jì)的精確度沒有貝葉斯估計(jì)的高，極大似然方法得出的估計(jì)值與給定的真實(shí)值相差較大，而貝葉斯方法得出的估計(jì)值與給定的真實(shí)值相差較小。因此可以知道用貝葉斯方法估計(jì)變點(diǎn)位置參數(shù)k的效果更好一些。根據(jù)表3中不同樣本量下的估計(jì)值，k可以知道，總樣本量n增加時(shí)，各個(gè)參數(shù)的估計(jì)值依然與真實(shí)值很接近，說明貝葉斯方法在小樣本和大樣本下都可以用來估計(jì)參數(shù)值，且估計(jì)效果比極大似然方法好。

并用文章使用的先驗(yàn)應(yīng)用于英國(guó)煤災(zāi)事故數(shù)據(jù)集，通過實(shí)證分析的結(jié)果可以看出，與前人的研究結(jié)果是基本一致的。綜上所述，用貝葉斯方法對(duì)變點(diǎn)位置參數(shù)k和尺度參數(shù)的估計(jì)效果都是非常不錯(cuò)的，因此對(duì)于推導(dǎo)出來的Poisson分布的尺度參數(shù)的無信息先驗(yàn)為是可行且精度較高的，使用此先驗(yàn)分布對(duì)Poisson分布的變點(diǎn)問題進(jìn)行研究是有效的。