李 勇,劉鶴飛

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

0 引言

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

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

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

1 模型描述

觀測變量yit來自0-1分布，其中i=1，2，...，N是觀測對象，t=1，2，...，T是觀測時點。Zit是系統(tǒng)觀測對象在第t個觀測時點的隱狀態(tài)，其取值范圍為{1，2，...，K}，稱為K個隱狀態(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是從前一個時點的隱狀態(tài)u向后一個時點的隱狀態(tài)s的轉(zhuǎn)移概率。全部可能的隱狀態(tài)轉(zhuǎn)移概率構(gòu)成一個矩陣，稱為隱狀態(tài)轉(zhuǎn)移概率矩陣，記為：

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

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

其中θj是第j個隱狀態(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)個數(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）將一個隱狀態(tài)分裂為兩個，或者將兩個隱狀態(tài)合并成一個。

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

2.2 先驗分布

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

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

θj～U（0，1）

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

2.3 全條件后驗分布

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

隱狀態(tài)Z的全條件后驗分布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），又因為A的每一行元素的先驗分布都是對稱的狄尼克萊分布D（α，α，...，α），所以每一行的全條件后驗分布為：

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

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

接下來推導(dǎo)0-1分布參數(shù)θ的全條件后驗分布。當(dāng)隱狀態(tài)Z確定時，隱狀態(tài)轉(zhuǎn)移概率矩陣A也隨之確定，此時第j個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的先驗分布是U（0，1），則其概率密度函數(shù)為：

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

Yj的邊際密度為：

利用貝葉斯公式可得θj的后驗分布：

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

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

2.4 可逆跳躍過程

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

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

M=似然比×先驗比×跳躍比×建議比×雅克比

對于隱馬爾可夫0-1分布來說，假設(shè)當(dāng)前的隱狀態(tài)個數(shù)為k，則有一個k階的隱狀態(tài)轉(zhuǎn)移概率矩陣，有k個在0～1之間的參數(shù)θj。根據(jù)可逆跳躍MCMC算法的要求，增加一個隱狀態(tài)，需要將k階的隱狀態(tài)轉(zhuǎn)移概率矩陣增加到k+1階，并且增加一個0～1之間的參數(shù)θ。將選中的待分解的隱狀態(tài)記為j*，將j*分解之后的兩個隱狀態(tài)分別記為j1和j2。下面分別介紹生成隨機書分解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）個自由元素，k+1階轉(zhuǎn)移概率矩陣有（k+1）k個自由元素。下面介紹如何生成k（k+1）-k（k-1）=2 個隨機數(shù)，對k階轉(zhuǎn)移概率矩陣進行分解。

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

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

再從Beta（r，s）中生成隨機數(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)介紹過了，最后來介紹隨機數(shù)t1的生成方法。

令：

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

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

最后來計算這種分解跳躍的接受概率。因為接受概率為min（1，M），且M=似然比×先驗比×跳躍比×建議比×雅克比。每部分的表達式分別為：

其中，z1t是分解之前第i個觀測對象在第t個觀測點的隱狀態(tài)，是分解之后第i個觀測對象在第t個觀測點的隱狀態(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 實證模擬

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

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

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

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

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

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

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

4 結(jié)束語

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

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