HMM算法-解碼問題

這篇文章記錄一下解決HMM三大問題的第二個問題的學習過程。回憶一下，第二個問題是什麼來着？給定HMM模型\(lambda\)和觀測序列O，求產生這個觀測序列機率最大的狀態序列是什麼？把這個問題叫作解碼問題，也是挺貼切的~

求解這個問題，有一個經典的算法，叫作Viterbi算法。Viterbi是個了不得的人物，數學之美第26就是講Viterbi和他的Viterbi算法。

Viterbi算法針對籬笆網絡有向圖（Lattice）的最短路徑問題提出，是一個特殊可是應用最普遍的動態規劃算法。凡是使用HMM的問題均可以用它來解碼，包括數字通訊、語音識別、機器翻譯、拼音轉漢字、分詞等。

Viterbi算法相似求解觀測序列機率時使用的forward方法，它也定義了一個變量:

\(\delta=\displaystyle\max_{q_1...q_{t-1}}P(q_1q_2...q_t=s_i,O_1O_2..O_t|\lambda)\)。

這是t時刻狀態\(q_t=s_i\)時的最優狀態序列和前t個觀察序列的聯合機率。\(\delta_i(t)\)和\(\delta_j(t+1)\)的關係是：

\(\delta_j(t+1)=[\displaystyle\max_{1 \leq i \leq N}\delta_i(t)a_{ij}]b_j(O_{t+1})\)

這個和forward方法很是的相似，forward方法是從t轉移到t+1時的全部N個可能的機率加和，而vebiter是這N個狀態的求最大。

整個解碼過程能夠歸納爲：

 

雖然如今看來維特比算法並非很複雜，可是當時提出來但是一件很是了不得的事！有些真理就是一旦發現就是如此簡單，可是發現它，可能要窮盡幾代人的努力。

下面用viterbi解決一個小的問題，問題是這樣的:假設手裏有三個不一樣的骰子。第一個骰子是咱們日常見的骰子(稱這個骰子爲D6), 6個面,每一個面(1,2,3,4,5,6)出現的機率是1/6。第二個骰子是個四面體(稱 這個骰子爲D4),每一個面(1,2,3,4)出現的機率是1/4。第三個骰子有八個面 (稱這個骰子爲D8),每一個面(1,2,3,4,5,6,7,8)出現的機率是1/8。 假設咱們開始擲骰子,咱們先從三個骰子裏挑一個,挑到每個骰子的機率都是1/3。 而後咱們擲骰子,獲得一個數字,1,2,3,4,5,6,7,8中的一個。不停的重複上 述過程,咱們會獲得一串數字,每一個數字都是1,2,3,4,5,6,7,8中的一個。例 如咱們可能獲得這麼一串數字(擲骰子10次):1 6 3 5 2 7 3 5 2 4 。如今求每次拋出去的是哪一個骰子？

python 代碼以下

def verbiter(S,K,A,B,pi,Obv):
    '''
    求解出現Obv機率最大的狀態序列
    :param S: array,狀態符號集合
    :param K: array,觀測符號集合
    :param A: matrix,轉移矩陣
    :param B: matrix,發射矩陣
    :param pi: 初始矩陣
    :param Obv: 觀測序列
    :return: 狀態序列
    '''
    N = len(S)
    M = len(K)
    Delta = []
    #初始化
    start_node = K.index(Obv[0])
    fst_row = [{"prob":pi[i]*B[i][start_node],"pre_node":-1} for i in range(0,N)]
    Delta.append(fst_row)


    for t in range(1,M):
        row = []
        for j in range(0,N):
            ob = K.index(Obv[t])
            ob_prop = B[j][ob]
            trans = [Delta[t-1][i]["prob"]*A[i][j] for i in range(0,N)]
            maxindex, maxvalue = max(enumerate(trans), key=lambda x: x[1])
            row.append({"prob":maxvalue*ob_prop,"pre_node":maxindex})
        Delta.append(row)

    #T時刻最大機率狀態
    max_s,max_prob = max(enumerate(Delta[M-1]),key=lambda x:x[1]["prob"])

    #最優路徑
    cur_node = max_s
    path = [max_s]
    for t in range(M-1,1,-1):
        pre_node = Delta[t][cur_node]["pre_node"]
        path.append(pre_node)
        cur_node = pre_node

    path.reverse()
    print ','.join(S[i] for i in path)
    print max_prob["prob"]



if __name__ == "__main__":
    S = ["D6","D4","D8"]
    K = [1,2,3,4,5,6,7,8]
    A = [[1/3.0,1/3.0,1/3.0],
         [1/3.0,1/3.0,1/3.0],
         [1/3.0,1/3.0,1/3.0]]
    B = [[1/6.0,1/6.0,1/6.0,1/6.0,1/6.0,1/6.0,0,0],
         [1/4.0,1/4.0,1/4.0,1/4.0,0,0,0,0],
         [1/8.0,1/8.0,1/8.0,1/8.0,1/8.0,1/8.0,1/8.0,1/8.0]]
    pi = [1/3.0,1/3.0,1/3.0]
    Obv = [1,6,3,5,2,7,3,5,2,4]
    verbiter(S,K,A,B,pi,Obv)