維特比演算法和隱馬爾可夫模型的解碼

一、概述

??維特比演算法是安德魯.維特比(Andrew Viterbi)于1967年為解決通信領域中的解碼問題而提出的，它同樣廣泛用于解決自然語言處理中的解碼問題，隱馬爾可夫模型的解碼是其中典型的代表，無論是通信中的解碼問題還是自然語言處理中的解碼問題，本質上都是要在一個籬笆網路中尋找得到一條最優路徑，
??所謂籬笆網路，指的是單向無環圖，呈層級連接，各層節點數可以不同，如圖是一個籬笆網路，連線上的數字是節點間概念上的距離（如間距、代價、概率等），現要找到一條從起始點到終點的最優路徑，

??在實際問題中，節點數和層數往往是大量的，因而采取遍歷所有的路徑計算其距離進行比較的方式是不可行的，維特比演算法正是通過動態規劃的方式高效求得這條最優路徑，

二、維特比演算法

1.演算法原理

??該問題具有這樣一個特性，對于最優（如最短距離）的路徑，任意一段子路徑一定是該段兩端點間所有可達路徑中最優的，如若不然，將該段中更優的子路徑接到兩端點便構成了另一個整體最優路徑，這是矛盾的，或者說，最優路徑中，從起始點到由近及遠的任一點的子路徑，一定是該段所有可達路徑中最優的，也即，整體最優，區域一定最優，

??該特性也就是說，對每個節點而言，如果最優路徑經過這一點，則一定是經過從起始點到這點的這段最優路徑，那么，只要從頭開始，由區域向整體推進，漸次地找到起始點到當前點的最優路徑，算至終點便得到了整體最優路徑，這樣的方式叫做動態規劃，是維特比演算法的基本思想，
維特比演算法求解籬笆網路最短路徑的程序為：
??從第一層開始，對每一層的每一個節點，計算出起始點到它的最短距離，并記錄下相應最短路徑下它的前一個節點，逐層遞推，算至終止點時便得到了整體最短距離，再依照節點記錄下的前置節點進行回溯，就得到了最短路徑的序列，對第\(i\)層第\(j\)個節點\(P_{i,j}\)，假設起始點到它的最短距離為\(D\left( P_{i,j} \right)\)，相應最短路徑下它的前一個節點為\(Pre\left( P_{i,j} \right)\)，則

也就是，對前一層的所有節點，計算每一個節點的記錄的最短距離D與它到當前節點的距離d的和，取其中最小的那個值（其中， \(d(A,B)\)表示A,B兩節點間的距離，）.

也就是，滿足距離最短的那條路徑上在前一層的節點.

2.示例

??在如下網路中，連線上的數字是節點間的距離，求S點到E點的最短距離和與之對應的路徑，

第一層：
對節點\(P_{1,1}\)，
起始點到它只有一條路徑，最短距離\(D(P_{1,1})=2\)，前一個節點\(Pre(P_{1,1}) = S\)；

對節點\(P_{1,2}\)，
起始點到它只有一條路徑，最短距離 \(D(P_{1,2}) = 1\)，前一個節點\(Pre(P_{1,2}) = S\)；

對節點\(P_{1,3}\)，
起始點到它只有一條路徑，最短距離\(D(P_{1,3}) = 3\)，前一個節點\(Pre(P_{1,3}) = S\) ，

第二層：
對節點 \(P_{2,1}\)，
\(D(P_{1,1}) + d(P_{1,1}，P_{2,1}) = 2 + 3 = 5\)，
\(D(P_{1,2}) + d(P_{1,2}，P_{2,1}) = 1 + 2 = 3\) ，
\(D(P_{1,3}) + d(P_{1,3}，P_{2,1}) = 3 + 9 = 12\)，
最短距離 \(D(P_{2,1}) = min\left\{ 5,3,12 \right\} = 3\)，前一個節點\(Pre(P_{2,1}) = P_{1,2}\) ;

對節點\(P_{2,2}\) ，
\(D(P_{1,1}) + d(P_{1,1}，P_{2,2}) = 2 + 6 = 8\)，
\(D(P_{1,2}) + d(P_{1,2}，P_{2,2}) = 1 + 5 = 6\)，
\(D(P_{1,3}) + d(P_{1,3}，P_{2,2}) = 3 + 2 = 5\)，
最短距離\(D(P_{2,2}) = min\left\{ 8,6,5 \right\} = 5\)，前一個節點\(Pre(P_{2,2}) = P_{1,3}\) ;

對節點\(P_{2,3}\)，
\(D(P_{1,1}) + d(P_{1,1}，P_{2,3}) = 2 + 4 = 6\) ，
\(D(P_{1,2}) + d(P_{1,2}，P_{2,3}) = 1 + 7 = 8\)，
\(D(P_{1,3}) + d(P_{1,3}，P_{2,3}) = 3 + 6 = 9\) ，
最短距離\(D(P_{2,3}) = min\left\{ 6,8,9 \right\} = 6\)，前一個節點\(Pre(P_{2,3}) = P_{1,1}\) ;

第三層：
對節點 \(P_{3,1}\)，
\(D(P_{2,1}) + d(P_{2,1}，P_{3,1}) = 3 + 9 = 12\)，
\(D(P_{2,2}) + d(P_{2,2}，P_{3,1}) = 5 + 2 = 7\)，
\(D(P_{2,3}) + d(P_{2,3}，P_{3,1}) = 6 + 6 = 12\)，
最短距離\(D(P_{3,1}) = min\left\{ 12,7,12 \right\} = 7\)，前一個節點\(Pre(P_{3,1}) = P_{2,2}\);

對節點\(P_{3,2}\)，
\(D(P_{2,1}) + d(P_{2,1}，P_{3,2}) = 3 + 3 = 6\)，
\(D(P_{2,2}) + d(P_{2,2}，P_{3,2}) = 5 + 6 = 11\)，
\(D(P_{2,3}) + d(P_{2,3}，P_{3,2}) = 6 + 3 = 9\)，
最短距離\(D(P_{3,2}) = min\left\{ 6,11,9 \right\} = 6\)，前一個節點\(Pre(P_{3,2}) = P_{2,1}\);

對節點\(P_{3,3}\)，
\(D(P_{2,1}) + d(P_{2,1}，P_{3,3}) = 3 + 8 = 11\)，
\(D(P_{2,2}) + d(P_{2,2}，P_{3,3}) = 5 + 7 = 12\)，
\(D(P_{2,3}) + d(P_{2,3}，P_{3,3}) = 6 + 4 = 10\)，
最短距離\(D(P_{3,3}) = min\left\{ 11,12,10 \right\} = 10\)，前一個節點\(Pre(P_{3,3}) = P_{2,3}\);

第四層（終點）：
對節點 \(E\)，
\(D(P_{3,1}) + d(P_{3,1}，E) = 7 + 3 = 10\)，
\(D(P_{3,2}) + d(P_{3,2}，E) = 6 + 7 = 13\)，
\(D(P_{3,3}) + d(P_{3,3}，E) = 10 + 6 = 16\)，
最短距離\(D(E) = min\left\{ 10,13,16 \right\} = 10\)，前一個節點\(Pre(E) = P_{3,1}\);

又\(Pre(E) = P_{3,1}\)，\(Pre(P_{3,1}) = P_{2,2}\)，\(Pre(P_{2,2}) = P_{1,3}\)，\(Pre(P_{1,3}) = S\)
故最短距離為10，與之對應的路徑為（\(S\)，\(P_{1,3}\)，\(P_{2,2}\)，\(P_{3,1}\)，\(E\)）.

三、隱馬爾可夫模型的解碼

1.問題描述

??隱馬爾可夫模型（HMM）的解碼問題指，給定模型和輸出序列，如何找出最有可能產生這個輸出的狀態序列，自然語言處理中，也即如何通過觀測信號確定最有可能對應的實際語意，在狀態序列上，每個狀態位是狀態集合中的元素之一，因此該問題等價于在狀態集合中的節點構成的有向網路（籬笆網路）中找出一條概率最大的路徑（最優路徑），如圖，該問題可以通過維特比演算法得到高效的解決，

2.演算法敘述

??假設 \(P(s_{t,j})\)表示從起始時刻到\(s_{t,j}\)的最優路徑的概率，\(Pre(s_{t,j})\)表示從起始時刻到 \(s_{t,j}\)的最優路徑上前一個節點，則隱馬爾可夫模型的維特比解碼演算法為：

輸入：隱馬爾可夫模型 \(\lambda=(\pi,A,B)\)和觀測 \(O=(o_1,o_2,...,o_T)\)；
輸出：最優狀態序列\(S^{\ast}=(s_{1}^{\ast},s_{2}^{\ast},...,s_{T}^{\ast})\).
(1)初始化
????\(P(s_{1,j})=\pi_{j}b_{j}(o_1)\)，
????\(Pre(s_{1,j})=None\)，\(j=1,2,...,N\)

(2)遞推
?對 \(t=2,3,...,T\)

，\(j=1,2,...,N\).

(3)遞推終止
?最大概率$$P^{\ast}=\max_{1\leq j \leq N}{P(s_{T,j})}$$
?最優路徑上的最后一個狀態$$s_{T}^{\ast}=arg\max_{1\leq j \leq N}{\left[ P(s_{T,j}) \right]}$$

(4)回溯路徑，確定最優狀態序列
????\(S^{\ast}=\left( s_{1}^{\ast},s_{2}^{\ast},...,s_{T-1}^{\ast},s_{T}^{\ast} \right)\)
?????\(=\left( Pre(s_{2}^{\ast}),Pre(s_{3}^{\ast}), ...,Pre(s_{T}^{\ast}),s_{T}^{\ast}\right)\)

3.示例

（參考自《統計學習方法》）
狀態集合 \(Q=\left\{ q_1, q_2, q_3 \right\}\)，觀測集合\(V=\left\{ 0,1 \right\}\)，模型 \(\lambda=\left( \pi,A,B \right)\) ，

\(A=\begin{bmatrix} 0.5 & 0.2 & 0.3 \\ 0.3 & 0.5 & 0.2 \\ 0.2 & 0.3 & 0.5 \\ \end{bmatrix}\) , \(B=\begin{bmatrix} 0.5 & 0.5 \\ 0.4 & 0.6 \\ 0.7 & 0.3 \end{bmatrix}\),\(\pi=\left( 0.2, 0.4, 0.4 \right)^{T}\)

已知觀測序列\(O=\left( 0, 1, 0 \right)\)，求最優狀態序列，

解：
(1)在t=1時（初始化），對每一個狀態，求觀測為0的最大概率
?\(P(s_{1,1})=0.2\times0.5=0.1\)，\(Pre(s_{1,1})=None\)
?\(P(s_{1,2})=0.4\times0.4=0.16\)，\(Pre(s_{1,2})=None\)
?\(P(s_{1,3})=0.4\times0.7=0.28\)，\(Pre(s_{1,3})=None\)

(2)在t=2時，對每一個狀態，求觀測為1的
?最大概率$$P(s_{2,j})=\max_{1 \leq k \leq 3}{\left[ P(s_{1,k})a_{kj} \right]b_{j}(1)}$$
?當前最優的前一個狀態$$Pre(s_{2,j})=arg\max_{1 \leq k \leq 3}{\left[ P(s_{1,k})a_{kj} \right]}$$，\(j=1,2,3.\)
\(P(s_{2,1})=max\left\{ 0.1\times0.5\times0.5, 0.16\times0.3\times0.5, 0.28\times0.2\times0.5 \right\}\)\(=0.028\)

\(Pre(s_{2,1})=s_{1,3}=q_3\)

\(P(s_{2,2})=max\left\{ 0.1\times0.2\times0.6,0.16\times0.5\times0.6,0.28\times0.3\times0.6 \right\}\)\(=0.0504\)

\(Pre(s_{2,2})=s_{1,3}=q_3\)

\(P(s_{2,3})=max\left\{ 0.1\times0.3\times0.3, 0.16\times0.2\times0.3,0.28\times0.5\times0.3 \right\}\)\(=0.042\)

\(Pre(s_{2,3})=s_{1,3}=q_3\)

(3)在t=3時，對每一個狀態，求觀測為0的
?最大概率 $$P(s_{3,j})=\max_{1 \leq k \leq 3}{\left[ P(s_{2,k})a_{kj} \right]b_{j}(0)}$$
?當前最優的前一個狀態$$Pre(s_{3,j})=arg\max_{1 \leq k \leq 3}\left[ P(s_{2,k})a_{kj} \right]$$，\(j=1,2,3.\)

\(P(s_{3,1})=max\left\{ 0.028\times0.5\times0.5, 0.0504\times0.3\times0.5,0.042\times0.2\times0.5 \right\}\)\(=0.00756\)

\(Pre(s_{3,1})=s_{2,2}=q_2\)

\(P(s_{3,2})=max\left\{ 0.028\times0.2\times0.4, 0.0504\times0.5\times0.4, 0.042\times0.3\times0.4 \right\}\)\(=0.01008\)

\(Pre(s_{3,2})=s_{2,2}=q_2\)

\(P(s_{3,3})=max\left\{ 0.028\times0.3\times0.7,0.0504\times0.2\times0.7,0.042\times0.5\times0.7 \right\}\)\(=0.0147\)

\(Pre(s_{3,3})=s_{2,3}=q_3\)

(4)得到結果.
?最大概率
??\(P^{\ast}=\max_{1 \leq j \leq 3}P\left( s_{3,j} \right)\)
???\(=max\left\{ 0.00756,0.01008,0.0147 \right\}\)
???\(=0.0147\)
?最優狀態序列
??\(S^{\ast}=\left( Pre(s_{2}^{\ast})t,Pre(s_{3}^{\ast}),s_{3}^{\ast} \right)\)
???\(=\left( s_{1,3},s_{2,3},s_{3,3} \right)\)
???\(=\left( q_3,q_3,q_3 \right)\)

4.python實作

對上述HMM解碼示例的python實作程式為

import numpy as np

def viterbi(pi, A, B, Q, V, obs_seq):
    '''
    :param pi:HMM初始狀態概率向量，list型別
    :param A:HMM狀態轉移概率矩陣，list型別
    :param B:HMM觀測生成概率矩陣，list型別
    :param Q:狀態集合，list型別
    :param V:觀測集合，list型別
    :param obs_seq:觀測序列，list型別
    :return:最優狀態序列的概率sta_pro，float型別；最優狀態序列sta_seq，list型別
    '''

    # HMM模型引數轉換為array型別
    pi = np.array(pi)
    A = np.array(A)
    B = np.array(B)

    # 1.定義動態計算結果存盤矩陣
    rowNum = len(Q)  # 行數，狀態數
    colNum = len(obs_seq)  # 列數，生成的觀測數，即時刻數

    # 存盤節點當前最大概率的矩陣
    probaMatrix = np.zeros((rowNum,colNum))

    # 存盤當前最優路徑下的前一個節點的矩陣
    preNodeMatrix = np.zeros((rowNum,colNum))

    # 2.初始化（第1時刻）
    probaMatrix[:,0] = pi*np.transpose(B[:,obs_seq[0]])
    preNodeMatrix[:,0] = [-1]*rowNum  # 第1時刻節點的前一個節點不存在，置為-1

    # 3.遞推，第2時刻至最后
    for t in range(1, colNum):
        list_pre_max_proba = []  # 節點最大前置概率串列
        list_pre_node = []  # 節點當前最優路徑中前一個節點串列
        for j in range(rowNum):
            pre_proba_list = list(np.array(probaMatrix[:,t-1])*np.transpose(A[:,j]))  # 前置概率串列，前一時刻的節點最大概率與到當前節點轉移概率的乘積
            '''
            注：因為計算機的二進制機制對小數的表達是有限的，所以對小數作運算將產生一定的誤差，
            在使用函式獲取pre_proba_list中的最大值和對應的索引時，為有效降低這種誤差，將資料放大后再進行操作，
            '''
            pre_proba_list = [x*pow(10,5) for x in pre_proba_list]  # 放大100000倍
            prePro = max(pre_proba_list)/pow(10,5)  # 最大前置概率
            preNodeIndexNum = pre_proba_list.index(max(pre_proba_list))  # 前置節點的索引號
            list_pre_max_proba.append(prePro)  # 最大前置概率加入串列
            list_pre_node.append(preNodeIndexNum)  # 前置節點的索引號加入串列

        probaMatrix[:,t] = np.array(list_pre_max_proba)*np.transpose(B[:,obs_seq[t]])  # 最大前置概率乘上觀測概率，即為當前最大概率
        preNodeMatrix[:,t] = list_pre_node  # 將該列前置節點索引號加入矩陣

    # 此時，得到了完整的probaMatrix和preNodeMatrix，對這兩個矩陣進行操作便可得到需要的結果
    # 4.得到最大概率
    maxPro = np.max(probaMatrix[:, colNum-1])  # 全域最大概率（即最后一列的最大值）

    # 5.得到最優狀態序列的狀態索引號串列
    lastStateIndexNum = np.argmax(probaMatrix[:, colNum-1])  # 最優狀態序列中最后一個狀態的索引號
    stateIndexList = []  # 定義最優狀態的索引號串列
    stateIndexList.append(lastStateIndexNum)

    # 回溯，完成狀態索引號串列
    currentIndex = lastStateIndexNum;
    for t in range(colNum-1, 0, -1):
        fls = preNodeMatrix[:, t].tolist()  # 矩陣中的數值是浮點型
        ls = list(map(int, fls))  # 轉為整型
        currentIndex = ls[currentIndex]
        stateIndexList.append(currentIndex)

    stateIndexList.reverse()  # 反轉串列

    # 6.由索引號序列得到最優狀態序列
    stateSeq = [Q[i] for i in stateIndexList]

    return maxPro,stateSeq

if __name__=='__main__':
    # 狀態集合
    Q = ["q1", "q2", "q3"]

    # 觀測集合
    V = [0, 1]

    # 初始狀態概率向量
    pi = [0.2, 0.4, 0.4]

    # 狀態轉移概率矩陣
    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]]

    # 觀測序列
    obs_seq = [0, 1, 0]

    maxPro, stateSeq = viterbi(pi, A, B, Q, V, obs_seq)

    print("最大概率為：", maxPro)
    print("最優狀態序列為：", stateSeq)

End.

參考

1. 吳軍. 數學之美(第二版). 人民郵電出版社.
2. 李航. 統計學習方法. 清華大學出版社.
3. https://www.cnblogs.com/zhibei/p/9391014.html

標籤：其他

上一篇：TCP/IP概述（網路互聯與TCP/IP）

下一篇：2021年在北京建一個網站一般需要多少錢呢？