隐马尔可夫模型的解码

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$
$$P(s_{t,j})=\underset{1\le k\le N}{\max }\left[P(s_{t-1,k})a_{kj}\right]b_j(o_t)$$
$$Pre(s_{t,j})=arg\underset{1\le k\le N}{\max }\left[P(s_{t-1,k})a_{kj}\right]$$，$j=1,2,...,N$.

(3)递推终止
 最大概率$$P^{\ast }=\underset{1\le j\le N}{\max }P(s_{T,j})$$
 最优路径上的最后一个状态$$s_T^{\ast }=arg\underset{1\le j\le N}{\max }\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={q1,q2,q3}Q=\left\{ q_1, q_2, q_3 \right\}Q={q1​,q2​,q3​}，观测集合V={0,1}V=\left\{ 0,1 \right\}V={0,1}，模型 $\lambda =\left(\pi ,A,B\right)$ ，

$A=\left[\begin{array}{ccc}0.5& 0.2& 0.3\\ 0.3& 0.5& 0.2\\ 0.2& 0.3& 0.5\end{array}\right]$ , $B=\left[\begin{array}{cc}0.5& 0.5\\ 0.4& 0.6\\ 0.7& 0.3\end{array}\right]$,$\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\times 0.5=0.1$，$Pre(s_{1,1})=None$
 $P(s_{1,2})=0.4\times 0.4=0.16$，$Pre(s_{1,2})=None$
 $P(s_{1,3})=0.4\times 0.7=0.28$，$Pre(s_{1,3})=None$

(2)在t=2时，对每一个状态，求观测为1的
 最大概率$$P(s_{2,j})=\underset{1\le k\le 3}{\max }\left[P(s_{1,k})a_{kj}\right]b_j(1)$$
 当前最优的前一个状态$$Pre(s_{2,j})=arg\underset{1\le k\le 3}{\max }\left[P(s_{1,k})a_{kj}\right]$$，$j=1,2,3.$
P(s2,1)=max{0.1×0.5×0.5,0.16×0.3×0.5,0.28×0.2×0.5}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\}P(s2,1​)=max{0.1×0.5×0.5,0.16×0.3×0.5,0.28×0.2×0.5}$=0.028$

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

P(s2,2)=max{0.1×0.2×0.6,0.16×0.5×0.6,0.28×0.3×0.6}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\}P(s2,2​)=max{0.1×0.2×0.6,0.16×0.5×0.6,0.28×0.3×0.6}$=0.0504$

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

P(s2,3)=max{0.1×0.3×0.3,0.16×0.2×0.3,0.28×0.5×0.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\}P(s2,3​)=max{0.1×0.3×0.3,0.16×0.2×0.3,0.28×0.5×0.3}$=0.042$

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

(3)在t=3时，对每一个状态，求观测为0的
 最大概率 $$P(s_{3,j})=\underset{1\le k\le 3}{\max }\left[P(s_{2,k})a_{kj}\right]b_j(0)$$
 当前最优的前一个状态$$Pre(s_{3,j})=arg\underset{1\le k\le 3}{\max }\left[P(s_{2,k})a_{kj}\right]$$，$j=1,2,3.$

P(s3,1)=max{0.028×0.5×0.5,0.0504×0.3×0.5,0.042×0.2×0.5}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\}P(s3,1​)=max{0.028×0.5×0.5,0.0504×0.3×0.5,0.042×0.2×0.5}$=0.00756$

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

P(s3,2)=max{0.028×0.2×0.4,0.0504×0.5×0.4,0.042×0.3×0.4}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\}P(s3,2​)=max{0.028×0.2×0.4,0.0504×0.5×0.4,0.042×0.3×0.4}$=0.01008$

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

P(s3,3)=max{0.028×0.3×0.7,0.0504×0.2×0.7,0.042×0.5×0.7}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\}P(s3,3​)=max{0.028×0.3×0.7,0.0504×0.2×0.7,0.042×0.5×0.7}$=0.0147$

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

(4)得到结果.
 最大概率
  $P^{\ast }={\max }_{1\le j\le 3}P\left(s_{3,j}\right)$
   =max{0.00756,0.01008,0.0147}=max\left\{ 0.00756,0.01008,0.0147 \right\}=max{0.00756,0.01008,0.0147}
   $=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