郭衛(wèi)娟

(湖北第二師范學(xué)院a.數(shù)學(xué)與經(jīng)濟學(xué)院；b.大數(shù)據(jù)建模與智能計算研究所，武漢 430205)

1 變點問題簡介

變點序列數(shù)據(jù)是數(shù)理統(tǒng)計中經(jīng)常遇見的一個序列，在該序列中，各個子部分的總體的分布并不是一樣的，對于這類問題，通常的處理方式是先識別該序列中的變點的位置，然后就可以利用相鄰的兩個變點之間的分布是相同的，進而來估計該部分的分布。

其一般模型如下：

(1)

2 狀態(tài)轉(zhuǎn)移概率

為此建立如下模型：在時刻，定義狀態(tài)(i=0,1,2,…,t-1)表示離t最近的前向變點位置在t-i位置上，記其概率為p(Ct=i|xt-1,xt-i+1,…xt-1)，意思即xt-i,xt-i+1,…xi-1這個觀測值是獨立同分布，例如(Ct=0|xt-i,xt-i+1,…xt)表示xi-1是變點，i=t-1表示該序列無變點。這與傳統(tǒng)的馬爾科夫鏈相比，就是將隱馬爾科夫鏈中有限個狀態(tài)改成成了與當(dāng)前時刻t相關(guān)的一個變量。這樣將會導(dǎo)致轉(zhuǎn)移概率矩陣維數(shù)無限增大，因此為了最大程度上簡化狀態(tài)轉(zhuǎn)移概率矩陣，為此筆者再假設(shè)模型(1)滿足如下特征：

p(Ct=i|xt-i,xt-i+1,…，xt-1)=p(Ct-k=i|xt-k-i,xt-k-i+1,…xt-k+1)，

也就是連續(xù)的i個觀測值是同一分布(大部分參考文獻稱該值為鏈長，用字母g表示)與該觀測值的起點位置無關(guān)，這樣，整個狀態(tài)概率概率就簡單的由鏈長的概率分布確定了?？紤]到本文是從當(dāng)前時刻開始，逐步向前查找最近的變點 ，若令p表示每個觀察值可能是變點的概率，即

p(xi是變點)=p

=p(1-p)i-1i=0,1,2,…，t-1。即此時鏈長g服從幾何分布Ge(p) 。

實際上為鏈長g可以為取值于i=0,1,2,…，t-1的任意離散型分布，同樣可以計算該分布的生存函數(shù)。利用生存函數(shù)可以計算出各個狀態(tài)之間的轉(zhuǎn)移概率。

則有各個狀態(tài)之間的轉(zhuǎn)移概率計算如下：

(1)若位置t-1是為變點，則此時離t最近的前向變點就是t-1，此時j=0,

(2)若位置t-1是為變點，則此時離t最近的前向變點就是t-1，此時j=i+1,

特別的，若鏈長g服從幾何分布Ge(p)，則利用(2)式可知其對應(yīng)的狀態(tài)轉(zhuǎn)移概率為：

TP=(t=j│t-1=i)

3 模擬和計算

因為該模型的重點是識別變點的位置，也就是識別當(dāng)前時刻該序列所處的狀態(tài)，因此按照隱馬爾科夫鏈模型，主要是學(xué)習(xí)該模型的參數(shù)，然后采用最大后驗概率進行模式(隱含狀態(tài))識別問題。實際上，本文的主要工作就是從最后一個觀察值，采用前向傳導(dǎo)算法找出該序列中所有的變點位置。也就是主要是求給定觀測值下鏈長的概率分布。為此采用貝葉斯方法。方法如下：

(1) 初始化令p(C1=0│x1) =0

(2)遞推公式：

上式p(Ct=i|x1,x2,…,xt)和p(Ct-1=i|x1,x2,…,xt-1)形式一致，因此可以建立二者之間的遞推關(guān)系，若記b(t,i)=p(Ct=i|x1,x2,…,xt)，則有遞推公式

由此可以計算出全部的b(t,i),i≤t≤n。

(3)隨機模擬方法：令T0=n，k=0,從b(T0,i)抽樣，得到Tk，然后令k=k+1，若Tk>0，則繼續(xù)從b(Tk-1,i)中抽樣，這樣就得到一序列變點位置，Tk-1,Tk-2,…T1。為提高精度，本文重復(fù)抽樣1000次，最后用均值估計Tk-1,Tk-2,…,T1。

顯然模擬數(shù)據(jù)一共是400個，有4個變點，位置分別位于第100,200,300,400處。因此假設(shè)p=0.04.

圖1 有4個變點的實際圖像

實際上，該方法主要問題是求出全部變點位置，而對參數(shù)估計值并未做出更多改進，因而，參數(shù)估計部分由于不同，必然會帶來一定的估計偏差，這是筆者下一步努力的方向?？偟膩碚f，該方法不失為多個變點位置估計的一種好方法。