隐马尔科夫模型（HMM）

前言

同时在LP 之分词技术概述中有提到 HMM 模型，虽然此方法在现代的作用和地位有所下降，但是依然是非常值得了解的学习机器学习经典算法。

本文的基本提纲如下：

基础知识
1. 定义
2. 适用场景
3. 隐马尔科夫模型的两个基本假设
4. 实例
HMM模型的三个基本问题
1. 前向后向算法评估观察序列概率
2. 鲍姆-韦尔奇算法求解HMM参数
3. 维特比算法解码隐藏状态序列
分词实战

基础知识

定义

直观定义
首先看一个图：

该图的文字描述是其分成两排，第一排是$y$序列，第二排是$x$序列。每个$x$都只有一个$y$指向它，每个$y$也都有另一个$y$指向它。

下面按照上述描述给出一个定义：

状态序列（上图中的$y$，下面的$I$）：隐藏的马尔科夫链随机生成的状态序列，称为状态序列（state sequence）
观测序列（上图中的$x$，下面的$O$）: 每个状态生成一个观测，而由此产生的观测的随机序列，称为观测序列（obeservation sequence）
隐马尔科夫模型：隐马尔科夫模型是关于时序的概率模型，描述由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程。

序列的每一个位置可以看做是一个时刻，HMM 由初始概率分布、状态转移概率分布和观测概率分布确定。

形式化的定义如下：

形式化定义

设$Q$是所有可能的状态的集合，$V$是所有可能的观测的集合，如下:

$Q={q_1,q_2,...,q_N},V={v_1,v_2,...,v_M}$

其中，$N$是可能的状态数，$M$是可能的观测数。

$I$是长度为$T$的状态序列，$O$是对应的观测序列。

$I=(i_1,i_2,...,i_T),O=(o_1,o_2,...,o_T)$

A是状态转移概率矩阵：

$A=[a_{ij}]_{N×N}$

在时刻$t$，处于$q_i$ 状态的条件下在时刻$t+1$转移到状态$q_j$ 的概率【齐次马尔科夫链假设,下面的第一个假设】：

$a_{ij}=P(i_{t+1}=q_j|i_t=q_i)$

B是观测概率矩阵:

$B=[b_j(k)]_{N×M}$

其中$k=1,2,…,M; j=1,2,…,N$
其中，在时刻$t$处于状态$q_j$ 的条件下生成观测$v_k$ 的概率【观测独立性假设，下面第二个假设】：

$b_j(k)=P(o_t=v_k|i_t=q_j)$

π是初始状态概率向量:

$π=(π_i)$

其中，$π_i=P(i_1=q_i)$
隐马尔科夫模型由初始状态概率向量$π$、状态转移概率矩阵A和观测概率矩阵$B$决定。$π$和$A$决定状态序列，$B$决定观测序列。因此，隐马尔科夫模型$λ$可以由三元符号表示，即：$λ=(A,B,π)$。$A,B,π$称为隐马尔科夫模型的三要素。

适用场景

有了上述的定义可以看出使用HMM模型时我们的问题一般有这两个特征：

我们的问题是基于序列的，比如时间序列，或者状态序列。
我们的问题中有两类数据，一类序列数据是可以观测到的，即观测序列；而另一类数据是不能观察到的，即隐藏状态序列，简称状态序列。

比如语言转文字的问题，发出的一串连续的声音就是观测序列，而实际要表达的一段话的文字就是状态序列，机器的任务，就是从这一串连续的声音中判断出最可能要表达的文字的内容。

隐马尔科夫模型的两个基本假设

设隐马尔科夫链在任意时刻$t$的状态只依赖于其前一时刻的状态，与其他时刻的状态及观测无关，也与时刻$t$无关。（齐次马尔科夫性假设）
假设任意时刻的观测只依赖于该时刻的马尔科夫链的状态，与其他观测和状态无关。（观测独立性假设）

一个例子

李航的《统计学习方法》的例子，假设我们有3个盒子，每个盒子里都有红色和白色两种球，这三个盒子里球的数量分别是：

盒子	1	2	3
红球数	5	4	7
白球数	5	6	3

按照下面的方法从盒子里抽球:

开始的时候，从第一个盒子抽球的概率是0.2，从第二个盒子抽球的概率是0.4，从第三个盒子抽球的概率是0.4。以这个概率抽一次球后，将球放回。
从当前盒子转移到下一个盒子进行抽球。规则是：
1. 如果当前抽球的盒子是第一个盒子，则以0.5的概率仍然留在第一个盒子继续抽球，以0.2的概率去第二个盒子抽球，以0.3的概率去第三个盒子抽球。
2. 如果当前抽球的盒子是第二个盒子，则以0.5的概率仍然留在第二个盒子继续抽球，以0.3的概率去第一个盒子抽球，以0.2的概率去第三个盒子抽球。
3. 如果当前抽球的盒子是第三个盒子，则以0.5的概率仍然留在第三个盒子继续抽球，以0.2的概率去第一个盒子抽球，以0.3的概率去第二个盒子抽球。
4. 如此下去，直到重复三次.

得到一个球的颜色的观测序列:

HMM观测序列的生成
从上面的例子，我们也可以抽象出HMM观测序列生成的过程。
输入的模型是$\lambda = (A, B, \Pi)$，观测序列长度是$T$.
输出是观测序列：$O =\{o_1,o_2,…o_T\}$

以下代码是 HMM 一个模拟，

HMM模型的三个基本问题

HMM在实际应用中，一般会遇上三种问题：

概率计算问题：给定模型$\lambda = (A,B,\pi)$ 和观测序列$O={o_1,o_2,…,o_T}$，计算在模型$\lambda$下观测序列$O$出现的概率$P(O|\lambda)$
学习问题：已知观测序列$O={o_1,o_2,\cdots,o_T}$，估计模型$\lambda= (A,B,\pi)$，使得$P(O|\lambda)$最大。即是用最大似然法的方法估计参数。
预测问题：（亦称为 decoding 解码问题）已知观测序列$O={o_1,o_2,…,o_T}$和模型$\lambda = (A,B,\pi)$，求给定的观测序列条件概率$P(I|O)$最大的状态序列$I = (i_1,i_2,\cdots,i_T)$，即是给定观测序列，求最有可能的状态序列。

回到上面的例子，这三个问题就是：

概率计算问题：如果给定模型参数，抽球结果某一系列观测的序列状态出现的概率
学习问题：根据抽球结果某一些列观测的状态，学习模型参数。
预测问题：根据学到的模型，预测抽球的序列。

前向后向算法评估观察序列概率

首先要解决的问题是上面的问题一，即是概率计算问题，在模型$\lambda$下观测序列$O$出现的概率$P(O|\lambda)$。

直接计算法

因为我们知道所有的隐藏状态之间的转移概率和所有从隐藏状态到观测状态生成概率，那么我们是可以暴力求解的。
我们可以列举出所有可能出现的长度为T的隐藏序列$I = (i_1,i_2,\cdots,i_T)$，分布求出这些隐藏序列与观测序列$O={o_1,o_2,…,o_T}$的联合概率分布$P(O,I|\lambda)$，这样我们就可以求出边缘分布$P(O|\lambda)$

算法步骤如下：

任意一个隐藏序列$I = (i_1,i_2,\cdots,i_T)$出现的概率是: $P(I|\lambda) = \pi_{i_1} a_{i_1i_2} a_{i_2i_3}... a_{i_{T-1}\;\;i_T}$
对上面这种状态序列，产生观测序列$O={o_1,o_2,…,o_T}$的概率是： $P(O|I,\lambda)=b_{i_1}(o_1)b_{i_2}(o_2)...b_{i_T}(o_{T})$
则O和I联合出现的概率是： $P(O,I|\lambda)=P(O|I,\lambda)P(I|\lambda)=\pi_{i_1}b_{i_1}(o_1)a_{i_1i_2}b_{i_2}(o_2)...a_{i_{T-1}i_T}b_{i_T}(o_{T})$
对所有可能的$I$求和，得到: $P(O|λ)=\sum_IP(O,I|\lambda)=\sum_{i_1,...,i_T}\pi_{i_1}b_{i_1}(o_1)a_{i_1i_2}b_{i_2}(o_2)...a_{i_{T-1}i_T}b_{i_T}(o_{T})$ 如果直接计算，时间复杂度太高，是$O(TN^T)$

因此对于一些隐藏状态数极少的模型，我们可以用暴力求解法来得到观测序列出现的概率，但是如果隐藏状态多，则上述算法太耗时，我们需要寻找其他简洁的算法。

前向后向算法就是来帮助我们在较低的时间复杂度情况下求解这个问题的。

前向算法（或者后向算法）

前向算法本质上属于动态规划的算法，也就是我们要通过找到局部状态递推的公式，这样一步步的从子问题的最优解拓展到整个问题的最优解。

在前向算法中，通过定义“前向概率”来定义动态规划的这个局部状态。

前向概率：定义时刻$t$时隐藏状态为$q_i$, 观测状态的序列为$o_1,o_2,\cdots,o_t$的概率为前向概率，记做：

$\alpha_t(i) = P(o_1,o_2,...o_t, i_t =q_i | \lambda)$

利用的思想是动态规划中的递推：

即是利用在时刻$t$掌握的前向概率信息推出$t+1$时刻各个隐藏状态的前向概率。具体步骤如下：

算法步骤如下：

输入：隐马模型$\lambda = (A, B, \Pi)$，观测序列$O=(o_1,o_2,…o_T)$;
输出：观测序列概率$P(O|λ)$.

step 1： 计算$t=1$时刻的各个隐藏状态初始前向概率：

$\alpha_1(i) = \pi_ib_i(o_1),\; i=1,2,...N$

step 2：递推时刻$2,3,\cdots,T$时刻的前向概率：

$\alpha_{t+1}(i) = \Big[\sum\limits_{j=1}^N\alpha_t(j)a_{ji}\Big]b_i(o_{t+1}),\; i=1,2,...N$

step 3：计算最终结果：

$P(O|\lambda) = \sum\limits_{i=1}^N\alpha_T(i)$

从递推公式可以看出，我们的算法时间复杂度是$O(TN^2)$

后向算法原理亦是如此。

实现代码如下:

# 接上面代码
def forward(obs_seq):
    """前向算法"""
    N = A.shape[0]
    T = len(obs_seq)
    
    # F保存前向概率矩阵
    F = np.zeros((N,T))
    F[:,0] = pi * B[:, obs_seq[0]]

    for t in range(1, T):
        for n in range(N):
            F[n,t] = np.dot(F[:,t-1], (A[:,n])) * B[n, obs_seq[t]]

    return F, sum(F[:,T-1])

def backward(obs_seq):
    """后向算法"""
    N = A.shape[0]
    T = len(obs_seq)
    # X保存后向概率矩阵
    X = np.zeros((N,T))
    X[:,-1:] = 1

    for t in reversed(range(T-1)):
        for n in range(N):
            X[n,t] = np.sum(X[:,t+1] * A[n,:] * B[:, obs_seq[t+1]])

    return X, sum(pi*B[:,0]*X[:,0])

结果如下：

鲍姆-韦尔奇算法求解HMM参数

HMM模型参数求解根据已知的条件可以分为两种情况。
第一种情况较为简单，就是我们已知$D$个长度为$T$的观测序列和对应的隐藏状态序列，$\{(O_1, I_1), (O_2, I_2), …(O_D, I_D)\}$是已知的，此时我们可以很容易的用最大似然来求解模型参数。

可见第一种情况下求解模型还是很简单的。但是在很多时候，我们无法得到HMM样本观察序列对应的隐藏序列，只有D个长度为T的观测序列，即是$\{(O_1), (O_2), …(O_D)\}$是已知的，此时我们能不能求出合适的HMM模型参数呢？这就是我们的第二种情况，它的解法最常用的是鲍姆-韦尔奇算法，其实就是基于EM算法的求解。

前言