matlab中的隱馬爾可夫模型(HMM)實現

原文連接：http://tecdat.cn/?p=7960

隱馬爾可夫模型（HMM）簡介

_隱馬爾可夫模型_（HMM）是一個在你觀察到的輸出順序，但不知道狀態序列模型產生輸出的過程。隱馬爾可夫模型的分析試圖從觀察到的數據中恢復狀態序列。

例如，考慮具備兩個狀態和六個可能輸出的馬爾可夫模型。該模型使用：

• 紅色骰子，有六個面，標記爲1到6。
• 一個綠色骰子，具備十二個側面，其中五個側面標記爲2到6，其他七個側面標記爲1。
• 加權的紅色硬幣，正面出現機率爲.9，背面出現機率爲1.。
• 加權綠色硬幣，其正面機率爲0.95，背面機率爲.05。

該模型使用如下規則從集合{一、二、三、四、五、6}中建立數字序列：

• 首先滾動紅色骰子，而後寫下出現的數字 。

• 投擲紅色硬幣並執行如下操做之一：

  • 若是結果爲正面，則滾動紅色骰子並記下結果。
  • 若是結果是反面，則滾動綠色骰子並記下結果。
• 在隨後的每一個步驟中，您翻轉與上一步中滾動的骰子顏色相同的顏色的硬幣。若是硬幣正面朝上，則與上一步驟滾動相同的骰子。若是硬幣出現反面，請切換到另外一個骰子。

該模型的狀態圖具備紅色和綠色兩種狀態，以下圖所示。

實現")

您能夠經過滾動具備與狀態相同顏色的骰子來肯定狀態的發射。您能夠經過翻轉與狀態相同顏色的硬幣來肯定到下一個狀態的過渡。

轉換矩陣爲：

T = [0.90.050.10.95]

輸出矩陣爲：

實現")

該模型不是隱藏的，由於您能夠從硬幣和骰子的顏色知道狀態的順序。可是，假設其餘人 沒有向您顯示骰子或硬幣。您所看到的只是輸出的順序。若是開始看到的數字比其餘數字多1，則可能會懷疑骰子處於綠色狀態，但因爲沒法看到要滾動的骰子的顏色，所以沒法肯定。

隱藏的馬爾可夫模型提出如下問題：

• 給定一系列輸出，最可能的狀態路徑是什麼？
• 給定一系列輸出，您如何估算模型的轉換和輸出機率？
• 什麼是_後驗機率_？

分析隱馬爾可夫模型

本節說明如何來分析隱馬爾可夫模型。

生成測試序列

TRANS = [.9 .1; .05 .95]; EMIS = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6;... 7/12, 1/12, 1/12, 1/12, 1/12, 1/12];

要從模型生成狀態和發射的隨機序列 ：

輸出seq是序列，輸出states是狀態序列。

估計狀態序列

likelystates是與長度相同的序列seq。

要測試的準確性hmmviterbi 。

sum(states==likelystates)/1000 ans = 0.8200

在這種狀況下，最有可能的狀態序列在82％的時間內與隨機序列一致。

===

估算轉移和輸出矩陣

 返回轉換矩陣和輸出矩陣的估計值：

您能夠將輸出與原始 矩陣進行比較， TRANS而且EMIS：

TRANS TRANS = 0.9000 0.1000 0.0500 0.9500 EMIS EMIS = 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.5833 0.0833 0.0833 0.0833 0.0833 0.0833

假設您對TRANS和 有如下初步猜想EMIS。

TRANS_GUESS = [.85 .15; .1 .9]; EMIS_GUESS = [.17 .16 .17 .16 .17 .17;.6 .08 .08 .08 .08 08];

您估計TRANS並EMIS以下：

TRANS_EST2 = 0.2286 0.7714 0.0032 0.9968 EMIS_EST2 = 0.1436 0.2348 0.1837 0.1963 0.2350 0.0066 0.4355 0.1089 0.1144 0.1082 0.1109 0.1220

 若是算法在最大迭代次數（默認值爲）內未能達到此容差100，則算法將暫停。 

若是算法未能達到所需的容差，請使用如下命令增長最大迭代次數的默認值：

其中，maxiter是算法執行的最大步驟數。

估計後驗狀態機率

輸出PSTATES爲_M_ × _L_矩陣，其中_M_爲狀態數，_L_爲的長度seq。PSTATES(i,j)是條件機率，該模型處於狀態i時，它產生j的 seq給出的是，seq 。

要返回序列機率的對數seq，請使用第二個輸出參數hmmdecode：

隨着序列長度的增長，序列的機率趨於0 。

更改初始狀態分佈

默認狀況下， 隱藏的Markov模型函數從狀態1開始。換句話說，初始狀態的分佈將其全部機率質量都集中在狀態1處。要分配不一樣的機率分佈，_p_ = [ _p _1，_p _2，...，_p __M_ ]，到_M個_初始狀態，執行如下操做：

若是轉換矩陣和發射矩陣分別爲TRANS和EMIS，則可使用如下命令來建立加強矩陣：

TRANS_HAT = [0 p; zeros(size(TRANS,1),1) TRANS]; EMIS_HAT = [zeros(1,size(EMIS,2)); EMIS];