这篇文章记录一下解决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)