隐马尔可夫模型HMM推导

机器学习-白板推导系列(十四)-隐马尔可夫模型HMM（Hidden Markov Model） 课程笔记

背景介绍

介绍一下频率派和贝叶斯派两大流派发展出的建模方式。

频率派

频率派逐渐发展成了统计机器学习，该流派通常将任务建模为一个优化问题，其中有三个关键：

• 模型（model），建立用于解决问题的模型，如 $f(w)=w^T+b$ ；
• 策略（strategy），指导模型优化的策略，比如计算损失函数
• 算法（algorithm），根据策略优化模型的算法，比如 GD、SGD、牛顿法、拟牛顿法

著作有李航老师的统计学习方法。

贝叶斯派

贝叶斯派逐渐发展出了概率图模型，最终是进行推断（inference），也就是求后验概率 $p(z|x)$ 及其期望、方差，最终是一个积分问题。这可以通过一系列数值积分方法（如 MCMC）来计算。

概率图模型

• 有向 —> 贝叶斯网络
• 无向 —> 马尔克夫随机场（马尔可夫网络）
• 概率图模型 + 时间序列 —> 动态模型（Dynamic Model）
  • 隐马儿可夫模型（HMM），隐状态是离散的；
  • 卡尔曼滤波（Kalman Filter），隐状态连续、线性；
  • 粒子滤波（Particle Filter），隐状态连续、非线性。

动态模型：在 GMM 等模型中，假设每个样本是独立同分布的，即不同样本之间没有任何关系。而在动态模型中，在此基础上有了时间序列。这里所谓的时间序列是广义的，可以是真正的时间序列，如语音信号，也可以是其他有序序列，如自然语言中的一个句子。在动态模型中，样本点之间是不满足独立同分布假设的。动态模型中，通常有隐藏状态和观测变量两种随机变量。如 HMM 示意图中，从横向来看，是时间序列（time），从纵向来看，是混合变量（mixture）。

HMM模型表示

HMM 模型的概率图表示如下所示，其中空白节点表示隐藏状态，用 $I=i_1,i_2,\dots ,i_t,\dots$ 来表示，其取值的可能为 $Q=q_1,q_2,\dots ,q_N$，阴影节点表示观测变量，用 $O=o_1,o_2,\dots ,o_t,\dots$ 来表示，其取值的可能为 $V=v_1,v_2,\dots ,v_M$。概率图模型的参数为 $\lambda =(\pi ,A,B)$ ，其中

• $\pi$ 是初始概率分布，即第一个隐状态 $i_1$ 的概率分布，作为隐状态，它的其取值可能有 $N$ 种
• $A=[a_{ij}]$ 是状态转移矩阵（transition matrix），矩阵中的元素 $a_{ij}$ 表示从隐状态 $i$ 到隐状态 $j$ 的转移概率，即 $a_{ij}=P(i_{t+1}=q_j|i_t=q_i)$
• $B=[b_j(k)]$ 是发射矩阵（emission matrix），矩阵中的元素 $b_j(k)$ 表示同一时间点从隐状态 $j$ 到观测变量 $k$ 的发射概率，即 $b_j(k)=P(o_t=v_k|i_t=q_j)$

两个假设

HMM 模型中有两个假设，分别是齐次马尔可夫假设和观测独立假设。

齐次马尔可夫假设，即无后效性，未来状态与过去状态无关，只与当前状态有关。即：
$$P(i_{t+1}|i_t,i_{t-1},\dots ,i_1,o_t,o_{t-1},\dots ,o_1)=P(i_{t+1}|i_t)$$
观测独立假设，某时刻的观测变量只与其自身对应的隐状态有关，与之前的其他变量都无关。即：
$$P(o_t|i_t,i_{t-1},\dots ,i_1,o_{t-1},\dots ,o_1)=P(o_t|i_t)$$

三个问题

HMM 模型中，通常有三种问题：evaluation、learning、decoding。

• evaluation：求值（观测变量）

  已知 λ={π,A,B}\lambda=\{\pi, A, B\}λ={π,A,B} ，求观测变量序列出现的概率，即求 $P(o|\lambda )$ 。通常用前向后向算法（forward-backward algorithm）求解。

• learning：求参数，参数估计

  已知观测序列，求 $\lambda$ ，通常用 EM 算法（baum welch algorithm）求解极大似然估计 $\hat{\lambda }=\arg {\max }_{\lambda }P(O|\lambda )$ 。

• decoding：求状态序列

  已知观测序列和模型参数，求使得该观测序列出现概率最大的隐状态序列，通常用维特比算法（viterbi algorithm）估计 $\hat{I}=\arg {\max }_{i}P(O|I)$ 。

  decoding 问题又可引申出两种问题：prediction 和 filtering。

  • prediction：已知过去和现在的观测序列，求下一时刻的隐状态，即求 $P(i_{t+1}|o_1,o_2,\dots ,o_t)$
  • filtering：已知过去和现在的观测序列，求当前时刻的隐状态，即求 $P(i_t|o_1,o_2,\dots ,o_t)$ ，HMM 一般不做这个问题。

evaluation问题

朴素解法

给定 HMM 模型的参数，求观测变量，即求 $P(O|\lambda )$ 。
$$P(O|\lambda )=\sum _{S}P(S,O|\lambda )=\sum _{S}P(O|I,\lambda )\cdot P(I|\lambda )$$
其中，后一项：
$$\begin{array}{rl}P(S|\lambda )& =P(i_1,i_2,\dots ,i_T|\lambda )\\ & =P(i_T|i_1,i_2,\dots ,i_{T-1},\lambda )\cdot P(i_1,i_2,\dots ,i_{T-1},\lambda )\\ & =P(i_T|i_{T-1})\cdot P(i_1,i_2,\dots ,i_{T-1},\lambda )齐次马尔可夫假设\\ & =P(i_T|i_{T-1})P(i_{T-1}|i_{T-2})\dots P(i_2|i_1)递归展开\\ & =a_{i_{T-1},i_T}{a}_{i_{T-2},i_{T-1}}\dots a_{12}\\ & =\prod _{t=2}^T{a}_{i_{t-1}{i}_t}\end{array}$$
同理，根据观测独立假设，前一项：
$$P(O|I,\lambda )=\prod _{t=1}^T{b}_{i_t}(o_t)$$
从而：
$$\begin{array}{rl}P(O|\lambda )& =\sum _I\prod _{t=2}^T{a}_{i_{t-1}{i}_t}\prod _{t=1}^T{b}_{i_t}(o_t)\\ & =\sum _{i_1}\sum _{i_2}\cdots \sum _{i_T}\prod _{t=2}^T{a}_{i_{t-1}{i}_t}\prod _{t=1}^T{b}_{i_t}(o_t)\end{array}$$
上式的时间复杂度为 $\mathcal{O}(T{N}^T)$ ，无疑是不能接受的。我们需要寻找时间复杂度更低的算法。

前向算法

推导前向后向算法的过程中用到的变换主要就是 HMM 的两个假设，以及条件概率的公式。

记 ${\alpha }_t(i)=P(o_1,\dots ,o_t,i_t=q_i|\lambda )$ ，则
$$P(O|\lambda )=\sum _{i=1}^{N}P(O,i_t=q_i|\lambda )=\sum _{i=1}^N{\alpha }_T(i)$$
要求观测变量 $P(O|\lambda )$， 现在我们只要求出 ${\alpha }_T(i)$ ，然后再 $i=1\dots N$ 上求和即可。

一个自然的想法就是找 ${\alpha }_t(i)$ 的递推式，即其与 ${\alpha }_{t+1}(i)$ 的关系。有：
$$\begin{array}{rl}{\alpha }_{t+1}(j)& =P(o_1,\dots ,o_t,o_{t+1},i_{t+1}=q_j|\lambda )\\ & =\sum _{i=1}^{N}P(o_1,\dots ,o_t,o_{t+1},i_t=q_i,i_{t+1}=q_j|\lambda )\\ & =\sum _{i=1}^{N}P(o_{t+1}|o_1,\dots ,o_t,i_t=q_i,i_{t+1}=q_j,\lambda )\cdot P(o_1,\dots ,o_t,i_t=q_i,i_{t+1}=q_j|\lambda )\\ & =\sum _{i=1}^{N}P(o_{i+1}|i_{t+1}=q_j)\cdot P(o_1,\dots ,o_t,i_t=q_i,i_{t+1}=q_j|\lambda )观测独立假设\\ & =\sum _{i=1}^{N}P(o_{i+1}|i_{t+1}=q_j)\cdot P(i_{t+1}=q_j|o_1,\dots ,o_t,i_t=q_i,\lambda )\cdot P(o_1,\dots ,o_t,i_t=q_i|\lambda )\\ & =\sum _{i=1}^{N}P(o_{i+1}|i_{t+1}=q_j)\cdot P(i_{t+1}=q_j|i_t=q_i,\lambda )\cdot P(o_1,\dots ,o_t,i_t=q_i|\lambda )齐次马尔可夫假设\\ & =\sum _{i=1}^{N}P(o_{i+1}|i_{t+1}=q_j)\cdot P(i_{t+1}=q_j|i_t=q_i,\lambda )\cdot {\alpha }_t(i)\\ & =\sum _{i=1}^N{b}_j(o_{t+1})\cdot a_{ij}\cdot {\alpha }_t(i)\end{array}$$
由此，我们就得到了递推式：
$${\alpha }_{t+1}(j)=\sum _{i=1}^N{b}_j(o_{t+1})\cdot a_{ij}\cdot {\alpha }_t(i)$$
从而可以前向（从前往后，从 $1$ 到 $T$）求得 ${\alpha }_T(i)$ ，从而求得 $P(O|\lambda )={\sum }_{i=1}^N{\alpha }_T(i)$ 。

后向算法

与前向算法类似，我们定义一个记号 ${\beta }_t(i)$ ：
$${\beta }_t(i)=P(o_t,\dots ,o_T|i_t=q_i,\lambda )$$
则：
$$\begin{array}{rl}P(O|\lambda )& =P(o_1,\dots ,o_T|\lambda )\\ & =\sum _{i=1}^{N}P(o_1,\dots ,o_T,i_1=q_i|\lambda )\\ & =\sum _{i=1}^{N}P(o_1,\dots ,o_T|i_1=q_i,\lambda )\cdot P(i_1=q_i)\\ & =\sum _{i=1}^{N}P(o_1,\dots ,o_T|i_1=q_i,\lambda )\cdot {\pi }_i\\ & =\sum _{i=1}^{N}P(o_1|o_2,\dots ,o_T,i_1=q_i,\lambda )\cdot P(o_2,\dots ,o_T|i_1=q_i,\lambda )\cdot {\pi }_i\\ & =\sum _{i=1}^{N}P(o_1|i_1=q_i,\lambda )\cdot {\beta }_1(i)\cdot {\pi }_i\\ & =\sum _{i=1}^N{b}_i(o_1)\cdot {\pi }_i\cdot {\beta }_1(i)\end{array}$$
$\lambda$ 为已知，下面就省略不写了。同样我们试图得到递推式，即 ${\beta }_{t+1}(j)$ 与 ${\beta }_t(i)$ 之间的关系。有：
$$\begin{array}{rl}{\beta }_t(i)& =P(o_{t+1},\dots ,o_T|i_t=q_i)\\ & =\sum _{j=1}^{N}P(o_{t+1},\dots ,o_T,i_{t+1}=q_j|i_t=q_i)\\ & =\sum _{j=1}^{N}P(o_{t+1},\dots ,o_T|i_{t+1}=q_j,i_t=q_i)\cdot P(i_{t+1}=q_j|i_t=q_i)\\ & =\sum _{j=1}^{N}P(o_{t+1},\dots ,o_T|i_{t+1}=q_j)\cdot a_{ij}前面拿掉i_t=q_i的证明暂略\\ & =\sum _{j=1}^{N}P(o_{t+1}|o_{t+2},\dots ,o_T,i_{t+1}=q_j)\cdot P(o_{t+2},\dots ,o_T|i_{t+1}=q_j)\cdot a_{ij}\\ & =\sum _{j=1}^{N}P(o_{t+1}|i_{t+1}=q_j)\cdot {\beta }_{t+1}(j)\cdot a_{ij}独立观测假设\\ & =\sum _{j=1}^N{b}_j(o_{t+1})\cdot a_{ij}\cdot {\beta }_{t+1}(j)\end{array}$$
由此，就得到了递推式：
$${\beta }_t(i)=\sum _{j=1}^N{b}_j(o_{t+1})\cdot a_{ij}\cdot {\beta }_{t+1}(j)$$
从而，可以后向（从后往前，从 $T$ 到 $1$）求得 ${\beta }_1(i)$ ，从而求得 $P(O|\lambda )={\sum }_{i=1}^N{b}_i(o_1)\cdot {\pi }_i\cdot {\beta }_1(i)$ 。

总结

	递推方向	设定	求解公式	复杂度
朴素解法	-	-	$P(O	\lambda)=\sum_{i_1}\sum_{i_2}\dots\sum_{i_T}\prod_{t=2}Ta_{i_{t-1}i_t}\prod_{t=1}Tb_{i_t}(o_t)$
前向算法	从前向后，从$1$到$T$	$\alpha_t(i)=P(o_1,\dots,o_t,i_t=q_i	\lambda)$	$P(O
后向算法	从后向前，从$T$到$1$	$\beta_t(i)=P(o_t,\dots,o_T	i_t=q_i,\lambda)$	$P(O

learning问题

HMM 的 learning 问题是要估计出其参数 $\lambda =(\pi ,A,B)$ 。如果直接用极大似然估计，即：
$${\lambda }_{MLE}=\arg \underset{\lambda }{\max }P(O|\lambda )$$
是无法直接解析的。因此，我们这里用迭代的 EM 算法来估计参数。EM 算法的迭代公式如下：
$${\theta }^{(t+1)}=\arg \underset{\theta }{\max }{\int }_Z\log P(X,Z|\theta )\cdot P(Z|X,{\theta }^{(t)})dZ$$
其中 $x$ 是观测变量，$z$ 是隐变量，$\theta$ 是参数。之前提到，HMM 除了可以横向看做是时间（time）序列之外，还可以纵向看做是混合（mixture）变量。将观测变量、隐变量和参数分别对应到 HMM 中符号之后，顺便做一步简化，把常数拿掉，有：
$$\begin{array}{rl}{\lambda }^{(t+1)}& =\arg \underset{\lambda }{\max }\sum _I\log P(O,I|\lambda )\cdot P(I|O,{\lambda }^{(t)})\\ & =\arg \underset{\lambda }{\max }\sum _I\log P(O,I|\lambda )\cdot \frac{P(O,I|{\lambda }^{(t)})}{P(O|{\lambda }^{(t)})}\\ & =\arg \underset{\lambda }{\max }\sum _I\log P(O,I|\lambda )\cdot P(O,I|{\lambda }^{(t)})\end{array}$$
记为 Q 函数 $Q(\lambda ,{\lambda }^{(t)})$ ，记 ${\lambda }^{(t)}=({\pi }^{(t)},A^{(t)},B^{(t)})$ ，又有 $P(O|\lambda )={\sum }_{I}P(O,I|\lambda )={\sum }_{i_1}\cdots {\sum }_{i_T}{\pi }_{i_1}{\prod }_{t=2}^T{a}_{i_{t-1},i_t}{\prod }_{t=1}^T{b}_{i_t}(o_t)$，则：
$$\begin{array}{rl}Q(\lambda ,{\lambda }^{(t+1)})& =\sum _I\log P(O,I|\lambda )\cdot P(O,I|{\lambda }^{(t)})\\ & =\sum _I[(\log {\pi }_{i_1}+\sum _{t=2}^T\log a_{i_{t-1},i_t}+\sum _{t=1}^T\log b_{i_t}(o_t))\cdot P(O,I|{\lambda }^{(t)})]\end{array}$$
这里就只介绍 $\pi$ 的参数估计，$A,B$ 类似。观察上面 Q 函数的展开结果，只有第一项与 $\pi$ 相关，有：
$$\begin{array}{rl}{\pi }^{(t+1)}& =\arg \underset{\pi }{\max }Q(\lambda ,{\lambda }^{(t)})\\ & =\arg \underset{\pi }{\max }\sum _I\log {\pi }_{i_1}\cdot P(O,I|{\lambda }^{(t)})\\ & =\arg \underset{\pi }{\max }\sum _{i_1}\cdots \sum _{i_T}[\log {\pi }_{i_1}\cdot P(O,i_1,\dots ,i_T|{\lambda }^{(t)})]\\ & =\arg \underset{\pi }{\max }\sum _{i=1}^N[\log {\pi }_{i}P(O,i_1=q_i|{\lambda }^{(t)})]\end{array}$$
由于 $\pi$ 是初始状态的概率分布，因此上述优化问题有约束：${\sum }_{i=1}^N{\pi }_i=1$ 。解决带约束的优化问题，这里用拉格朗日乘子法。

令：
$$\mathcal{L}(\pi ,\eta )=\sum _{i=1}^N\log {\pi }_{i}P(O,i_1=q_i|{\lambda }^{(t)})+\eta (\sum _{i=1}^N{\pi }_i-1)$$
求其对 ${\pi }_i$ 的偏导，并令其为 0：
$$\begin{gathered} \frac{\partial \mathcal{L}}{\partial {\pi }_i}=\frac{1}{{\pi }_i}P(O,i_1=q_i|{\lambda }^{(t)})+\eta =0 \\ \sum _{i=1}^N[P(O,i_1=q_i|{\lambda }^{(t)})+{\pi }_i\eta ]=0 \\ P(O|{\lambda }^{(t)})+\eta =0 \\ \eta =-P(O|{\lambda }^{(t)}) \end{gathered}$$
将其带回，有：
$${\pi }_i^{(t+1)}=-\frac{1}{\eta }P(O|i_1=q_i|{\lambda }^{(t)})=\frac{P(O,i_1=q_i|{\lambda }^{(t)})}{P(O|{\lambda }^{(t)})}$$
至此，我们就得到了由 ${\pi }_i^{(t)}$ 迭代计算 ${\pi }_i^{(t+1)}$ 的公式，从而，就可以估计 ${\pi }^{(t+1)}$：
$${\pi }^{(t+1)}=({\pi }_1^{(t+1)},\dots ,{\pi }_N^{(t+1)})$$
这就完成了 HMM 参数中 $\pi$ 的 learning 问题，对于另外两个参数 $A,B$ ，方法类似。这就是 baum-welch 算法，可以看做是一种特殊形式的 EM 算法。

decoding问题

decoding 问题是已知观测序列和模型参数的情况下，求使得该观测序列出现概率最大的隐状态序列，通常用维特比算法（viterbi algorithm）估计 $\hat{I}=\arg {\max }_{i}P(O|I)$ 。

在之前的题目设定中有：隐状态 $i_t$ 取值的集合为 $Q=q_1,q_2,\dots ,q_N$ ，每个时间点的隐状态 $i_t$ 都有可能取这些值，我们就是要每一步 $i_t$ 选取一个 $q_i$ ，最终找到一组 $I=i_1,\dots ,i_T$ 的取值（图中蓝色虚线），使得观测变量序列 $O=o_1,\dots ,o_T$ 出现的概率最大。

实际中，我们使用动态规划的思想来解决该问题。首先定义：${\delta }_t(i)={\max }_{i_1,i_2,\dots ,i_{t-1}}P(o_1,\dots ,o_t.i_1,\dots ,i_{t-1},i_t=q_i)$ 表示到达第 $i_t$ 个位置取值为 $i_t=q_i$ 的最大概率值。接下来我们要找到动态规划中的状态转移方程，即 ${\delta }_{t+1}(j)$ 与 ${\delta }_t(i)$ 的关系。其实就是乘上从 $i$ 到 $j$ 对应的转移概率 $a_{ij}$ 和 $t+1$ 时刻的发射概率 $b_j(o_{t+1})$ 和 ${\delta }_t(i)$ 乘积的最大者，即有：
$$\begin{array}{rl}{\delta }_{t+1}(j)& =\underset{i_1,\dots ,i_t}{\max }P(o_1,\dots ,o_{t+1},i_1,\dots ,i_t,i_{t+1}=q_j)\\ & =\underset{1\le i\le N}{\max }{\delta }_t(i)a_{ij}{b}_j(o_{t+1})\end{array}$$
这样，我们就求出了每一步的最大概率值，注意现在只是求得了最大概率值，没有记录路径，而我们的目标是要求最大概率值对应的隐状态序列路径。

令
$${\varphi }_{t+1}(j)=\arg \underset{1\le i\le N}{\max }{\delta }_t(i)\cdot a_{ij}$$
即可将每一步所选的 $q_i$ 记录下来。这就是求解 HMM decoding 问题的维特比算法。

总结

HMM 是一种动态模型，即含有时间序列这一维度。更广义地来说，它是一种状态空间模型（State Space Model），状态空间模型可以不包含时间维度，如卡尔曼滤波、粒子滤波等，它们的概率图都类似（如图 1 所示）。

记隐变量为 $Z$，观测变量为 $X$，模型参数为 $\theta$ ，这一类模型的问题可以归纳为下面几种：

• Learning，根据观测变量估计模型参数，如由极大似然来估计参数 $\hat{\theta }=\arg {\max }_{\theta }P(X|\theta )$，可由 Baum Welch（EM）算法求解；
• Inference，根据模型参数进行分析、求值、预测等，以下 $\theta$ 均已知
  • Decoding，根据观测变量解码出隐变量，即求 $\hat{Z}=\arg {\max }_{Z}P(Z|X)$，可由维特比算法（动态规划）求解；
  • Evaluation，已知模型参数计算观测变量，即求 prob of evidence $P(X)$ 。通常如果使用朴素方法直接计算，会有较高的复杂度，因此对该问题的求解也有对应的算法，如前向算法和后向算法；
  • Filtering，根据历史观测变量求当前时刻的隐变量，即求 $P(z_t|x_1,x_2,\dots ,x_t)$ 。比起只根据当前时刻观测变量求隐变量 $P(z_t|x_t)$ ，滤波能够根据更多历史信息，滤除噪声，故称为 Filtering。常用于在线（online）分析。也可由前向算法求解。
  • Smoothing，根据全部观测变量，求某点的隐变量，即求 $P(z_t|x_1,x_2,\dots ,x_T)$ ，常用于离线（offline）分析。可由前向后向算法求解。
  • Prediction，根据历史观测变量，求下一个或下两个时刻的隐变量或观测变量，即求 $P(z_{t+1}|x_1,x_2,\dots ,x_t)$ 或 $P(x_{t+1}|x_1,x_2,\dots ,x_t)$ 。比如 NLP 中根据 ”我爱中国共“ 这 5 个历史观测变量，预测出下两个观测变量为 ”产党“。可用前向算法求解。