李 勇,劉鶴飛

（曲靖師范學(xué)院a.信息工程學(xué)院；b.數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，云南曲靖 655011）

0 引言

隱馬爾可夫模型是一種基于隨機(jī)過程的統(tǒng)計(jì)學(xué)模型。隱馬爾可夫模型由兩個(gè)隨機(jī)過程構(gòu)成，一個(gè)是狀態(tài)轉(zhuǎn)移序列、一條單純的馬爾可夫鏈；另一個(gè)是與狀態(tài)對應(yīng)的觀測序列。在實(shí)際問題的研究中，只能看到觀測序列的集合，而不能看到狀態(tài)序列集合。也就是說，模型的狀態(tài)隱藏在觀測序列中，因此稱之為“隱”馬爾可夫模型[1]。隱馬爾可夫模型的研究內(nèi)容就是通過可觀測的變量序列集合去推斷不可觀測的狀態(tài)序列轉(zhuǎn)移特征及每個(gè)狀態(tài)下的分布信息。

近年來，馬爾可夫模型在許多領(lǐng)域都有了極大的發(fā)展和廣泛的應(yīng)用。例如，在圖像處理領(lǐng)域，沈杰等[2]利用隱馬爾可夫模型改進(jìn)了人臉識別技術(shù)；在語音識別領(lǐng)域，魏明哲[3]通過隱馬爾可夫模型的分類識別提高了語音識別的效果；在生物醫(yī)學(xué)領(lǐng)域，劉青等[4]通過隱馬爾可夫模型實(shí)現(xiàn)基因識別系統(tǒng)的設(shè)計(jì)。在經(jīng)濟(jì)金融領(lǐng)域，孫鵬飛等[5]利用隱馬爾可夫模型來擬合美元LIBOR相對美聯(lián)儲基準(zhǔn)利率的變動(dòng)情況。

雖然近年來關(guān)于隱馬爾科夫模型的研究成果比較豐富，但是由于隱馬爾可夫模型的復(fù)雜性，大多數(shù)研究都是關(guān)于隱馬爾科夫模型的應(yīng)用性研究，對于隱馬爾科夫模型的純理論性研究非常少。本文在已有的關(guān)于隱馬爾科夫模型理論研究成果的基礎(chǔ)上，研究隱狀態(tài)個(gè)數(shù)未知條件下的隱馬爾可夫0-1分布模型。首先利用可逆跳躍MCMC[6]方法對未知的隱狀態(tài)個(gè)數(shù)進(jìn)行貝葉斯模型選擇；確定隱狀態(tài)個(gè)數(shù)之后再用MCMC算法估計(jì)隱馬爾可夫0-1分布的參數(shù)；最后生成觀測數(shù)據(jù)集，實(shí)證模擬檢驗(yàn)該方法的推斷效果。

1 模型描述

觀測變量yit來自0-1分布，其中i=1，2，...，N是觀測對象，t=1，2，...，T是觀測時(shí)點(diǎn)。Zit是系統(tǒng)觀測對象在第t個(gè)觀測時(shí)點(diǎn)的隱狀態(tài)，其取值范圍為{1，2，...，K}，稱為K個(gè)隱狀態(tài)的隱馬爾可夫模型。

假設(shè)模型不可觀測的隱式鏈滿足以下馬爾可夫鏈條件：P（Zit=s│Zi1，Zi2，...，Zi（t-1）=u）=P（Zit=s│Zi（t-1）=u）=aus

其中：

u=1，2，...，K;s=1，2，...，K;i=1，2，...，N;t=2，3，...，T。

則aus是從前一個(gè)時(shí)點(diǎn)的隱狀態(tài)u向后一個(gè)時(shí)點(diǎn)的隱狀態(tài)s的轉(zhuǎn)移概率。全部可能的隱狀態(tài)轉(zhuǎn)移概率構(gòu)成一個(gè)矩陣，稱為隱狀態(tài)轉(zhuǎn)移概率矩陣，記為：

隱狀態(tài)轉(zhuǎn)移概率矩陣的每一行之和為1。

在給定隱狀態(tài)的條件下，隱馬爾可夫0-1分布的觀測變量定義為：

其中θj是第j個(gè)隱狀態(tài)下0-1分布的參數(shù)，yit是0-1分布的觀測度量，其取值僅為0或1。

2 貝葉斯推斷

2.1 推斷原理

記θ為不同隱狀態(tài)下0-1分布的參數(shù)集合，即θ=（θ1，θ1，...，θk）。 則 該 模 型 的 貝 葉 斯 推 斷 問 題 為P（K，Z，A，θ|Y），就是要根據(jù)觀測度量的信息推斷出模型的因狀態(tài)個(gè)數(shù)K，隱狀態(tài)集合信息Z，轉(zhuǎn)移概率矩陣A，以及不同隱狀態(tài)下0-1分布的參數(shù)θ。利用可逆跳躍馬爾可夫鏈蒙特卡羅算法推斷的具體執(zhí)行步驟為：

（1）更新隱狀態(tài)Z；

（2）更新隱狀態(tài)轉(zhuǎn)移概率矩陣A；

（3）更新0-1分布的參數(shù)θ；

（4）將一個(gè)隱狀態(tài)分裂為兩個(gè)，或者將兩個(gè)隱狀態(tài)合并成一個(gè)。

步驟（1）至步驟（3）中每一步更新參數(shù)都是普通的馬爾可夫鏈蒙特卡羅算法，主要利用Gibbs抽樣[7]和MH算法[8]從相關(guān)參數(shù)的全條件后分布抽樣即可。步驟（4）是可逆跳躍步，在這一步中，允許以pk和qk=1-pk的概率增加或者減少一個(gè)隱狀態(tài)。一般情況下，使增加一個(gè)隱狀態(tài)和減少一個(gè)隱狀態(tài)的概率相等，即pk=qk=0.5，但是q2=pmax=0除外。其中，q2=0表示隱狀態(tài)個(gè)數(shù)為2時(shí)，再減少隱狀態(tài)的概率為0；pmax=0表示隱狀態(tài)達(dá)到設(shè)定的最大值時(shí)，再增加隱狀態(tài)的概率為0；當(dāng)隱狀態(tài)為其他值時(shí)，增加一個(gè)隱狀態(tài)和減少一個(gè)隱狀態(tài)的概率是相等的，都是0.5。

2.2 先驗(yàn)分布

貝葉斯推斷中，所有參數(shù)都看成是隨機(jī)變量，所以，在進(jìn)行貝葉斯推斷時(shí)，要先為每一個(gè)參數(shù)選擇一個(gè)合適的先驗(yàn)分布。本文借鑒已有的研究經(jīng)驗(yàn)[9]，為各個(gè)參數(shù)選取的先驗(yàn)分布如下：

Ak～D（α，α，...，α）

θj～U（0，1）

其中，Ak表示隱狀態(tài)轉(zhuǎn)移概率矩陣的第k行，D表示狄尼克萊分布，α是狄尼克萊分布的超參數(shù)。θj表示第j個(gè)隱狀態(tài)下0-1分布的參數(shù)，為其選擇的先驗(yàn)分布是（0，1）上的均勻分布。

2.3 全條件后驗(yàn)分布

貝葉斯后驗(yàn)推斷主要是利用樣本信息和參數(shù)的先驗(yàn)信息對模型進(jìn)行推斷。對于本文來說，即P（K，Z，A，θ|Y）。由于這個(gè)后驗(yàn)分布的復(fù)雜性，本文利用馬爾可夫鏈蒙特卡羅方法來對后驗(yàn)分布進(jìn)行模擬，這就需要推導(dǎo)相關(guān)參數(shù)的全條件后驗(yàn)分布。

隱狀態(tài)Z的全條件后驗(yàn)分布P（Z|A，θ，Y）為：

也即：

其中，fk（Yt|θk）是觀測變量在隱狀態(tài)k下的似然函數(shù)。

由于隱狀態(tài)轉(zhuǎn)移概率只與隱狀態(tài)的取值情況有關(guān)，所以P（A|Y，Z，θ）=P（A|Z），又因?yàn)锳的每一行元素的先驗(yàn)分布都是對稱的狄尼克萊分布D（α，α，...，α），所以每一行的全條件后驗(yàn)分布為：

P（Ak|Z）～D（α+nk1，...，α+nkK）

其中nki是前一個(gè)隱狀態(tài)為k，后一個(gè)隱狀態(tài)為i的樣本總個(gè)數(shù)。

接下來推導(dǎo)0-1分布參數(shù)θ的全條件后驗(yàn)分布。當(dāng)隱狀態(tài)Z確定時(shí)，隱狀態(tài)轉(zhuǎn)移概率矩陣A也隨之確定，此時(shí)第j個(gè)0-1分布的參數(shù)θj只與觀測變量Yj有關(guān)，即P（θj|Y，A，Z）=P（θj|Yj），其中，Yj是隱狀態(tài)為j的觀測變量集合。

由于Yj1，Yj2，...，Yjn～b（1，θj），則Yj～b（nj，θj），也即：

由于θj的先驗(yàn)分布是U（0，1），則其概率密度函數(shù)為：

則樣本Yj與參數(shù)θj的聯(lián)合概率密度函數(shù)為：。

Yj的邊際密度為：

利用貝葉斯公式可得θj的后驗(yàn)分布：

此式是參數(shù)為χ+1，nj-χ+1的貝塔分布，即θj的全條件后驗(yàn)分布為：

θj～Beta（χ+1，nj-χ+1）

2.4 可逆跳躍過程

1995年，Green[7]提出了可逆跳躍馬爾可夫鏈蒙特卡羅算法，該方法通過在可變維參數(shù)空間中跳躍式抽樣，實(shí)現(xiàn)貝葉斯模型選擇的目的。1997年，Richardson[10]和Green利用可逆跳躍MCMC算法對正態(tài)混合模型中未知的混合個(gè)數(shù)進(jìn)行選擇。2000年，Robert[11]利用可逆跳躍MCMC算法對隱馬爾可夫正態(tài)分布中未知的隱狀態(tài)個(gè)數(shù)進(jìn)行選擇。本文借鑒以上研究方法，利用可逆跳躍MCMC算法對隱馬爾可夫0-1分布中未知的隱狀態(tài)個(gè)數(shù)進(jìn)行模型選擇。

可逆跳躍MCMC算法與普通的MCMC算法的主要區(qū)別是，在可逆跳躍步，允許未知的隱狀態(tài)個(gè)數(shù)隨機(jī)地增加或者減少一個(gè)，模型的參數(shù)發(fā)生相應(yīng)的變化，最后以一定的概率接受或者拒接這種跳躍。以增加一個(gè)隱狀態(tài)為例，可逆跳躍步的具體執(zhí)行方法為：（1）從當(dāng)前的k個(gè)隱狀態(tài)中等概率隨機(jī)的選擇一個(gè)；（2）從預(yù)先給定的分布中生成一定數(shù)量的隨機(jī)數(shù)；（3）利用生成的隨機(jī)數(shù)，根據(jù)相應(yīng)的分解方式，將選中的隱狀態(tài)分解成兩個(gè)；（4）以概率min（1，M）接受這種跳躍，其中：

M=似然比×先驗(yàn)比×跳躍比×建議比×雅克比

對于隱馬爾可夫0-1分布來說，假設(shè)當(dāng)前的隱狀態(tài)個(gè)數(shù)為k，則有一個(gè)k階的隱狀態(tài)轉(zhuǎn)移概率矩陣，有k個(gè)在0～1之間的參數(shù)θj。根據(jù)可逆跳躍MCMC算法的要求，增加一個(gè)隱狀態(tài)，需要將k階的隱狀態(tài)轉(zhuǎn)移概率矩陣增加到k+1階，并且增加一個(gè)0～1之間的參數(shù)θ。將選中的待分解的隱狀態(tài)記為j*，將j*分解之后的兩個(gè)隱狀態(tài)分別記為j1和j2。下面分別介紹生成隨機(jī)書分解A和θ的具體方法。

Robert[12]根據(jù)矩陣的特征值，提出了一種能保留轉(zhuǎn)移狀態(tài)特征的轉(zhuǎn)移概率矩陣的分解方法。以k階隱狀態(tài)轉(zhuǎn)移概率矩陣分解為k+1階為例，由于轉(zhuǎn)移概率矩陣的每一行元素之和為1，因此k階轉(zhuǎn)移概率矩陣有k（k-1）個(gè)自由元素，k+1階轉(zhuǎn)移概率矩陣有（k+1）k個(gè)自由元素。下面介紹如何生成k（k+1）-k（k-1）=2 個(gè)隨機(jī)數(shù)，對k階轉(zhuǎn)移概率矩陣進(jìn)行分解。

首先生成t0～Beta（2，2），然后根據(jù)以下公式計(jì)算r和s:

其中，c2是Beta（r，s）的方差，根據(jù)經(jīng)驗(yàn)，選擇c2=0.3。

再從Beta（r，s）中生成隨機(jī)數(shù)，使得：

uj～Beta（r，s），j=1，2，...，k且j≠j1，j2

vj～Beta（r，s），j=1，2，...，k且j≠j1，j2

記分裂前的轉(zhuǎn)移概率矩陣A=aij，i，j=1，2，...，k；分裂后的轉(zhuǎn)移概率矩陣為B=bij，i，j=1，2，...，k+1,則有：

上式中，t0、uj、vj的生成方法已經(jīng)介紹過了，最后來介紹隨機(jī)數(shù)t1的生成方法。

令：

如果＜，則沒有符合條件的t1，隱狀態(tài)轉(zhuǎn)移概率矩陣分解不成功，分解跳躍過程立即停止。如果＜，則在區(qū)間[，]上，隨機(jī)地選擇一個(gè)值作為t1。

以上就是通過生成隨機(jī)數(shù)，將k階隱狀態(tài)轉(zhuǎn)移概率矩陣分解成k+1階的詳細(xì)過程，這種分解方式不僅能保證每行元素之和為1，而且是可逆的，更重要的是能最大限度保留原來的轉(zhuǎn)移狀態(tài)特征。生成隨機(jī)數(shù)，將0-1分布的參數(shù)θj*分解成θj1和θj2的方法比較簡單，只需令ε～U（0，1），θj2=ε*θj*，θj2=（1-ε）*θj*即可。

最后來計(jì)算這種分解跳躍的接受概率。因?yàn)榻邮芨怕蕿閙in（1，M），且M=似然比×先驗(yàn)比×跳躍比×建議比×雅克比。每部分的表達(dá)式分別為：

其中，z1t是分解之前第i個(gè)觀測對象在第t個(gè)觀測點(diǎn)的隱狀態(tài)，是分解之后第i個(gè)觀測對象在第t個(gè)觀測點(diǎn)的隱狀態(tài)。

建議比 =[B1，2（t0）ΠjBr，s（uj）ΠjBr，s（vj）]-1

其中，B2，2表示Beta（2,2）分布的密度函數(shù)；Br，s表示Beta（r,s）分布的密度函數(shù)。

3 實(shí)證模擬

根據(jù)隱馬爾可夫0-1分布的數(shù)學(xué)定義，生成一個(gè)觀測變量數(shù)據(jù)集。首先利用可逆跳躍MCMC算法對未知的隱狀態(tài)個(gè)數(shù)進(jìn)行模型選擇；然后在隱狀態(tài)個(gè)數(shù)已知的情況下利用MCMC算法對模型的參數(shù)進(jìn)行貝葉斯估計(jì)；最后將估計(jì)結(jié)果與模型參數(shù)的真實(shí)值進(jìn)行比較，觀察所給方法的推斷效果。

取真實(shí)的隱狀態(tài)個(gè)數(shù)k=2，觀測對象個(gè)數(shù)N=500，觀測時(shí)間T=10，轉(zhuǎn)移概率矩陣兩個(gè)0-1分布的參數(shù)分別為θ1=0.4，θ2=0.7。

取狄尼克萊分布中的超參數(shù)α=1，在使用可逆跳躍MCMC算法時(shí)，取迭代總次數(shù)為5萬次，去掉前面的3萬次，用后面2萬次的結(jié)果來計(jì)算未知隱狀態(tài)個(gè)數(shù)的后驗(yàn)概率。圖1是5萬次迭代過程隱狀態(tài)個(gè)數(shù)的迭代路徑圖，表1（見下頁）是隱狀態(tài)個(gè)數(shù)的后驗(yàn)概率。

圖1 隱狀態(tài)個(gè)數(shù)迭代路徑圖

表1 隱狀態(tài)后驗(yàn)概率

根據(jù)隱狀態(tài)個(gè)數(shù)的后驗(yàn)概率結(jié)果，本文選擇K=2作為未知的隱狀態(tài)個(gè)數(shù)。然后再用MCMC算法取估計(jì)隱狀態(tài)的轉(zhuǎn)移概率和每個(gè)隱狀態(tài)下的0-1分布參數(shù)值。超參數(shù)取值不變，MCMC算法的迭代總次數(shù)為8千次，去掉前面4千次，用后面4千次的后驗(yàn)均值作為未知參數(shù)的估計(jì)值，估計(jì)結(jié)果見表2。

表2 參數(shù)估計(jì)結(jié)果

4 結(jié)束語

本文研究了隱狀態(tài)個(gè)數(shù)未知情況下的隱馬爾可夫0-1分布的貝葉斯推斷。首先為模型中的所有參數(shù)選擇了合適的先驗(yàn)分布，并推導(dǎo)出了各參數(shù)的全條件后驗(yàn)分布，然后詳細(xì)介紹了利用可逆跳躍MCMC進(jìn)行模型選擇時(shí)，可逆跳躍步的具體執(zhí)行過程，包括隨機(jī)數(shù)的生成方式、參數(shù)的分解方式、跳躍的接受概率計(jì)算公式等。隱狀態(tài)個(gè)數(shù)確定之后，再用普通的MCMC算法估計(jì)模型中未知參數(shù)的值。最后，將模型的貝葉斯推斷結(jié)果與設(shè)定的真實(shí)模型進(jìn)行比較，發(fā)現(xiàn)利用可逆跳躍MCMC算法選出的隱狀態(tài)個(gè)數(shù)與設(shè)定的真實(shí)隱狀態(tài)個(gè)數(shù)一致，并且利用MCMC算法計(jì)算出的參數(shù)的后驗(yàn)均值也與模型設(shè)定的參數(shù)的真實(shí)值非常接近，這說明，該模型的貝葉斯推斷效果良好。

本文的模型中考慮的是最簡單最常用的0-1分布，0-1分布比較簡單，分布中只含有一個(gè)未知參數(shù)，但其觀測變量的取值只有0和1兩個(gè)值，也導(dǎo)致觀測變量的信息比較單調(diào)，所以該模型的推斷也有一定的困難。今后的研究，可以在此基礎(chǔ)上考慮更復(fù)雜的分布。