HMM算法

14. HMM

1. HMM Background

机器学习大致分为两个派别，也就是频率派和贝叶斯派。

1. 频率派，称为一门独立学科叫统计机器学习，核心问题是优化问题。主要有 3 个步骤：

• model，如 $f(w) = w^T + b$ 等。
• strategy，loss function。
• algorithm：GD/SGD/ 牛顿法 / 拟牛顿法。

1. 贝叶斯派，衍生出概率图模型 PGM，最终是推断，即求解后验分布 $$p(z|x)$$ ，因为贝叶斯公式中，分母需要对所有状态空间进行积分，非常复杂，必须要解决这个数值积分问题，才能让贝叶斯网络能实际运用。这类方法代表有 MCMC，MCMC 发展出来的一系列数值积分方法把贝叶斯理论引入实际应用中。

Hidden Markov Model 隐马尔可夫模型 在概率图模型中的背景：

离散动态模型中 Kalman filter 和 Particle filter 将会在后续章节介绍。

2. HMM 算法整体介绍

1 个基本模型

隐马尔可夫模型是关于 时序 的概率模型 ，描述一个 隐藏的马尔可夫链 随机 生成不可观测的状态随机序列 ，再由各个 状态 生成一个观测从而产生 观测随机序列的过程。HMM 结构图如下：

我们引入些符号来定义模型：

• 状态变量定义为：$I ={i_1, i_2, \cdots, i_t}$，如图中 I 这个序列。
• 观测变量定义为：$O={o_1, o_2, \cdots, o_t}$, 如图中 O 这个序列。
• 状态变量值的集合用 $Q = {q_1, q_2, \cdots, q_N}$, 就是说状态变量 $i_t$ 的值域，每个观测变量 $i_t$ 有 N 个状态可以取。
• 观测变量值的集合用 $V = {v_1, v_2, \cdots, v_M}$, 同样 V 是观测变量 $o_t$ 的值域。
• 状态转移矩阵 $A = [a_{ij}]{N\times N}$, $a{ij} = P(i_{t+1}=q_j| i_t=q_i)$
  这里表示的就是状态变量从前一时刻 t 到下一时刻 t + 1 的概率值，就是 I 表示的空心圈状态序列从前往后变化的概率，根据上面定义的符号这是个 $N \times N$ 的矩阵。
• 发射矩阵(观测概率矩阵)$B =[b_j(k)]_{N \times M}, \quad b_j(k) = P(o_t =v_k|i_t=q_j)$, 这里表示的是在 t 时刻，从状态变量 $i_t=q_j$ 转移到观测变量 $o_t=v_k$ 的概率值，它是 $N \times M$ 的。
• 还有个初始状态概率 $\pi = {\pi_1, \pi_2, \cdots, \pi_N}$. 在 POS 中定义为每个词性出现在句首的概率。

这样，我们可以将 HMM 表示为一个三元组 $\lambda = (\pi, A, B)$。

2 个基本假设

另外，我们对于 HMM 有两个基本假设 李航——统计学习方法 P194：

1. 齐次马尔可夫性假设 ，就是假设隐藏状态马尔可夫链在任意时刻只与前一时刻状态有关。
   $$P(i_t|i_{t-1}, o_{t-1}, \cdots, i_1, o_1) = P(i_t|i_{t-1}) \tag{1}$$
   就是满足该假设情况下，我们不用考虑隐状态和观察状态关系，也不需要考虑隐状态之前所有状态关系，只用计算前一时刻。
2. 观测独立性假设：即任意时刻的观测只依赖该时刻马尔科夫链的状态，与其他观测和其他状态无关。

$$P(o_t|i_t, o_{t-1}, i_{t-1}, \cdots, o_2, i_2, o_1, i_1) = P(o_t|i_t) \tag{2}$$

从图上看，就是计算 t 时刻的观测状态概率只需要考虑该时刻的隐状态，不需要其它状态和节点。

3 个基本问题

1. 概率计算问题：上面部分我们给定的 HMM $\lambda = (\pi, A, B)$ 和观测序列 $O$, 计算其在模型 $\lambda$ 下，序列 $O$ 出现的概率 $P(O|\lambda)$。具体来说，就是前向后向算法（接下来会具体介绍）。
2. 学习问题: 已知观测序列 $O$，估计模型 $\lambda = (\pi, A, B)$ 参数，使得在该模型下观测序列概率 $P(O|\lambda)$ 最大。如 POS 中就是已知：

• 在有词性的初始状态(句首某个词性的概率)$\pi$，
• 当前词性转移到下一个词性的转移矩阵 $A$，
• 以及从该时刻的词性发射出具体某个词的发射矩阵 $B$ 的情况下，有具体的输出序列 $O$，因为已经生成了序列 $O$, 那么有观测序列概率 $P(O|\lambda)$ 最大。
• 具体求解 $\lambda$ 算法有：EM 算法和 Baum Welch 算法。

1. Decoding 问题：求解后的模型 $\lambda = (\pi, A, B)$ 和观测序列 $O={o_1, o_2, \cdots, o_t}$, 求给定观测序列条件概率 $P(I|O)$ 最大的状态序列 $I = {i_1, i_2, \cdots, i_t}$。简单来说，就是在给定模型和观测序列后，求解隐状态序列。这分为预测和滤波：

• 预测问题：$P(i_{t+1}|o_1, \cdots, o_t)$, 也就是在已知当前观测变量情况下预测下一个状态，具体就是 Viterbi 算法。
• 滤波问题：$P(i_t|o_1, \cdots, o_t)$, 已知观测变量序列求 t 时刻隐变量。

接下来，我们逐个介绍 3 个问题的具体求解过程。

3. HMM 概率计算问题

由上节，我们知道概率计算问题就是在给定模型参数 $\lambda$ 情况下，求解 $P(O|\lambda)$，总的来说一下几种计算方法：

1. 直接计算法 我们利用边缘概率和联合概率乘法公式有：

$$P(O|\lambda) = \sum_I P(O, I | \lambda) = \sum_I P(O |I, \lambda)P(I |\lambda) \tag{3}$$
其中 $\sum_I$ 表示所有隐状态序列求和，$P(O|I, \lambda)$ 表示对应隐状态产生观测状态序列的概率，$P(I|\lambda)$ 是某个隐状态出现概率。

并且利用链式法则有：
$$\begin{aligned} P(I|\lambda)& = P(i_1, i_2\cdots , i_{T}|\lambda) \\ &=P(i_T|i_1, i_2\cdots , i_{T-1},\lambda) \cdots P(i_2|i_1, \lambda)P(i_1|\lambda) \end{aligned} \tag{4}$$
而根据齐次马尔可夫假设式 1 有：
$$P(i_T|i_1, i_2\cdots , i_{T-1},\lambda) = P(i_T|i_{T-1})=a_{i_{T-1}i_T }\tag{5}$$
那么式 4 根据式 5 依次代入替换有：
$$\begin{align} P(I|\lambda) &=P(i_T|i_1, i_2\cdots , i_{T-1},\lambda) \cdots P(i_2|i_1, \lambda)P(i_1|\lambda) \\ &=a_{i_{T-1}i_T }a_{i_{T-2}i_{T-1}}\cdots a_{i_{1}i_{2}} \pi_{i_1}\\ &=\pi_{i_1}\prod_{t=2}^T a_{i_{T-1}i_T } \end{align} \tag{6}$$
对于公式 3 我们只要求计算出 $$P(O|I, \lambda)$$ 就可以了，利用观测独立性假设及联合概率有：
$$\begin{align} P(O|I, \lambda) &= P(o_1, o_2, \cdots, o_T|I, \lambda)\\ &= \prod_{i=1}^TP(o_t|I, \lambda)\\ &= \prod_{i=1}^T b_{i_t}(o_t) \end{align} \tag{7}$$
现在公式 3 就等于：
$$\begin{align} P(O|\lambda) &= \sum_I P(O |I, \lambda)P(I |\lambda)\\ & = \sum_I \pi_{i_1}\prod_{t=2}^T a_{i_{T-1}i_T } \prod_{i=1}^T b_{i_t}(o_t) \\ & = \sum_{i_1} \sum_{i_2} \cdots \sum_{i_T} \pi_{i_1}\prod_{t=2}^T a_{i_{T-1}i_T } \prod_{i=1}^T b_{i_t}(o_t) \end{align} \tag{8}$$
公式 8 最后一步只是为了说明这将会有 $O(N^T)$ 个求和，总共计算复杂度为 $O(TN^T)$。在实际中是不可行的，为此我们引入前向 - 后向算法（forward-backward algorithm）。

2. 前向 - 后向算法

前向算法

如上图橙色阴影所示，到时刻 t，隐状态变量为 $i_t=q_i$, 观测状态序列为 $o_1, o_2, \cdots, o_t$ 的概率，为前向概率，我们引入记号 $\alpha_t(i) = P(o_1, o_2, \cdots, o_t，i_t=q_i|\lambda)$ 来表示，更普遍一点表示为 $\alpha_T(i) = P(O,i_T= q_{T}|\lambda)$。

那么 $P(O|\lambda) = \sum_{i=1}^NP(O, i_t=q_i|\lambda)$ , 这表示时刻 t 时隐状态为 $q_i$, 观测状态序列为 $O$ 的前向概率。现在我们利用上述记号来推导:

$$\begin{align} \alpha_{t+1}(j) &= P(o_1, \cdots, o_t, o_{t+1}, i_{t+1}=q_j|\lambda)\\ &=\sum_{i=1}^N P(o_1, \cdots, o_t, o_{t+1}, i_{t+1}=q_j, i_t=q_i|\lambda)\\ &=\sum_{i=1}^N P(o_{t+1}|o_1, \cdots, o_t,i_{t+1}=q_j, i_t=q_i, \lambda) \cdot P(o_1, \cdots, o_t,i_{t+1}=q_j, i_t=q_i, \lambda) \\ &= \sum_{i=1}^N P(o_{t+1}|i_{t+1}=q_j) \cdot P(o_1, \cdots, o_t,i_{t+1}=q_j, i_t=q_i, \lambda) \quad \quad\text{式 2：观测独立性假设}\\ & =\sum_{i=1}^N P(o_{t+1}|i_{t+1}=q_j)\cdot P(i_{t+1}=q_j|o_1, \cdots, o_t, i_t=q_i, \lambda)\cdot P(o_1, \cdots, o_t, i_t=q_i| \lambda)\\ & =\sum_{i=1}^N b_j(o_{t+1}) a_{ij} \alpha_t(i) \end{align} \tag{9}$$

这表示，t 时刻隐状态为 $q_i$, 观测状态为 $O$ 乘以从隐状态 $i$ 转移到 $j$ 的概率再乘以由隐状态变量 $i_{t+1}=q_j$ 观测到 $o_{t+1}$ 的概率再求和。

下面借用 14 隐马尔可夫模型.pdf笔记中图来加强理解，对于前向算法来说，每个时刻 t 的输出状态都等于上一个隐状态分别计算概率再求和 。这实际上跟神经网络非常类似，就是前一状态值乘以权重。因为隐状态可取值为 $N$ 个，乘以对应概率，在对总共的序列长度 $T$ 求和，因此总的时间复杂度为 $O(TN^2)$，减少了大量计算。还可以看看 统计机器学习 P200 例 10.2。

后向算法

如上图所示，我们记 $\beta_t(i) = P(o_{t+1}, \cdots, o_T|i_t=q_i, \lambda)$, 就是 t 时刻在隐状态 $i_t=q_i$ 条件下，从 $t+1$ 到 $T$ 部分观测序列为 $o_{t+1}, \cdots, o_T$ 的概率，这就是后向概率。那么 $\beta_t(1) = P(o_{2}, \cdots, o_T|i_1=q_i, \lambda)$。若要计算 $P(O|\lambda)$，那么有：

$$\begin{aligned} P(O|\lambda) &= P(o_1, o_2, \cdots, o_N|\lambda)\\&= \sum_{i=1}^NP(o_1, o_2, \cdots, o_N, i_1=q_i|\lambda)\quad \quad\text{边缘概率}\\ &= \sum_{i=1}^NP(o_1| o_2, \cdots, o_N, i_1=q_i,\lambda){\underbrace{P(i_1=q_i|\lambda)}{\pi_i} }\quad \text{条件概率}\\ &= \sum{i=1}^NP(o_1| o_2, \cdots, o_N, i_1=q_i,\lambda)\cdot P(o_2, \cdots, o_N, i_1=q_i|\lambda) \cdot \pi_i\\&=\sum_{i=1}^NP(o_1|i_1=q_i, \lambda) \cdot\beta_1(i)\cdot \pi_i\ \\&=\sum_{i=1}^Nb_i(o_1)\cdot\beta_1(i)\cdot \pi_i\\ &=\sum_{i=1}^N\pi_i\cdot b_i(o_1)\cdot \beta_1(i) \end{aligned} \tag{10}$$
其中，$\pi_i$ 为某个初始状态的概率，$b_i(o_1)$ 第一个隐状态条件下生成第一个观测状态概率，$\beta_1(i)$ 第 1 个观测状态确定后，生成后面观测序列的概率。另外，对于上述推导过程中，可以将 $\lambda$ 省去再思考推导过程。依旧借用 14 隐马尔可夫模型.pdf 中图如下，要确定输出观测序列，我们要知道初始隐状态概率，该隐状态转化为第一个观测状态概率，以及第一个观测条件下，生成后续观测序列的可能性只要将其相乘就可以了，分别对应下图蓝色，绿色，红色部分，并且可以看作一个信息传递过程。

现在，我们要推广到 $\beta_t(i)$ 情况，这里省去 $\lambda$，方便推导表达式。
$$\begin{align} \beta_t(i) &= P(O|\lambda) \\&= P(o_{t+1}, o_{t+2}, \cdots, o_T|i_t=q_i)\\ &= \sum_{j=1}^N P(o_{t+1}, o_{t+2}, \cdots, o_T| i_{t+1}=q_j,i_t=q_i)P(o_{t+1}, o_{t+2}, \cdots, o_T, i_{t+1}=q_j|i_t=q_i) \quad{链式法则}\\ &= \sum_{j=1}^N P(o_{t+1}, o_{t+2}, \cdots, o_T| i_{t+1}=q_j,i_t=q_i)\cdot P(i_{t+1}=q_j|i_t=q_i) \quad{齐次假设} \\ &= \sum_{j=1}^N P(o_{t+1}, o_{t+2}, \cdots, o_T| i_{t+1}=q_j)\cdot a_{ij} \\ &= \sum_{j=1}^N P(o_{t+1}| o_{t+2}, \cdots, o_T, i_{t+1}=q_j)\cdot P(o_{t+2}, \cdots, o_T|i_{t+1}=q_j) \cdot a_{ij} \\ &= \sum_{j=1}^N b_{j}(o_{t+1}) \cdot \beta_{t+1}(j) \cdot a_{ij}\\ &= \sum_{j=1}^N a_{ij}\cdot b_{j}(o_{t+1}) \cdot \beta_{t+1}(j) \end{align}\tag{11}$$
跟上面一样的递推过程，依旧如下图的后向算法计算过程。

例：考虑盒子和球模型 $\lambda = (\pi,A,B)$，状态集合 $Q={1,2, 3}$，观测集合 $V={红, 白}$。(李航《统计学方法》P200)

$$A = \begin{bmatrix}0.5 & 0.2 & 0.3\ 0.3& 0.5 & 0.2\ 0.2 & 0.3 & 0.5 \end{bmatrix}，B = \begin{bmatrix}0.5 & 0.5\ 0.4& 0.6\ 0.7 & 0.3 \end{bmatrix}，\pi = (0.2,0.4,0.4)^{T}$$

设观测序列长度 $T=3$，观测序列为：$O=(红, 白, 红)$，求：$P(O|\lambda)$。代码来自 HMM-Evaluation.

import numpy as np
def forward(lamda, O):
    '''
HMM 模型：概率计算问题——前向算法，李航《统计学方法》P177
step1:
    t= 0 时刻，计算初值：alpha[t=1][j] = pai(j) * B[j][o1]
step2:
    递推：t=1, 2, ..., T-1，计算 alpha[t][j]：alpha[t][j] = ∑(alpha[t-1][i] * A[i][j]) * B[j][ot]
step3:
    计算：P(O|λ) = ∑alpha[T][j]
'''
    A, B, pai = lamda  # HMM 的参数：A, B, pai
    N, M = B.shape  # 状态集合 Q 的长度 N，观测集合 V 的长度 M
    T = len(O)  # 由观测序列得到 T 个时刻
    '''HMM 的前向概率矩阵：T*N'''
    alphas = np.zeros((T, N))
    '''step1：t= 0 时刻，计算初值：alpha[t=0][j] = pai(j) * B[j][o1]'''
    for j in range(N):  # j 是当前状态
        k = V.index(O[0])  # 第 Ot 个观测对应的 V 观测集合的位置 k:O[0]=V[k]
        alphas[0][j] = pai[j] * B[j][k]
    '''step2：递推：t=1, 2, ..., T-1，计算：alpha[t][j] = ∑(alpha[t-1][i] * A[i][j]) * B[j][Ot]'''
    for t in range(1, T):  # j 是当前状态
        k = V.index(O[t])  # 第 Ot 个观测对应的 V 观测集合的位置 k:O[t]=V[k]
        for j in range(N):
            temp = 0
            for i in range(N):
                temp += alphas[t - 1][i] * A[i][j]  # 计算：∑(alpha[t-1][i] * A[i][j])
            alphas[t][j] = temp * B[j][k]  # 计算：∑(alpha[t-1][i] * A[i][j]) * B[j][Ot]
    '''step3：计算：P(O|λ) = ∑alpha[T][j]'''
    prob = 0  # 求解 P(O|λ)
    for j in range(N):  # j 是当前隐状态
        prob += alphas[T - 1][j]
    return prob
if __name__ == '__main__':
    '''状态集合：Q = {q1, q2, ..., qN}'''
    Q = ['盒子 1', '盒子 2', '盒子 3']
    '''观测集合：V = {v1, v2, ..., vM}'''
    V = ['红', '白']
    '''
        状态转移概率矩阵：A = [aij]，其中 aij = P(it+1=qj | it=qi)
        A: shape N * N (盒子之间互相转换)
    '''
    A = np.array([[0.5, 0.2, 0.3],
        [0.3, 0.5, 0.2],
        [0.2, 0.3, 0.5]])
    '''
        观测概率矩阵：B=[bj(k)]，其中 bj(k) = P(ot=vk | it=qj)
        B: shape N * M (N 个盒子 M 个观测)
    '''
    B = np.array([[0.5, 0.5],
        [0.4, 0.6],
        [0.7, 0.3]])
    '''初始状态概率向量：pai = P(i1=qi)：第一次出现每个盒子的概率分布'''
    pai = np.array([0.2, 0.4, 0.4])
    '''随机序列 1——观测序列：O = {o1, o2, ..., oT}'''
    '''随机序列 2——状态序列（隐序列，叠加态的黑盒子）：I = {i1, i2, ..., iT}'''
    O = np.array(['红', '白', '红'])
    '''前向算法计算：P(O|λ)'''
    P_O_or_lamda = forward(lamda=[A, B, pai], O=O)
    print('通过 HMM 前向算法计算概率：\n 观测序列：{} 出现的概率 P(O|λ) = {:.4f}'.format(O, P_O_or_lamda))

4. Viterbi 算法

隐马尔可夫模型中的维特比算法，用于解码给定的观测序列和模型参数，得到最有可能的隐藏状态序列。

该代码实现了 HMM 模型的三个基本问题：

1. 概率计算问题：给定模型参数和观测序列，计算这个序列出现的概率。
2. 学习问题：给定观测序列，估计模型参数，使得这个观测序列出现的概率最大。
3. 解码问题：给定模型参数和观测序列，找到最有可能的隐藏状态序列。

具体来说，该代码实现了维特比算法（Viterbi algorithm），用于解码问题，即在给定模型参数和观测序列的情况下，找到最有可能的隐藏状态序列。

代码中的 viterbi 函数接受两个参数：一个观测序列 obs，一个 HMM 模型 model。model 包括三个参数：初始状态概率向量 model['start_prob']、状态转移概率矩阵 model['trans_prob']、观测概率矩阵 model['emit_prob']。

函数中的主要步骤是：

1. 初始化：将起始概率向量 model['start_prob'] 和第一个观测值对应的观测概率向量 model['emit_prob'][:, obs[0]] 相乘，得到一个长度为状态数的向量 probs，并将其作为维特比算法的第一个状态。
2. 递推：对于观测序列中的每个观测值，根据当前状态向量 probs 和状态转移概率矩阵 model['trans_prob']，计算出下一个状态向量 next_probs，并将其作为下一次迭代的当前状态向量。
3. 回溯：在所有观测值都处理完之后，根据每个状态向量的最大值以及指向最大值的箭头，得到最有可能的隐藏状态序列。

import numpy as np
 
# Define the Viterbi algorithm
def viterbi(obs, states, obs_symbols, initial_prob, transitions, emissions):
    """
    viterbi algorithm to get best_path and prob_best_path
    :param obs: 要已经划分完毕的标注的句子, split/tokenized
    :param states: 词性的所有状态
    :param obs_symbols: 所有词的列表
    :param initial_prob: 初始概率矩阵
    :param transitions: 转移矩阵
    :param emissions: 发射矩阵
    :return:
    """
    # Define the T x N matrix for probabilities
    T = len(obs) #要标注的句子长度
    N = len(states) #pos tag 长度
    probs = np.zeros((T, N)) #保存句子每个词到指定 tag 的概率
 
    # Define the T x N matrix for pointers to backtrack
    pointers = np.zeros((T, N), dtype=int)
 
    # Calculate probabilities and pointers for the first observation
    # 初始化第一个词对应词性 prob
    obs_index = obs_symbols.index(obs[0])
    probs[0] = np.multiply(initial_prob, emissions[:, obs_index])
 
    # Iterate over the remaining observations
    for t in range(1, T):
        #输入句子的词对应在观察序列列表的索引
        obs_index = obs_symbols.index(obs[t])
        for s in range(N):
            #第 t 个词前一个词的概率乘以转移概率. 这里从 1 循环到 T 句子长度的过程中都是拿前面一个
            #来获得当前词性的概率, 已经初始化得到 0，那么 1 就是 probs[0] *transitions[:, s]
            #迭代获取所有当前词性的概率
            temp_probs = np.multiply(probs[t - 1], transitions[:, s])
            #取最大概率的 index
            max_prob_index = np.argmax(temp_probs)
            #那么对应最大概率 = 发射矩阵矩阵中词 * 当前词性最大概率
            probs[t, s] = emissions[s, obs_index] * temp_probs[max_prob_index]
            #为解码取路径保存索引来 这个索引表示 t 时刻对应上一时刻的最大的 prob 连接的索引
            pointers[t, s] = max_prob_index
 
 
    # Find the best path
    #最后一行的最大的索引
    best_path = [np.argmax(probs[-1])]
    for t in range(T - 1, 0, -1):
        best_path.append(pointers[t, best_path[-1]])
    best_path.reverse()
 
    # Calculate the probability of the best path
    prob_best_path = np.max(probs[-1])
 
    return best_path, prob_best_path
 
# Example usage
obs = ['walk', 'shop', 'clean']
 
# Define states
states = ('Rainy', 'Sunny')
 
# Define observation symbols
obs_symbols = ('walk', 'shop', 'clean')
 
# Define initial probabilities
initial_prob = np.array([0.6, 0.4])
 
# Define transition matrix
transitions = np.array([[0.7, 0.3],
                        [0.4, 0.6]])
 
# Define emission matrix
emissions = np.array([[0.1, 0.4, 0.5],
                      [0.6, 0.3, 0.1]])
 
best_path, prob_best_path = viterbi(obs, states, obs_symbols, initial_prob, transitions, emissions)
print("Best path:", [states[i] for i in best_path])
print("Probability of the best path:", prob_best_path)

Inference

[1] HMM 前向与后向算法

[2] 基于 HMM 模型的词性标注器 Python 实现

[3] HMM-Evaluation