HMM-維特比算法理解與實現（python）

時間 2020-05-14

標籤 hmm 維特算法理解實現 python 欄目 Python 简体版

原文原文鏈接

HMM-前向後向算法理解與實現（python）
HMM-維特比算法理解與實現（python）

解碼問題

給定觀測序列 $O=O_1O_2...O_T$，模型 $\lambda (A,B,\pi)$，找到最可能的狀態序列 $I^∗=\{i^∗_1,i^∗_2,...i^∗_T\}$

近似算法

在每一個時刻 $t$ 選擇最可能的狀態，獲得對應的狀態序列

根據HMM-前向後向算法計算時刻 $t$ 處於狀態 $i^*_t$ 的機率：html

\[i^∗_t=argmax[\gamma_t(i)],t=1,2,...T\\ \gamma_t(i) = \frac{\alpha_{i}(t) \beta_{i}(t)}{\sum_{i=1}^{N} \alpha_{i}(t) \beta_{i}(t)} \]

可是沒法保證獲得的解是全局最優解java

維特比算法

維特比算法的基礎能夠歸納爲下面三點（來源於吳軍：數學之美）：python

若是機率最大的路徑通過籬笆網絡的某點，則從起始點到該點的子路徑也必定是從開始到該點路徑中機率最大的。算法
假定第 t 時刻有 k 個狀態，從開始到 t 時刻的 k 個狀態有 k 條最短路徑，而最終的最短路徑必然通過其中的一條。網絡
根據上述性質，在計算第 t+1 時刻的最短路徑時，只須要考慮從開始到當前的k個狀態值的最短路徑和當前狀態值到第 t+1 時刻的最短路徑便可。如求t=3時的最短路徑，等於求t=2時,從起點到當前時刻的全部狀態結點的最短路徑加上t=2到t=3的各節點的最短路徑。spa

通俗理解維特比算法，對上面三點加深理解.net

假如你從S和E之間找一條最短的路徑，最簡單的方法就是列出全部可能的路徑 ($O(T^N)$)，選出最小的，顯然時間複雜度過高。怎麼辦？（摘自[3]）code

使用維特比算法htm

S到A列的路徑有三種可能：S-A1,S-A2,S-A3，以下圖blog

S-A1,S-A2,S-A3 中一定有一個屬於全局最短路徑。繼續往右，到了B列

對B1：

會產生3條路徑：

S-A1-B1，S-A2-B1，S-A3-B1

假設S-A3-B1是最短的一條，刪掉其餘兩條。獲得

對B2:

會產生3條路徑：

S-A1-B2，S-A2-B2，S-A3-B2

假設S-A1-B2是最短的一條，刪掉其餘兩條。獲得

對B3：

會產生3條路徑：

S-A1-B3，S-A2-B3，S-A3-B3

假設S-A2-B3是最短的一條，刪掉其餘兩條。獲得

如今咱們看看對B列的每一個節點有哪些，回顧維特比算法第二點

假定第 t 時刻有 k 個狀態，從開始到 t 時刻的 k 個狀態有 k 條最短路徑，而最終的最短路徑必然通過其中的一條

B列有三個節點，因此會有三條最短路徑，最終的最短路徑必定會通過其中一條。以下圖

同理，對C列，會獲得三條最短路徑，以下圖

到目前爲止，仍然沒法肯定哪條屬於全局最短。最後，咱們繼續看E節點

最終發現最短路徑爲S-A1-B2-C3-E

數學描述

在上述過程當中，對每一列（每一個時刻）會獲得對應狀態數的最短路徑。在數學上如何表達？記錄路徑的最大機率值 $ \delta_t(i)$ 和對應路徑通過的節點 $\psi_t(i)$。

定義在時刻 $t$ 狀態爲 $i$ 的全部單條路徑中機率最大值爲

\[\delta_{t}(i)=\max _{i_{1}, i_{2}, \ldots, i_{t-1}} P\left(i_{t}=i, i_{t-1}, \ldots, i_{1}, o_{t}, \ldots, o_{1} | \lambda\right), i=1,2, \ldots, N \]

遞推公式

\[\begin{aligned} \delta_{t+1}(i) &=\max _{i_{1}, i_{2}, \ldots, i_{t}} P\left(i_{t+1}=i, i_{t}, \ldots, i_{1}, o_{t+1}, \ldots, o_{1} | \lambda\right) \\ &=\max _{1 \leq j \leq N}\left[\delta_{t}(j) a_{j i}\right] b_{i}\left(o_{t+1}\right), i=1,2, \ldots, N ; t=1,2, \ldots, T-1 \end{aligned} \]

定義在時刻 $t$ 狀態爲 $i$ 的全部單條路徑中，機率最大路徑的第 $t - 1$ 個節點爲

\[\psi_{t}(i)=\arg \max _{1 \leq j \leq N}\left[\delta_{t-1}(j) a_{j i}\right], i=1,2, \ldots, N \]

維特比算法步驟：

step1：初始化

\[\begin{aligned}&\delta_{1}(i)=\pi_{i} b_{i}\left(o_{1}\right), i=1,2, \ldots, N\\&\psi_{1}(i)=0, i=1,2, \ldots, N\\\end{aligned} \]

step2：遞推，對 $t=2,3,...,T$

\[\delta_{t}(i)=\max _{1 \leq j \leq N}\left[\delta_{t-1}(j) a_{j i}\right] b_{i}\left(o_{t}\right), i=1,2, \ldots, N \\\psi_{t}(i)=\arg \max _{1 \leq j \leq N}\left[\delta_{t-1}(j) a_{j i}\right], i=1,2, \ldots, N \\ \]

step3：計算時刻 $T$ 最大的 $ \delta _T(i)$,即爲最可能隱藏狀態序列出現的機率。計算時刻$T$最大的$\psi_T(i)$,即爲時刻$T$最可能的隱藏狀態。

\[P^{*}=\max _{1 \leq i \leq N} \delta_{T}(i) \quad i_{T}^{*}=\arg \max _{1 \leq i \leq N} \delta_{T}(i) \]

step4：最優路徑回溯，對$t=T-1,...,1$

\[i_{t}^{*}=\psi_{t+1}\left(i_{t+1}^{*}\right)\\I^*=(i_{1}^{*},i_{2}^{*},...,i_{T}^{*}) \]

代碼實現

假設從三個袋子 {1,2,3}中取出 4 個球 O={red,white,red,white}，模型參數$\lambda = (A,B,\pi)$ 以下，計算狀態序列，即取出的球來自哪一個袋子

#狀態 1 2 3
A = [[0.5,0.2,0.3],
	 [0.3,0.5,0.2],
	 [0.2,0.3,0.5]]

pi = [0.2,0.4,0.4]

# red white
B = [[0.5,0.5],
	 [0.4,0.6],
	 [0.7,0.3]]

def hmm_viterbi(A,B,pi,O):
    T = len(O)
    N = len(A[0])
    
    delta = [[0]*N for _ in range(T)]
    psi = [[0]*N for _ in range(T)]
    
    #step1: init
    for i in range(N):
        delta[0][i] = pi[i]*B[i][O[0]]
        psi[0][i] = 0
        
    #step2: iter
    for t in range(1,T):
        for i in range(N):
            temp,maxindex = 0,0
            for j in range(N):
                res = delta[t-1][j]*A[j][i]
                if res>temp:
                    temp = res
                    maxindex = j

            delta[t][i] = temp*B[i][O[t]]#delta
            psi[t][i] = maxindex

    #step3: end
    p = max(delta[-1])
    for i in range(N):
        if delta[-1][i] == p:
            i_T = i

    #step4：backtrack
    path = [0]*T
    i_t = i_T
    for t in reversed(range(T-1)):
        i_t = psi[t+1][i_t]
        path[t] = i_t
    path[-1] = i_T
    
    return delta,psi,path

A = [[0.5,0.2,0.3],[0.3,0.5,0.2],[0.2,0.3,0.5]]
B = [[0.5,0.5],[0.4,0.6],[0.7,0.3]]
pi = [0.2,0.4,0.4]
O = [0,1,0,1]
hmm_viterbi(A,B,pi,O)

結果

references：

[1]https://www.cnblogs.com/kaituorensheng/archive/2012/12/04/2802140.html

[2] https://blog.csdn.net/hudashi/java/article/details/87875259

[3] https://www.zhihu.com/question/20136144