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)