【机器学习】隐马尔可夫模型及其三个基本问题（四）状态序列预测算法及python实现

一、维特比算法
二、python实现
参考资料

隐马尔可夫模型状态序列预测问题是指给定模型 λ=[A,B,∏]\lambda = \left[ {A,B,\prod } \right]λ=[A,B,∏]和观测序列 X={x1,x2,⋯ ,xn}X = \left\{ {{x_1},{x_2}, \cdots ,{x_n}} \right\}X={x1,x2,⋯,xn}，求最可能出现的对应状态序列。本篇博文介绍状态序列预测算法中的维特比算法（Viterbi algorithm）[参考资料1]。

一、维特比算法

维特比算法实际是用动态规划解隐马尔科夫模型的预测问题，首先导入两个变量δ\deltaδ和ψ\psiψ。

定义在时刻ttt状态为qiq_iqi的所有单个路径(y1,y2,⋯ ,yt)({y_1},{y_2}, \cdots ,{y_t})(y1,y2,⋯,yt)中概率最大值为：
δt(qi)=max⁡y1,y2,⋯ ,yt−1P(yt=qi,yt−1,⋯ ,y1,xt,⋯ ,x1∣λ),i=1,2,⋯ ,N{\delta _t}\left( {{q_i}} \right) = \mathop {\max }\limits_{{y_1},{y_2}, \cdots ,{y_{t - 1}}} P\left( {{y_t} = {q_i},{y_{t - 1}}, \cdots ,{y_1},{x_t}, \cdots ,{x_1}\left| \lambda \right.} \right),i = 1,2, \cdots ,Nδt(qi)=y1,y2,⋯,yt−1maxP(yt=qi,yt−1,⋯,y1,xt,⋯,x1∣λ),i=1,2,⋯,N

则变量δ\deltaδ的递推公式为：
δt+1(qi)=max⁡y1,y2,⋯ ,ytP(yt+1=qi,yt,⋯ ,y1,xt+1,⋯ ,x1∣λ)=max⁡1≤j≤N[δt(j)aji]bi(xt+1),i=1,2,⋯ ,N;t=1,2,⋯ ,n−1{\delta _{t + 1}}\left( {{q_i}} \right) = \mathop {\max }\limits_{{y_1},{y_2}, \cdots ,{y_t}} P\left( {{y_{t + 1}} = {q_i},{y_t}, \cdots ,{y_1},{x_{t + 1}}, \cdots ,{x_1}\left| \lambda \right.} \right) = \mathop {\max }\limits_{1 \le j \le N} \left[ {{\delta _t}\left( j \right){a_{ji}}} \right]{b_i}({x_{t + 1}}),i = 1,2, \cdots ,N;t = 1,2, \cdots ,n - 1δt+1(qi)=y1,y2,⋯,ytmaxP(yt+1=qi,yt,⋯,y1,xt+1,⋯,x1∣λ)=1≤j≤Nmax[δt(j)aji]bi(xt+1),i=1,2,⋯,N;t=1,2,⋯,n−1

定义在时刻ttt状态为qiq_iqi的所有单个路径(y1,y2,⋯ ,yt−1,qi)({y_1},{y_2}, \cdots ,{y_{t - 1}},q_i)(y1,y2,⋯,yt−1,qi)中概率最大的路径的第t−1t-1t−1个节点为：
ψt(qi)=arg⁡max⁡1≤j≤N[δt−1(j)aji],i=1,2,⋯ ,N{\psi _t}(q_i) = \arg \mathop {\max }\limits_{1 \le j \le N} \left[ {{\delta _{t - 1}}\left( j \right){a_{ji}}} \right],i = 1,2, \cdots ,Nψt(qi)=arg1≤j≤Nmax[δt−1(j)aji],i=1,2,⋯,N

维特比算法步骤：
输入：模型λ=[A,B,∏]\lambda = [A,B,\prod ]λ=[A,B,∏]和观测序列X={x1,x2,⋯ ,xn}X = \left\{ {{x_1},{x_2}, \cdots ,{x_n}} \right\}X={x1,x2,⋯,xn}
输出：最优状态序列Y={y1,y2,⋯ ,yn}Y = \left\{ {{y_1},{y_2}, \cdots ,{y_n}} \right\}Y={y1,y2,⋯,yn}
（1）初始化：

δ1(qi)=πibi(x1),i=1,2,⋯ ,Nψ1(qi)=0,i=1,2,⋯ ,N\begin{array}{l} {\delta _1}({q_i}) = {\pi _{_i}}{b_i}\left( {{x_1}} \right),i = 1,2, \cdots ,N\\ {\psi _1}\left( {{q_i}} \right) = 0,i = 1,2, \cdots ,N \end{array}δ1(qi)=πibi(x1),i=1,2,⋯,Nψ1(qi)=0,i=1,2,⋯,N

（2）递推：对于t=2,3,⋯ ,nt = 2,3, \cdots ,nt=2,3,⋯,n

δt(qi)=max⁡1≤j≤N[δt−1(j)aji]bi(xt),i=1,2,⋯ ,Nψt(qi)=arg⁡max⁡1≤j≤N[δt−1(j)aji],i=1,2,⋯ ,N\begin{array}{l} {\delta _t}\left( {{q_i}} \right) = \mathop {\max }\limits_{1 \le j \le N} \left[ {{\delta _{t - 1}}\left( j \right){a_{ji}}} \right]{b_i}({x_t}),i = 1,2, \cdots ,N\\ {\psi _t}\left( {{q_i}} \right) = \arg \mathop {\max }\limits_{1 \le j \le N} \left[ {{\delta _{t - 1}}\left( j \right){a_{ji}}} \right],i = 1,2, \cdots ,N \end{array}δt(qi)=1≤j≤Nmax[δt−1(j)aji]bi(xt),i=1,2,⋯,Nψt(qi)=arg1≤j≤Nmax[δt−1(j)aji],i=1,2,⋯,N

（3）终止：

yn=arg⁡max⁡1≤i≤N[δn(qi)]{y_n} = \arg \mathop {\max }\limits_{1 \le i \le N} \left[ {{\delta _n}\left( {{q_i}} \right)} \right]yn=arg1≤i≤Nmax[δn(qi)]

（4）最优路径回朔：对于t=n−1,n−2,⋯ ,1t = n - 1,n - 2, \cdots ,1t=n−1,n−2,⋯,1

yt=ψt+1(yt+1){y_t} = {\psi _{t + 1}}\left( {{y_{t + 1}}} \right)yt=ψt+1(yt+1)

二、python实现

代码参考资料2

import numpy as npclass HMM:def __init__(self, A, B, Pi, O):self.A = np.array(A)     # 状态转移概率矩阵self.B = np.array(B)     # 观测概率矩阵self.Pi = np.array(Pi)   # 初始状态概率矩阵self.O = np.array(O)     # 观测序列self.N = self.A.shape[0] # 状态取值个数self.M = self.B.shape[1] # 观测值取值个数def viterbi(self):n = len(self.O)  #观测序列和状态序列的长度Y = np.zeros(n)  #初始化状态序列delta = np.zeros((n,self.N))psi = np.zeros((n,self.N))#初始化for i in range(self.N):delta[0,i] = self.Pi[i] * self.B[i,self.O[0]]psi[0,i] = 0#递推for t in range(1,n):for i in range(self.N):delta[t,i] = self.B[i,self.O[t]] * np.array([delta[t-1,j] * self.A[j,i] for j in range(self.N)]).max()psi[t,i] = np.array([delta[t-1,j] * self.A[j,i] for j in range(self.N)]).argmax()print(psi)#终止Y[n-1] = delta[n-1,:].argmax()print(Y[n-1])#回朔for t in range(n-2,-1,-1):Y[t] = psi[int(t+1),int(Y[t+1])]return Yif __name__ == '__main__':print('----------------1.初始化HMM模型参数------------------')A = [[0,1,0,0],[0.4,0,0.6,0],[0,0.4,0,0.6],[0,0,0.5,0.5]]B = [[0.5,0.5],[0.3,0.7],[0.6,0.4],[0.8,0.2]]Pi = [0.25,0.25,0.25,0.25]print('----------------2.观测序列---------------------------')O = [0,1,0]print('----------------3.viterbi 算法-----------------------')hmm = HMM(A,B,Pi,O)Y = hmm.viterbi()

参考资料

1、《统计学习方法》李航
2、https://www.cnblogs.com/d-roger/articles/5719979.html