非齊次隱馬爾可夫模型及其參數估計

徐妍

摘要：隱馬爾可夫模型是現今被廣泛使用的統計模型之一。本文在現有的對隱馬爾可夫模型研究的基礎上，通過在轉移模型的馬爾可夫鏈中增加協變量，構建了非齊次隱馬爾可夫模型。在對模型進行參數估計時，首先用k-means聚類分析確定了觀測數據的隱狀態，然后用極大似然估計對模型的轉移參數進行估計。在數值模擬時，以非齊次隱馬爾可夫多元正態分布和非齊次隱馬爾可夫多元線性回歸為例，利用文章所介紹的方法對模型的參數進行估計，驗證了估計的可靠性。

關鍵詞：非齊次隱馬爾可夫模型 k-means聚類 極大似然估計

一、引言

隱馬爾可夫模型最早是由Leonard E.Baum等在1966年和1970年的兩篇論文中提出的[1，2]。隱馬爾可夫模型由兩部分組成，一部分是轉移模型，用來描述模型隱狀態之間的轉移關系，在實際應用中隱狀態是不能被直接觀測到的;另一部分是發射模型，即和隱狀態一一對應的觀測變量，每個發射模型來自一個特定的分布。20世紀70年代，隱馬爾可夫模型首先被應用在語音識別領域[3]，之后，又成為分析生物遺傳信息的有用工具[4]。隨著隱馬爾可夫模型的不斷發展，對隱馬爾可夫模型的研究也越來越多。Dempster等提出了EM算法，用來對包含隱狀態的模型進行極大似然估計[5]。Jelinek等將Baum-Welch算法與EM算法的關系進行了完整的描述，說明了Baum-Welch算法是EM算法在隱馬爾可夫模型中的具體應用[6]。國內也有很多關于隱馬爾可夫模型的研究，夏葉茂等研究了隱馬爾可夫因子分析模型的半參數貝葉斯分析 [7]。劉鶴飛等研究了隱狀態個數未知的隱馬爾可夫多元正態分布的貝葉斯推斷，利用可逆跳躍MCMC算法對隱狀態個數進行模型選擇，再對模型參數進行估計[8]。王坤等將隱馬爾可夫模型與結構方程模型相結合，并通過貝葉斯方法對模型的未知參數進行估計[9]。

隨著大數據的發展，隱馬爾可夫模型在生活中也得到了更廣泛的應用。金融領域，隱馬爾可夫模型可以用來對股市收益率波動和狀態轉換進行建模和預測[10];互聯網領域，隱馬爾可夫模型可以用來進行網絡環境監管;安全領域，隱馬爾可夫模型可以和貝葉斯神經網絡結合，通過分析一些先前發生的事件來預測恐怖分子在未來一段時間內可能發動恐怖活動的概率，用來預防可能發生的恐怖活動[11]。

目前國內外期刊發表的論文主要集中在齊次隱馬爾可夫模型的研究，即在轉移模型中只考慮了隱狀態從時刻的狀態轉移到時刻的狀態的概率，沒有考慮觀測變量自身的某些因素對狀態之間轉移概率的影響。而我們發現，在許多實際應用中，模型在時刻的狀態往往不僅僅受時刻狀態的影響，還可能會受時刻自身某些協變量的影響。為此，本文提出了一種新的轉移模型用來描述隱狀態個數已知的情況下，隱狀態之間的非齊次概率轉移關系，稱為非齊次隱馬爾可夫模型，這是本文的第一個創新點。本文的第二個創新點是提出了利用k-means聚類方法確定非齊次隱馬爾可夫模型的隱狀態，在隱狀態確定的條件下，再分別估計轉移模型和發射模型中的未知參數。

二、模型描述

以上描述的就是非齊次隱馬爾可夫模型的狀態轉移過程，我們將這種形式的轉移模型稱作“連續比對數轉移模型”。

在隱狀態確定的條件下，觀測向量為服從特定分布的向量，其中p是觀測向量的維度。即時，對非齊次隱馬爾可夫模型來說，向量為模型的觀測變量，為影響轉移效果的固定協變量，q是固定協變量的維度。本文要研究的就是如何利用可觀測到的信息和去估計非齊次隱馬爾可夫模型中的轉移模型和每個隱狀態下發射模型中的未知參數。

三、估計原理

（一）隱狀態向量

隱馬爾可夫模型的觀測變量來自不同的隱狀態，然而在實際應用中，隱狀態無法直接觀測，這也是隱馬爾可夫模型推斷的困難之處。已有的許多研究，都試圖利用觀測變量中隱含的信息首先對觀測變量的隱狀態進行判定，在確定隱狀態向量之后，再對隱馬爾可夫模型中的轉移參數和每個隱狀態下特定分布的未知參數進行估計。其中，最有代表性的就是向前向后遞歸算法。然而向前向后遞歸算法的理論性很強，計算方法非常復雜，要求使用者具有一定的統計學理論背景和計算機編程能力，這極大地限制了隱馬爾可夫模型在實際生產生活中的推廣和應用。

本文利用k-means聚類方法，來確定隱馬爾可夫模型觀測變量的隱狀態，在隱狀態確定的條件下對模型中的未知參數進行估計。

（二）標簽交換

在對非齊次隱馬爾可夫模型進行研究時，首先要解決標簽交換問題。如果沒有解決這一問題，模型的參數估計結果就不具有可解釋性，甚至會發生混淆。Richardson和 Green 在文章中通過比較均值來解決標簽交換問題[12]。本文借鑒這一經驗，在發射模型為多元正態分布時，通過比較多元正態分布中均值向量第一分量的方法來解決這一問題;在發射模型為多元線性回歸時，通過比較回歸系數的方法來解決這一問題。即在參數估計之后，重新確定觀測變量隱狀態的標簽，然后再根據隱狀態標簽結果，重新確定每個隱狀態下參數的估計結果。

（三）轉移模型中未知參數的估計

利用極大似然估計對非齊次隱馬爾可夫模型中轉移模型的未知參數進行估計。未知參數的似然函數如下：

最大似然估計就是要找到使得似然函數取最大值時未知參數和的值。本文中似然函數含有個未知參數。在實際應用中，可以用統計軟件中求極值的函數得到未知參數的數值解。比如，本文就是通過R語言BB程序包中的fun函數，近似求解似然函數的極大值點。

（四）發射模型中未知參數的估計

隱馬爾可夫模型的發射模型可以來自各種不同的特定分布。本文選擇比較經典和常用的多元正態分布和多元線性回歸作為兩個模擬實驗的發射模型。

1.多元正態分布。假設聚類分析確定隱狀態向量后，第k個隱狀態下的觀測集合為，。則觀測模型中的待估參數為每個隱狀態下的均值向量和協方差矩陣。

用極大似然估計對發射模型中的待估參數進行估計，結果如下：

2.多元線性回歸。假設聚類分析確定隱狀態向量后，第k個隱狀態下的觀測向量為。則觀測模型中的待估參數為每個隱狀態下自變量的系數向量。

用最小二乘估計對觀測模型中的待估參數進行估計，結果如下：

四、實證分析

本實例的數據來源于曲靖師范學院數學與統計學院2014級數學與應用數學專業兩個班91位同學八個學期的綜合測評成績。每名學生每學期的綜合測評成績分為文化知識成績和創新發展成績兩類，文化知識成績和創新發展成績均采用百分制計算。我們將每名學生每學期的文化知識成績和創新發展成績作為發射模型的觀測數據，將它們近似的看作一個二維的正態分布。將每名學生的性別作為固定協變量，男生設置為1，女生設置為0;將學生是否獲得獎學金作為模型的隱狀態，綜合測評分數高的同學獲得獎學金。用k-means聚類分析確定每個觀測變量的隱狀態，再利用文章中介紹的標簽交換的方法，聚類后，表示學生獲得獎學金，表示學生沒有獲得獎學金。

根據文中介紹的參數估計方法，得到各參數的估計值：

獲得獎學金的同學，綜合測評成績服從：

沒有獲得獎學金的同學，綜合測評成績服從：

根據發射模型參數估計結果可知：獲得獎學金的同學綜合測評成績明顯高于沒有獲得獎學金的同學。具體來說，獲得獎學金的同學，文化知識成績平均分為86.24，創新發展成績平均分為86.40;沒有獲得獎學金的同學，文化知識成績平均分為78.48，創新發展成績平均分為77.36。

根據隱狀態轉移概率公式可以知道，當學生性別為女時，如果前一個學期該學生獲得了獎學金，那么后一個學期該學生獲得獎學金的概率為0.80，不能獲得獎學金的概率為0.20;如果前一個學期該學生沒有獲得獎學金，那么后一個學期該學生可以獲得獎學金的概率為0.45，不能獲得獎學金的概率為0.55。當學生性別為男生時，如果前一個學期該學生獲得獎學金，那么后一個學期該學生可以獲得獎學金的概率為0.51，不能獲得獎學金的概率為0.49;如果前一個學期該學生沒有獲得獎學金，那么后一個學期該學生可以獲得獎學金的概率為0.18，不能獲得獎學金的概率為0.82。

五、結論

本文在現有的對齊次隱馬爾可夫模型研究的基礎上，通過在轉移模型的馬爾可夫鏈中增加協變量，提出了非齊次隱馬爾可夫模型，解決了當模型的觀測變量存在固定協變量時，對隱狀態轉移關系進行建模的問題，這是本文的第一個創新點。本文采用了“連續比對數轉移模型”來描述隱狀態之間的轉移關系，因為模型中未知參數較多，且不同情況下未知參數的個數也不同，所以今后可以對非齊次轉移模型進行相關的改進研究，使得轉移模型形式更簡單，使用更方便。

本文的第二個創新點是利用較簡單的k-means聚類分析確定模型觀測變量的隱狀態，在隱狀態確定的情況下對非齊次隱馬爾可夫模型中的未知參數進行估計。在今后的研究中，可以嘗試利用系統聚類方法對隱狀態個數未知的隱馬爾可夫模型進行隱狀態個數的模型選擇，或者嘗試利用神經網絡、決策樹等聚類方法來確定觀測變量的隱狀態。

參考文獻：

[1]Baum L E，Petrie T.Statistical inference for probabilistic functions of finite state Markov chains[J].Annals of Mathematical Statistics，1966，37（1）：1554-1563.

[2]Baum L E，Petrie T，Soules G，Weiss N.A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains[J].Annals of Mathematical Statistical，1970，41（1）：164-171.

[3]Jelinek，F.，Bahl，L.，Mercer，R.Design of a linguistic statistical decoder for the recognition of continuous speech[J].IEEE Transactions on Information Theory，1975 2（3）：250-256.

[4]Bishop，Martin J，Thompson，Elizabeth A.Maximum Likelihood alignment of DNA sequences[J].Journal of Molecular Biology，1986，190（2）：159-165.

[5]Dempster，A.P.Maximum likelihood from incomplete data via the EM algorithm[J].Journal of Royal Statistical Society B，1977，39（1）：1-38.

[6]Jelinek，Frederick.Statistical methods for speech recognition[M].1997.

[7]夏業茂，勾建偉，劉應安.隱馬爾可夫因子分析模型的半參數貝葉斯分析[J].高校應用數學學報，2015，30（1）：17-30.

[8]劉鶴飛，王坤，蔣成飛.隱狀態個數未知的隱馬爾可夫多元正態分布的貝葉斯推斷[J].統計研究，2017，34（12）：119-125.

[9]王坤，劉鶴飛，蔣成飛.隱馬爾可夫結構方程及其貝葉斯估計[J].數理統計與管理，2018，37（2）：273-279.

[10]劉金全，李楠，鄭挺國.隨機波動模型的馬爾可夫鏈—蒙特卡羅模擬方法—在滬市收益率序列上的應用[J].數理統計與管理，2010，29（6）：1026-1035.

[11]戰兵，韓銳.基于隱馬爾可夫的恐怖事件預測模型[J].解放軍理工大學學報，2015，16（4）：386-393.

[12]Richardson S，Green P J.On Bayesian analysis of mixtures with an unknown number of components[J].Journal of the Royal Statistical Society（Series B），1997，59（4）：731—758.

（作者單位：首都經濟貿易大學統計學院）