基于RJMCMC算法的Gamma分布形狀參數(shù)多變點(diǎn)檢測(cè)

程慧慧, 許淑月

(華北水利水電大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，河南 鄭州 450046)

0 引言

所謂變點(diǎn)是指觀測(cè)值在某一個(gè)位置或時(shí)間點(diǎn)發(fā)生了分布或數(shù)字特征的突然變化，這個(gè)位置或時(shí)間點(diǎn)稱為變點(diǎn)。不考慮可能的變化點(diǎn)就進(jìn)行統(tǒng)計(jì)分析很大可能會(huì)得到具有誤導(dǎo)性的結(jié)果，因此變點(diǎn)問題的研究在金融學(xué)、醫(yī)學(xué)、氣象學(xué)和計(jì)算機(jī)科學(xué)等領(lǐng)域有著廣泛的應(yīng)用價(jià)值。眾多學(xué)者對(duì)檢測(cè)變點(diǎn)的方法進(jìn)行了研究，如JAMES B J等[1]通過似然比方法檢驗(yàn)多元正態(tài)分布變點(diǎn)是否存在；李拂曉等[2]使用Pearson卡方統(tǒng)計(jì)量的二元分割方法檢驗(yàn)了多元Logistic回歸模型中存在的變點(diǎn)；陳睿軒等[3]利用非參數(shù)極大似然方法對(duì)金融時(shí)間序列中的變點(diǎn)進(jìn)行估計(jì)；CHERNOFF H等[4]應(yīng)用貝葉斯方法檢驗(yàn)正態(tài)分布中變點(diǎn)是否存在。

與經(jīng)典統(tǒng)計(jì)學(xué)相比，貝葉斯方法不僅能對(duì)參數(shù)做出有效的估計(jì)，還使得變點(diǎn)檢測(cè)變得更加簡(jiǎn)便，文獻(xiàn)[5-9]都是基于貝葉斯方法解決多變點(diǎn)問題。盡管如此，貝葉斯變點(diǎn)問題中的計(jì)算仍是一大難點(diǎn)，尤其是在高維情形下與后驗(yàn)分布相關(guān)的一些積分很難用數(shù)值方法計(jì)算。對(duì)此，馬爾可夫鏈蒙特卡洛(Markov Chain Monte Carlo，MCMC)算法作為一種重要的貝葉斯計(jì)算方法，將貝葉斯統(tǒng)計(jì)中復(fù)雜的計(jì)算簡(jiǎn)單化，在變點(diǎn)問題中得到廣泛應(yīng)用。如張晗等[10]在艾拉姆咖分布單變點(diǎn)模型中運(yùn)用MCMC算法對(duì)變點(diǎn)位置做出了有效估計(jì)；石凱等[11]采用MCMC算法為多維混合分布數(shù)據(jù)的參數(shù)估計(jì)和識(shí)別提供了一種有效的解決途徑。由于進(jìn)行變點(diǎn)識(shí)別時(shí)貝葉斯后驗(yàn)計(jì)算會(huì)產(chǎn)生不同維數(shù)的多參數(shù)子空間，1995年GREEN P J[12]提出了可逆跳躍馬爾可夫鏈蒙特卡洛(Reversible Jump Markov Chain Monte Carlo，簡(jiǎn)記為RJMCMC)算法。該采樣器實(shí)現(xiàn)了抽樣過程在不同維數(shù)的參數(shù)子空間之間跳躍，十分適用于多變點(diǎn)檢測(cè)，如ZHAO X 等[13]在層次貝葉斯框架下利用RJMCMC算法識(shí)別極端事件序列中的多個(gè)突變狀態(tài)；石永亮[14]利用RJMCMC算法對(duì)線性回歸模型的異常點(diǎn)進(jìn)行識(shí)別；劉鶴飛等[15]運(yùn)用RJMCMC算法對(duì)隱狀態(tài)個(gè)數(shù)未知的隱馬爾可夫多元正態(tài)分布進(jìn)行貝葉斯推斷；李勇等[16]利用RJMCMC確定對(duì)隱狀態(tài)個(gè)數(shù)未知的隱馬爾可夫0-1分布進(jìn)行貝葉斯推斷；范元靜等[17]利用RJMCMC算法確定泊松分布參數(shù)多變點(diǎn)模型中變點(diǎn)的個(gè)數(shù)，并在變點(diǎn)個(gè)數(shù)確定的基礎(chǔ)上，結(jié)合MCMC方法進(jìn)行參數(shù)估計(jì)。

Gamma分布經(jīng)常出現(xiàn)在各種應(yīng)用中，特別是在可靠性、生存分析和收入分配模型中。有關(guān)伽馬分布形狀參數(shù)ν的重要性也是毋庸置疑的，正如RAMANAYAKE A[18]所述：伽馬分布是降低故障率、常數(shù)，還是增加故障率，都由ν-1的值所決定，并且形狀參數(shù)ν在更新理論中也起著重要作用。關(guān)于Gamma分布參數(shù)變點(diǎn)的研究已有一些結(jié)果,如文獻(xiàn)[19-22]考慮的是典型的伽馬分布序列中的單變點(diǎn)問題；胡俊迎[23]在變點(diǎn)個(gè)數(shù)未知的情形下，利用RJMCMC算法對(duì)Gamma分布尺度參數(shù)進(jìn)行多變點(diǎn)檢測(cè)，并通過仿真模擬證明了在變點(diǎn)檢測(cè)中RJMCMC算法顯著優(yōu)于S-N方法；但是有關(guān)Gamma分布的形狀參數(shù)發(fā)生變化，以及形狀參數(shù)和尺度參數(shù)同時(shí)發(fā)生變化的研究仍然較少。因此，在Gamma分布參數(shù)多變點(diǎn)模型中找到多個(gè)變點(diǎn)的數(shù)量及位置是一個(gè)重要且具有挑戰(zhàn)性的統(tǒng)計(jì)問題。對(duì)比以往關(guān)于Gamma分布序列變點(diǎn)問題的研究，本文擬著重研究基于RJMCMC算法下的Gamma分布序列中形狀參數(shù)的多變點(diǎn)檢測(cè)，先利用算法確定變點(diǎn)個(gè)數(shù)，再通過MCMC算法估計(jì)變點(diǎn)位置和形狀參數(shù)。

1 Gamma分布形狀參數(shù)多變點(diǎn)模型

現(xiàn)有一組含有k個(gè)變點(diǎn)的觀測(cè)值序列y1,y2,…,yn,yi獨(dú)立且均服從Gamma分布。c1,c2,…,ck,將此觀測(cè)值序列分成了k+1段，其中的第j段ycj-1+1,ycj-1+2,…,ycj,服從形狀參數(shù)νj>0，尺度參數(shù)λj>0，j=1,2,…,k+1的Gamma分布?？蓪⑦@個(gè)多變點(diǎn)模型表示為

(1)

為了方便，規(guī)定c0=0，ck+1=n。式(1)中，變點(diǎn)個(gè)數(shù)為k，變點(diǎn)位置分別在c1,c2,…,ck?，F(xiàn)在先假定較為簡(jiǎn)單的一種情形：λj始終不發(fā)生變化，那么擬研究的變點(diǎn)模型(1)就是關(guān)于Gamma分布形狀參數(shù)的多變點(diǎn)檢測(cè)，需要估計(jì)的參數(shù)有2(k+1)個(gè)，分別是變點(diǎn)個(gè)數(shù)k；變點(diǎn)位置c1,c2,…,ck；形狀參數(shù)ν1,ν2,…,νk+1。

2 推斷原理及推斷過程

2.1 貝葉斯估計(jì)原理

貝葉斯統(tǒng)計(jì)學(xué)是基于總體信息、樣本信息、先驗(yàn)信息進(jìn)行的統(tǒng)計(jì)推斷。設(shè)參數(shù)θ的先驗(yàn)信息分布為π(θ)，隨機(jī)變量θ給定值時(shí)，總體的條件概率函數(shù)為p(x|θ)。樣本X和參數(shù)θ的聯(lián)合分布h(X,θ)=p(X|θ)π(θ)，利用貝葉斯公式

(2)

但是在實(shí)際問題中，上述后驗(yàn)密度分布(2)通常是比較復(fù)雜的未知形式，馬爾可夫鏈蒙特卡洛(MCMC)方法可以解決這一難題，它以目標(biāo)后驗(yàn)分布作為平穩(wěn)分布的馬爾可夫鏈生成隨機(jī)數(shù)，代替從后驗(yàn)分布中直接抽取樣本。

2.2 推斷過程

在貝葉斯框架下的變點(diǎn)分析，需要確定各參數(shù)的先驗(yàn)分布，并推導(dǎo)出其后驗(yàn)分布。根據(jù)文獻(xiàn)[12]可考慮各參數(shù)的先驗(yàn)分布。

1)變點(diǎn)個(gè)數(shù)k服從截?cái)嗟牟此煞植迹?/p>

2)從離散的均勻分布{0,1,2,…,n}上產(chǎn)生2k+1個(gè)順序統(tǒng)計(jì)量，c1,c2，…,ck作為其中的偶數(shù)階統(tǒng)計(jì)量,

其中0

3)取形狀參數(shù){ν1,ν2，…，νk+1}獨(dú)立同分布于形狀參數(shù)a和尺度參數(shù)b的Gamma分布且均與變點(diǎn)位置相互獨(dú)立，則

νj～Gamma(a,b),j=1,2,…,k+1。

由貝葉斯分層理論，可得所有未知參數(shù)的聯(lián)合先驗(yàn)分布

p(k,c1c2…ck,ν1ν2…νk+1)=p(k)p(c1c2…ck|k)p(ν1ν2…νk+1|k)，

再結(jié)合貝葉斯公式可得到樣本與參數(shù)聯(lián)合分布的密度函數(shù)

(3)

用π(ck|·)表示變點(diǎn)位置ck的滿條件分布，同理π(νj|·)表示形狀參數(shù)νj的滿條件分布。由(3)式可得各參數(shù)的滿條件分布

3 基于RJMCMC算法的參數(shù)估計(jì)

3.1 RJMCMC算法

為獲取可逆的蒙特卡洛樣本，需要設(shè)計(jì)合適的移動(dòng)類型來達(dá)到維數(shù)的平衡。當(dāng)形狀參數(shù)是一個(gè)階梯函數(shù)時(shí)，本文可基于變點(diǎn)模型設(shè)計(jì)4種移動(dòng)類型改變馬爾可夫鏈的狀態(tài){k,c1,c2,…,ck,ν1,ν2,…,νk+1},

(a)任意改變一個(gè)形狀參數(shù)值；

(b)任意改變一個(gè)變點(diǎn)的位置；

(c)在{1,2,…，n}{c1,c1,…，ck}上任意選擇新增加一個(gè)；

(d)在{c1,c1,…，ck}中任意選擇減少一個(gè)。

顯然，(c)、(d)的移動(dòng)會(huì)使得馬氏鏈的狀態(tài)空間維數(shù)隨著變點(diǎn)數(shù)量的改變而變化，因此傳統(tǒng)的蒙特卡洛馬爾可夫鏈理論無法適用。用{(a),(b),(c),(d)}表示這些移動(dòng)的集合，且用m表示移動(dòng)類型。假設(shè)當(dāng)前狀態(tài)為x，建議移動(dòng)類型為m∈{(c),(d)}，新狀態(tài)為x′，通過獲得一個(gè)獨(dú)立于x的連續(xù)隨機(jī)變量u且滿足dim(x,u)=dim(x′)，并利用可逆的確定性函數(shù)x′(x,u)設(shè)置x′后，可實(shí)現(xiàn)狀態(tài)空間的可逆跳躍。

1995年，GREEN P J探索了4種移動(dòng)類型的接受概率，并且證明了每個(gè)移動(dòng)類型都能達(dá)到馬爾可夫鏈細(xì)致平衡的既定目標(biāo)。當(dāng)移動(dòng)類型m∈{(a),(b)}的接受概率

(4)

其中π(x)為狀態(tài)x時(shí)參數(shù)后驗(yàn)分布的密度函數(shù)，qm(x′,x)為建議密度函數(shù)，為簡(jiǎn)化計(jì)算，(4)式中后驗(yàn)分布可分別用似然函數(shù)與先驗(yàn)分布的乘積項(xiàng)代替或由相應(yīng)參數(shù)的滿條件分布代替。當(dāng)移動(dòng)類型m∈{(c),(d)}的接受概率

(5)

下面探討Gamma分布形狀參數(shù)多變點(diǎn)模型在每種移動(dòng)下的接受概率。

Pallow=min{1,A1}，

這里

Pallow=min{1,A2}，

這里，當(dāng)c′j

當(dāng)c′j>cj時(shí)

對(duì)于m=(c)，不失一般性，我們假設(shè)在區(qū)間(cj-1,cj)上增加一個(gè)變點(diǎn)c*，則在區(qū)間(cj-1,c*)和(c*，cj)上會(huì)產(chǎn)生新的參數(shù)ν′j，ν′j+1，且νj在ν′j，ν′j+1之間，其關(guān)系用權(quán)重方式表示為

(c*-cj-1)logν′j+(cj-c*)logν′j+1=(cj-cj-1)logνj。

在此種移動(dòng)類型的可接受概率表示為

Pallow=min{1,(l.r.b)×(p.r.b)×(pro.r.b)×|(Jaco.b)|}，

其中l(wèi).r.b,p.r.b,pro.r.b,Jaco.b分別表示似然比、先驗(yàn)比、建議比、雅可比行列式。

經(jīng)計(jì)算，似然比可直接表示為

先驗(yàn)比

建議比

雅克比行列式

因此隨機(jī)增加一個(gè)新變點(diǎn)c*的接受概率為Pallow=min{1,A3}，這里

針對(duì)m=(d)，同樣不失一般性，假設(shè)隨機(jī)選擇被減去的變點(diǎn)為cj，則區(qū)間(cj-1,cj,cj+1)變?yōu)?cj-1,cj+1)。假設(shè)ν′j，ν′j+1為區(qū)間(cj-1,cj,cj+1)上的2個(gè)舊參數(shù)，νj為區(qū)間(cj-1,cj+1)上的新參數(shù)。同理可得，隨機(jī)減少變點(diǎn)cj的可接受概率為Pallow=min{1,A4}，這里，

3.2 參數(shù)估計(jì)

首先利用RJMCMC算法確定模型中變點(diǎn)的個(gè)數(shù)。由于變點(diǎn)位置及形狀參數(shù)的后驗(yàn)分布比較復(fù)雜，所以采用Metroplis-Hastings算法對(duì)其抽樣，再根據(jù)前面推導(dǎo)出的相應(yīng)移動(dòng)的接受概率，得到模型參數(shù)的更新值，具體步驟如下：

1)初始化，給出參數(shù)k,ck,λ的初始值，給定超參數(shù)a,b,kmax,α的值。

2)利用Metroplis-Hastings算法更新參數(shù){c1,c2,…,ck}。

3)利用Metroplis-Hastings算法更新參數(shù){ν1,ν2,…,νk+1}。

4)每次以bk或dk隨機(jī)增加或減少一個(gè)變點(diǎn)，參數(shù)發(fā)生相應(yīng)的改變，最后以概率Pallow接受變點(diǎn)個(gè)數(shù)的跳躍，實(shí)現(xiàn)可逆跳躍過程，其中當(dāng)0

重復(fù)以上迭代步驟，直到最大迭代次數(shù)。

4 實(shí)證研究

4.1 仿真模擬

隨機(jī)生成含有300個(gè)數(shù)據(jù)的Gamma分布序列，分為3段: 1～120,121～200,201～300，數(shù)據(jù)分別服從Gamma(0.5,2)，Gamma(2,2)，Gamma(5,2)。根據(jù)3段形狀參數(shù)ν不一致，可見存在2個(gè)變點(diǎn)，分別在120處和200處，3段Gamma分布參數(shù)分別為(0.5,2),(2,2),(5,2)。300個(gè)隨機(jī)數(shù)據(jù)模擬圖如圖1所示。設(shè)定參數(shù)的初始值k=3，c1=50，c2=150，c3=180，ν1,ν2，…，ν4取初始分段上的均值，超參數(shù)kmax=10，α=5a=25/4，b=5/4。迭代20 000次算法，去掉前10 000次，用后10 000次的結(jié)果計(jì)算變點(diǎn)個(gè)數(shù)的后驗(yàn)概率，得出的變點(diǎn)個(gè)數(shù)估計(jì)如圖2所示。

圖1 兩變點(diǎn)的Gamma分布數(shù)據(jù)模擬圖

圖2 變點(diǎn)個(gè)數(shù)估計(jì)直方圖

由圖2可知，變點(diǎn)個(gè)數(shù)為2的后驗(yàn)概率最大，因此確定300個(gè)Gamma分布序列的變點(diǎn)個(gè)數(shù)為2。進(jìn)一步在變點(diǎn)個(gè)數(shù)的基礎(chǔ)上利用MCMC算法估計(jì)變點(diǎn)位置參數(shù)和Gamma分布形狀參數(shù)。通過R軟件實(shí)現(xiàn)模擬，在模擬的過程中進(jìn)行50 000次迭代抽樣，為保證參數(shù)的收斂性，舍棄前30 000次抽樣，根據(jù)后20 000次結(jié)果進(jìn)行統(tǒng)計(jì)分析，20 000次變點(diǎn)位置參數(shù)和形狀參數(shù)迭代過程如圖3和圖4所示。

圖3 變點(diǎn)位置c1的迭代抽樣過程

圖4 變點(diǎn)位置c2的迭代抽樣過程

兩變點(diǎn)位置的后驗(yàn)密度分布如圖5所示，兩個(gè)變點(diǎn)下的形狀參數(shù)的后驗(yàn)分布如圖6所示。由圖5可知，變點(diǎn)位置的后驗(yàn)密度分布有2個(gè)峰，分別在120，200附近;由圖6可知，形狀參數(shù)的后驗(yàn)密度分布有3個(gè)峰，分別在0.5，2，5附近。由眾數(shù)后驗(yàn)估計(jì)可得變點(diǎn)位置c1,c2為120、200，3個(gè)形狀參數(shù)ν1,ν2,ν3近似為0.5、2、5。這與真實(shí)的變點(diǎn)位置及Gamma分布形狀參數(shù)相符，說明RJMCMC算法對(duì)Gamma分布形狀參數(shù)多變點(diǎn)檢測(cè)的有效性。

圖5 兩變點(diǎn)位置的后驗(yàn)密度分布

圖6 兩變點(diǎn)形狀參數(shù)的后驗(yàn)密度分布

4.2 英國(guó)煤礦災(zāi)害

MAGUIRE B A等[24-25]對(duì)英國(guó)在1875年12月6日至1951年5月29日期間發(fā)生的造成10人以上死亡的煤礦爆炸的時(shí)間間隔進(jìn)行了研究(如圖7所示)，得出結(jié)論：這組數(shù)據(jù)可能存在一個(gè)變點(diǎn)在1890年。隨后RAMANAYAKE A[18]又驗(yàn)證了這109個(gè)數(shù)據(jù)服從Gamma分布，進(jìn)一步證實(shí)了在第45個(gè)位置(即對(duì)應(yīng)1890年)發(fā)生了突變并估計(jì)出突變前的形狀參數(shù)為0.629，突變后的形狀參數(shù)為1.036。

圖7 英國(guó)煤礦爆炸數(shù)據(jù)折線圖

現(xiàn)在根據(jù)RJMCMC算法對(duì)這組數(shù)據(jù)的變點(diǎn)個(gè)數(shù)進(jìn)行檢測(cè)，變點(diǎn)個(gè)數(shù)的估計(jì)如圖8所示，109個(gè)數(shù)據(jù)中存在1個(gè)變點(diǎn)的后驗(yàn)概率最大，表明1875年至1951年間發(fā)生了一次明顯的變化。接下來在變點(diǎn)個(gè)數(shù)為1的基礎(chǔ)上估計(jì)變點(diǎn)位置c1及變點(diǎn)前后的形狀參數(shù)ν1、ν2。

圖8 變點(diǎn)個(gè)數(shù)估計(jì)直方圖

變點(diǎn)位置c1的迭代抽樣過程如圖9所示，變點(diǎn)位置的后驗(yàn)密度分布如圖10所示，變點(diǎn)前后的形狀參數(shù)ν1、ν2的后驗(yàn)密度分布如圖11所示。由圖10和圖11，利用最大后驗(yàn)估計(jì)法進(jìn)行分析，可知變點(diǎn)的位置在44附近(即在1890年附近)，變點(diǎn)前后的形狀參數(shù)ν1在0.63附近且ν2在1附近。這與 RAMANAYAKE A[18]的研究結(jié)果基本相符，進(jìn)一步表明了該方法對(duì)伽馬分布參數(shù)變點(diǎn)檢測(cè)的有效性。

圖9 變點(diǎn)位置的迭代抽樣圖

圖10 變點(diǎn)位置c1的后驗(yàn)分布圖

圖11 形狀參數(shù)ν1、ν2的后驗(yàn)分布圖

5 總結(jié)

本文基于貝葉斯方法和RJMCMC算法對(duì)Gamma分布形狀參數(shù)多變點(diǎn)模型進(jìn)行變點(diǎn)檢測(cè)。首先給出了Gamma分布形狀參數(shù)多變點(diǎn)模型，確定各參數(shù)的先驗(yàn)分布，由貝葉斯分層理論將樣本與參數(shù)的聯(lián)合分布以乘積項(xiàng)的形式表示。然后針對(duì)變點(diǎn)模型設(shè)計(jì)4種移動(dòng)，并分別得到其接受概率的表達(dá)式，建立RJMCMC算法。最后利用該算法和MCMC算法得到變點(diǎn)個(gè)數(shù)的確定以及變點(diǎn)位置和形狀參數(shù)的后驗(yàn)估計(jì)。仿真模擬和實(shí)例的結(jié)果都說明了RJMCMC算法對(duì)Gamma分布形狀參數(shù)多變點(diǎn)模型變點(diǎn)個(gè)數(shù)確定的有效性。

目前，研究還存在一些不足之處：①?zèng)]有進(jìn)一步研究RJMCMC算法的收斂性，可能會(huì)出現(xiàn)由于計(jì)算量較大，算法耗時(shí)較長(zhǎng)的情形；②本文考慮的Gamma分布多變點(diǎn)模型有尺度參數(shù)不發(fā)生變化的假定，算法僅實(shí)現(xiàn)了變點(diǎn)位置和形狀參數(shù)的估計(jì)，而實(shí)際數(shù)據(jù)往往與假定條件存在著些許偏差，所以如何利用RJMCMC算法更好地結(jié)合伽馬分布的兩個(gè)參數(shù)進(jìn)行的多變點(diǎn)檢測(cè)有待進(jìn)一步研究。