程慧慧, 許淑月

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

0 引言

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

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

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

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

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

(1)

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

2 推斷原理及推斷過程

2.1 貝葉斯估計原理

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

(2)

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

2.2 推斷過程

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

1)變點個數(shù)k服從截斷的泊松分布，

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

其中0

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

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

由貝葉斯分層理論，可得所有未知參數(shù)的聯(liá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|·)表示變點位置ck的滿條件分布，同理π(νj|·)表示形狀參數(shù)νj的滿條件分布。由(3)式可得各參數(shù)的滿條件分布

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

3.1 RJMCMC算法

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

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

(b)任意改變一個變點的位置；

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

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

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

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

(4)

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

(5)

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

Pallow=min{1,A1}，

這里

Pallow=min{1,A2}，

這里，當c′j

當c′j>cj時

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

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

在此種移動類型的可接受概率表示為

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分別表示似然比、先驗比、建議比、雅可比行列式。

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

先驗比

建議比

雅克比行列式

因此隨機增加一個新變點c*的接受概率為Pallow=min{1,A3}，這里

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

3.2 參數(shù)估計

首先利用RJMCMC算法確定模型中變點的個數(shù)。由于變點位置及形狀參數(shù)的后驗分布比較復雜，所以采用Metroplis-Hastings算法對其抽樣，再根據(jù)前面推導出的相應移動的接受概率，得到模型參數(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隨機增加或減少一個變點，參數(shù)發(fā)生相應的改變，最后以概率Pallow接受變點個數(shù)的跳躍，實現(xiàn)可逆跳躍過程，其中當0

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

4 實證研究

4.1 仿真模擬

隨機生成含有300個數(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個變點，分別在120處和200處，3段Gamma分布參數(shù)分別為(0.5,2),(2,2),(5,2)。300個隨機數(shù)據(jù)模擬圖如圖1所示。設定參數(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é)果計算變點個數(shù)的后驗概率，得出的變點個數(shù)估計如圖2所示。

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

圖2 變點個數(shù)估計直方圖

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

圖3 變點位置c1的迭代抽樣過程

圖4 變點位置c2的迭代抽樣過程

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

圖5 兩變點位置的后驗密度分布

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

4.2 英國煤礦災害

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

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

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

圖8 變點個數(shù)估計直方圖

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

圖9 變點位置的迭代抽樣圖

圖10 變點位置c1的后驗分布圖

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

5 總結(jié)

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

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