机器学习之隐马尔科夫模型(HMM)原理及Python实现 (大章节)

HMM

隐马尔可夫模型（hidden Markov model, HMM）是可用于标注问题的统计学模型，是生成模型。

本章节内容参考李航博士的《统计学习方法》
本章节添加了一些结论性结果的推导过程。

1. 从一个天然语言处理的例子开始

例若有三个个句子：
句子一：我/名词 看见/动词 猫/名词
句子二：猫/名词 是/动词 可爱的/形容词
句子三：我/名词 是/动词 可爱的/形容词
通常只能观察到具体的词，因此像"我 看见 猫 …"是观测集合，而词性如"名词 动词 形容词 …"是状态序列

设 $Q$是全部可能的状态集合， $V$是全部可能的观测集合：

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

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

例如： $Q=\{名词，动词，形容词 \}，V=\{我， 看见， 猫， 是，可爱的\}，N=3, M=5$

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

$I = \{i_1, i_2,..., i_T \}, O=\{o_1, o_2,..., o_T\}$

例如： $I=(名词，动词，名词)， O=(我，看见，猫)$

$A$是状态转移矩阵：

$A=[a_{ij}]_{N*N} \tag1$

其中，

$a_{ij} = p(i_{t+1}=q_j|i_t=q_i), i=1,2,...,N; j=1,2,...,N \tag2$

例如：

转态转移几率	名词	动词	形容词
名词	0	1	0
动词	1/3	0	2/3
形容词	1/3	1/3	1/3

$B$是观测几率矩阵，也就是发射矩阵：

$B=[b_j(k)]_{N*M} \tag3$

其中，

$b_j(k) = p(o_t=v_k|i_t=q_j), k=1,2,...,M; j=1,2,...,N \tag4$

例如：

观测矩阵几率	我	看见	猫	是	可爱的
名词	1	0	1	0	0
动词	0	1	0	1	0
形容词	0	0	0	0	1

$\pi$是初始状态几率向量：

$\pi = (\pi_i) \tag5$

其中，

$\pi_i = p(i_1 = q_i), i = 1,2,...,N \tag6$

$A,B$和 $\pi$是HMM的参数，用 $\lambda$表示：

$\lambda = (A,B,\pi) \tag7$

例如：

名词	动词	形容词
1	0	0

隐马尔可夫的三个基本问题
1.几率计算问题。给定模型 $\lambda=(A,B,\pi)$和观测序列 $O=(o_1,o_2,...,o_T)$,计算在已知模型参数的状况下，观测序列的几率，即 $p(O|\lambda)$。
2.学习问题。已知观测序列 $O=(o_1,o_2,...,o_T)$,估计模型参数 $\lambda=(A,B,\pi)$，使 $p(O|\lambda)$最大。
3.预测问题，也称解码问题。已知模型 $\lambda=(A,B,\pi)$和 $O=(o_1,o_2,...,o_T)$，求条件几率最大 $p(I|O)$最大的状态序列 $I=(i_1,i_2,...,i_T)$。

2. 几率预测问题

几率问题预测用直接计算法，计算复杂度高，能够采用动态规划形式的前向和后向算法下降计算复杂度。
为了表示方便，记：

$(o_{1:t} )= (o_1,o_2,...,o_n); (o_{t_:T})=(o_t,o_{t+1},...,o_T)$

2.1 前向算法

接下来就是解前向几率 $p(i_t,o_{1:t}|\lambda)$：

$\begin{aligned} p(i_t,o_{1:t}|\lambda) &=\sum_{i_{t-1}} p(i_{t-1},i_t,o_{1:t-1},o_t|\lambda) \\ &=\sum_{i_{t-1}} p(o_t|i_{t-1},i_t,o_{1:t-1},\lambda)p(i_t|i_{t-1},o_{1:t-1},\lambda)p(i_{t-1},o_{1:t-1}|\lambda) \end{aligned}$

由隐马尔科夫的条件独立性假设可得：

$p(o_t|i_{t-1},i_t,o_{1:t-1},\lambda) = p(o_t|i_t,\lambda)$

$p(i_t|i_{t-1},o_{1:t-1},\lambda)=p(i_t|i_{t-1},\lambda)$

故

$p(i_t,o_{1:t}|\lambda)=\sum_{i_{t-1}} p(o_t|i_t,\lambda) p(i_t|i_{t-1},\lambda)p(i_{t-1},o_{1:t-1}|\lambda)=[\sum_{i_{t-1} } p(i_{t-1},o_{1:t-1}|\lambda) p(i_t|i_{t-1},\lambda)] p(o_t|i_t,\lambda)$

设：

$\alpha_{t+1}(i) = p(o_{1:t+1},i_{t+1}=q_i|\lambda) \tag8$

且：

$p(i_{t+1}=q_i|i_t=q_j,\lambda)] = a_{ji}$

$p(o_{t+1}|i_{t+1},\lambda)=b_i(o_{t+1})$

则：

$\alpha_{t+1}(i)=[\sum_{j=1}^N \alpha_t(j)a_{ji}]b_i(o_{t+1}) \tag9$

因此前向算法就可迭代进行。

前向算法：
1.初值

$\alpha_1(i) = \pi_ib_i(o_1)$

2.递推 $t=1,2,...,T-1$

$\alpha_{t+1}(i)=[\sum_{j=1}^N \alpha_t(j)a_{ji}]b_i(o_{t+1})$

3.终止
$p(O|\lambda) = \sum_{i=1}^N \alpha_T(i)$

2.2 后向算法

后向算法解决后向几率 $p(o_{t+1:T}|i_t, \lambda)$:

$\begin{aligned} p(o_{t+1:T}|i_t, \lambda) &= \sum_{i_{t+1}} p(i_{t+1},o_{t+1},o_{t+2:T} | i_t, \lambda) \\ &= \sum_{i_{t+1}} p(o_{t+2:T}|i_{t+1}, i_t, o_{t+1}, \lambda) p(o_{t+1}|i_{t+1}, i_t, \lambda) p(i_{t+1}|i_t,\lambda)\\ \end{aligned}$

由隐马尔科夫的条件独立假设得：

$p(o_{t+2:T}|i_{t+1}, i_t, o_{t+1}, \lambda)=p(o_{t+2:T}|i_{t+1}, \lambda)$

$p(o_{t+1}|i_{t+1}, i_t, \lambda) = p(o_{t+1}|i_{t+1}, \lambda)$

设：

$\beta_t(i) = p(o_{t+1:T}|i_t=q_i, \lambda) \tag{10}$

又：

$p(i_{t+1}=q_j|i_t=q_i,\lambda) = a_{ij}$

$p(o_{t+1}|i_{t+1}=q_j, \lambda) = b_j(o_{t+1})$

则：

$\beta_t(i) = \sum_{j=1}^N a_{ij} b_j(o_{t+1}) \beta_{t+1}(i) \tag{11}$

后向算法：
(1)

$\beta_T (i) = 1$

(2) 对t=T-1,T-2,…,1

$\beta_t(i) = \sum_{j=1}^N a_{ij} b_j(o_{t+1}) \beta_{t+1}(i)$

(3)

$p(O|\lambda) = \sum_{i=1}^N \pi_i b_i(o_1) \beta_1(i)$

2.3 一些几率与指望值

这两个指望值都是后面EM算法用到的中间参量
1.计算 $t$时刻处于状态 $q_i$的几率。
几率计算问题是计算 $p(O|\lambda)$，则有：

$p(O|\lambda)=\sum_{i_t}p(O,i_t|\lambda)$

依据隐马尔科夫的独立性假设：

$p(o_{t+1:T}|i_t,o_{1:t}, \lambda) = p(o_{t+1:T}|i_t, \lambda)$

因此：

$\begin{aligned} p(O|\lambda) &=\sum_{i_t}p(O,i_t|\lambda) \\ &=\sum_{i_t} p(o_{t+1:T}|i_t,o_{1:t}, \lambda) p(i_t,o_{1:t}|\lambda) \\ &=\sum_{i_t} p(o_{t+1:T}|i_t, \lambda) p(i_t,o_{1:t}|\lambda) \\ \end{aligned}$

又有：

$\alpha_t(i) = p(o_{1:t},i_t=q_i|\lambda) \tag{12}$

$\beta_t(i) = p(o_{t+1:T}|i_t=q_i, \lambda) \tag{13}$

故：

$p(O,i_t=q_i|\lambda) = p(o_{t+1:T}|i_t=q_i, \lambda) p(i_t=q_i,o_{1:t}|\lambda) = \alpha_t(i) \beta_t(i)$

$p(O|\lambda) = \sum_{i_t} \alpha_t(i) \beta_t(i)$

设：

$\gamma_t(i) = p(i_t=q_i|O,\lambda)$

因而能够获得：

$\gamma_t(i) = p(i_t=q_i|O,\lambda) = \frac {p(i_t=q_i,O|\lambda)}{p(O|\lambda)} = \frac {\alpha_t(i) \beta_t(i)}{\sum_{j=1}^N \alpha_t(j) \beta_t(j)} \tag{14}$

2.计算计算 $t$时刻处于状态 $q_i$且计算 $t+1$时刻处于状态 $q_j$的几率

$\begin{aligned} p(O|\lambda) &=\sum_{i_t} \sum_{i_{t+1}} p(O,i_t, i_{t+1}|\lambda) \\ &=\sum_{i_t} \sum_{i_{t+1}} p(o_{1:t},o_{t+1},o_{t+2:T},i_t, i_{t+1}|\lambda) \\ &=\sum_{i_t} \sum_{i_{t+1}} p(o_{t+2:T}|o_{1:t},o_{t+1},i_t, i_{t+1},\lambda)p(o_{t+1}|o_{1:t},i_t,i_{t+1},\lambda) p(i_{t+1}|i_t,o_{1:t},\lambda) p(i_t,o_{1:t}|\lambda) \\ \end{aligned}$

由隐马尔科夫的独立性假设可得：

$p(O|\lambda) = \sum_{i_t} \sum_{i_{t+1}} p(o_{t+2:T}| i_{t+1},\lambda)p(o_{t+1}|i_{t+1},\lambda) p(i_{t+1}|i_t,\lambda) p(i_t,o_{1:t}|\lambda)$

设：

$\xi_t(i,j)=p(i_t=q_i,i_{t+1}=q_j|O,\lambda)$

又有公式(2)(4)(12)(13)

得：

$\xi_t(i,j) = \frac {p(i_t=q_i,i_{t+1}=q_j|O,\lambda)}{p(O|\lambda)} =\frac {\alpha_t(i) a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)} {\sum_{i=1}^N \sum_{j=1}^N \alpha_t(i) a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)} \tag{15}$

3. 学习问题

3.1 监督学习

若是有标记好状态序列的样本，那就太好办了，直接将接个矩阵统计的各个维度定义后进行统计就能够了。统计过程当中注意几率之和为一的约束。

3.2 无监督学习

若是没有标记状态序列的样本，能够用Baum-Welch算法(EM算法)实现。

已知：包含 $S$个长度为 $T$的观测序列的观测序列 $\{O_1,O_2,...,O_S \}$
目标：学习隐马尔可夫模型的参数 $\lambda=(A,B,\pi)$

记观测数据 $O$,隐数据 $I$，那么隐马尔可夫模型能够表示为：

$p(O|\lambda) = \sum_I p(O|I,\lambda) p(I|\lambda)$

E步：

由于对 $\lambda$而言， $1/p(O| \overline \lambda)$是常数项，因此

$\begin{aligned} Q(\lambda,\overline \lambda) &= E_I[\log p(O,I|\lambda)|O, \overline \lambda] \\ &= \sum_I \log p(O,I|\lambda) p(I|O,\overline \lambda) \\ &= \sum_I \log p(O,I|\lambda) \frac {p(I,O|\overline \lambda)}{p(O| \overline \lambda)} \\ &= \sum_I \log p(O,I|\lambda) p(I,O|\overline \lambda) \\ \end{aligned}$

将几率计算问题2.1小姐中前向算法的递归公式展开就能够获得：

$p(O,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_{iT}(o_T) = \pi_{i_1} [\prod_{t=1}^{T-1} a_{i_ti_{t+1}}][\prod_{t=1}^T b_{it}(o_t)]$

因而：

$Q(\lambda, \overline \lambda) = \sum_I \log \pi_{i_1} p(O, I| \overline \lambda) + \sum_I (\sum_{t=1}^{T-1} a_{i_ti_{t+1}}) p(O, I| \overline \lambda) + \sum_I (\sum_{t=1}^T b_{it}(o_t)) p(O, I| \overline \lambda) \tag{16}$

特此说明隐变量
隐马尔可夫模型的隐变量就是观测序列对应的状态序列，因此隐变量能够用（14）式的变量表示
后面在M步中更新模型参数的时候也用到了（15）式，是否是就说明隐变量是两个，其实不是的，这儿只是为了表示的方便和算法的方便。
也就是在E步中，用 $\gamma$和 $\xi$表示隐变量，只是为了编程和表示的便利，这两个变量在E步中信息是重复的。

M步：

1.求解 $\pi_i$
由(15)式可得：

$L(\pi_{i_1}) = \sum_I \log \pi_{i_1} p(O, I| \overline \lambda) = \sum_{i}^N \log \pi_{i_1} p(O, i_1=i| \overline \lambda)$

又由于 $\pi_i$知足约束条件 $\sum_{i=1}^N \pi_{i_1}=1$，利用拉格朗日乘子法，写出拉格朗日函数：

$\sum_{i=1}^N \log \pi_{i} p(O, i_1=i| \overline \lambda) + \gamma(\sum_{i=1}^N \pi_{i} - 1)$

对其求偏导而且令其结果为0得：

$\frac {\partial} {\partial \pi_i} [\sum_{i=1}^N \log \pi_{i} p(O, i=i| \overline \lambda) + \gamma(\sum_{i_1=1}^N \pi_{i} - 1)]=0 \tag{17}$

得：

$p(O, i_1=i| \overline \lambda) + \gamma \pi_i=0$

获得：

$\pi_i = \frac {p(O, i_1=i| \overline \lambda)} {-\lambda}$

带入 $\sum_{i=1}^N \pi_{i_1}=1$的：

$-\lambda = \sum_{i=1}^N p(O, i_1=i| \overline \lambda) = p(o|\overline \lambda)$

求得并有公式(14)：

$\pi_i = \frac {p(O, i_1=i| \overline \lambda)} {p(o|\overline \lambda)} = \gamma_1(i) \tag{18}$

2.求解 $a_{ij}$:

$L(a_{ij})=\sum_I (\sum_{t=1}^{T-1} a_{i_ti_{t+1}}) p(O, I| \overline \lambda) = \sum_{i=1}^N (\sum_{t=1}^{T-1} a_{i_ti_{t+1}}) ( \sum_{j=1}^N p(O, i_t=i, i_{t+1}=j| \overline \lambda) ) \\ = \sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T-1} a_{ij} p(O, i_t=i, i_{t+1}=j| \overline \lambda)$

应用约束条件 $\sum_{j=1}^N a_{ij} = 1$,用拉格朗日乘子法能够求出：

$\sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T-1} a_{ij} p(O, i_t=i, i_{t+1}=j| \overline \lambda) + \lambda(\sum_{j=1}^N a_{ij} - 1)$

对上式求骗到并等于0获得：

$\frac {\partial}{\partial a_{ij}} [\sum_{i=1}^N \sum_{j=1}^N \sum_{t=1}^{T-1} a_{ij} p(O, i_t=i, i_{t+1}=j| \overline \lambda) + \lambda(\sum_{j=1}^N a_{ij} - 1)] = 0$

获得：

$\sum_{t=1}^{T-1} p(O, i_t=i, i_{t+1}=j| \overline \lambda) + \lambda a_{ij} = 0$

因此：

$a_{ij} = \frac {\sum_{t=1}^{T-1} p(O, i_t=i, i_{t+1}=j| \overline \lambda)}{- \lambda}$

将上式带入 $\sum_{j=1}^N a_{ij} = 1$：

$- \lambda = \sum_{j=1}^N \sum_{t=1}^{T-1} p(O, i_t=i, i_{t+1}=j| \overline \lambda) = \sum_{t=1}^{T-1} p(O, i_t=i| \overline \lambda)$

故得：

$a_{ij} = \frac {\sum_{t=1}^{T-1} p(O, i_t=i, i_{t+1}=j| \overline \lambda)}{\sum_{t=1}^{T-1} p(O, i_t=i| \overline \lambda) } = \frac {\sum_{t=1}^{T-1} p(O, i_t=i, i_{t+1}=j| \overline \lambda) / p(o|\overline \lambda)} {\sum_{t=1}^{T-1} p(O, i_t=i| \overline \lambda) / p(o|\overline \lambda) }$

将（14）和（15）带入的：

$a_{ij} = \frac {\sum_{t=1}^{T-1} \xi_t(i,j)} {\sum_{t=1}^{T-1} \gamma_t(i) } \tag{19}$

3.求解 $b_j{k}$:

$L(b_j{k}) = \sum_I (\sum_{t=1}^T b_{it}(o_t)) p(O, I| \overline \lambda) = \sum_{j=1}^N \sum_{t=1}^T b_{j}(o_t) p(O, i_t=j| \overline \lambda)$

在约束条件 $\sum_{k=1}^M b_j(k) = 1$的拉格朗日乘子法：

$\sum_{j=1}^N \sum_{t=1}^T b_{j}(o_t) p(O, i_t=j| \overline \lambda) + \lambda(\sum_{k=1}^M b_j(k) - 1)$

对其求偏导得：

$\frac {\partial}{\partial b_j(k)} [\sum_{j=1}^N \sum_{t=1}^T b_{j}(o_t) p(O, i_t=j| \overline \lambda) + \lambda(\sum_{k=1}^M b_j(k) - 1)] = 0$

由于只有在 $o_t=v_k$时偏导才不会等于0，以 $I(o_t=v_k)$表示，则：

$\sum_{t=1}^T p(O, i_t=j| \overline \lambda) I(o_t=v_k) + \lambda b_{j}(o_t)I(o_t=v_k) = 0$

$b_{j}(o_t)I(o_t=v_k)$能够写做 $b_{j}(k)$，故：

$b_{j}(k) = \frac {\sum_{t=1}^T p(O, i_t=j| \overline \lambda) I(o_t=v_k)} {- \lambda}$

将上式带入 $\sum_{k=1}^M b_j(k) = 1$得：

$- \lambda = \sum_{k=1}^M \sum_{t=1}^T p(O, i_t=j| \overline \lambda) I(o_t=v_k) = \sum_{t=1}^T p(O, i_t=j| \overline \lambda)$

获得：

$b_{j}(k) = \frac {\sum_{t=1}^T p(O, i_t=j| \overline \lambda) I(o_t=v_k)} {\sum_{t=1}^T p(O, i_t=j| \overline \lambda)}$

又有（14）式可得：

$b_{j}(k) = \frac {\sum_{t=1,o_t=v_k}^T \gamma_t(j)} {\sum_{t=1}^T \gamma_t(j)} \tag{20}$

EM算法总结：
E步：

$\gamma_t(i) = p(i_t=q_i|O,\lambda) = \frac {p(i_t=q_i,O|\lambda)}{p(O|\lambda)} = \frac {\alpha_t(i) \beta_t(i)}{\sum_{j=1}^N \alpha_t(j) \beta_t(j)}$

$\xi_t(i,j) = \frac {p(i_t=q_i,i_{t+1}=q_j|O,\lambda)}{p(O|\lambda)} =\frac {\alpha_t(i) a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)} {\sum_{i=1}^N \sum_{j=1}^N \alpha_t(i) a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)}$

M步：
$\pi_i = \frac {p(O, i_1=i| \overline \lambda)} {p(o|\overline \lambda)} = \gamma_1(i)$

$a_{ij} = \frac {\sum_{t=1}^{T-1} \xi_t(i,j)} {\sum_{t=1}^{T-1} \gamma_t(i) }$

$b_{j}(k) = \frac {\sum_{t=1,o_t=v_k}^T \gamma_t(j)} {\sum_{t=1}^T \gamma_t(j)}$

4. 预测问题（解码问题）

用维特比算法进行求解：
已知：模型 $\lambda=(A,B,\pi)$和 $O=(o_1,o_2,...,o_T)$
求：条件几率最大 $p(I|O,\lambda)$最大的状态序列 $I=(i_1,i_2,...,i_T)$
由于 $p(O)$是一个定值，因此：

$\max_I p(I|O,\lambda) = \max_I p(I, O|\lambda) / p(O|\lambda) = \max_I p(I, O|\lambda)$

定义在时刻 $t$状态为 $i$的全部单个路径 $(i_1,i_2,...,i_t)$中几率最大值为：

$\delta_t(i) = \max_{i_1,i_2,...,i_{t-1}} p(i_t=i, i_{t-1:i_1},o_{t:1}|\lambda)$

递推推导：

$\begin{aligned} &p(i_{t+1}=i,i_{t:1},o_{t+1:1}| \lambda) \\ &=p(i_{t+1}=i,i_t,i_{t-1:1},o_{t+1},o_{t:1}| \lambda) \\ &= p(o_{t+1}|i_{t+1}=i,i_t,o_{t:1},\lambda) p(i_{t+1}=i|i_t,i_{t-1:1},o_{t:1}, \lambda) p(i_t,i_{t-1:1},o_{t:1}|\lambda) \\ &= p(o_{t+1}|i_{t+1}=i,\lambda) p(i_{t+1}=i|i_t,\lambda) p(i_t,i_{t-1:1},o_{t:1}|\lambda) \\ \end{aligned}$

故：

$\delta_{t+1}(i) = \max_{i_1,i_2,...,i_{t-1}} p(i_{t+1}=i,i_{t:1},o_{t+1:1}| \lambda) = \max_{1 \le j \le N} [\delta _t(j) a_{ji}] b_i(o_{t+1}) \tag{21}$

定义在时刻 $t$状态为 $i$的全部单个路径 $(i_1,i_2,...,i_{t-1})$中几率最大的第 $t-1$