搞懂HMM

1 概述

本文是B站上机器学习-白板推导系列(十四)-隐马尔可夫模型HMM的学习笔记，UP主讲得实在是太清楚了，赶紧记录下来，以防之后忘记。

某些细节上根据个人理解做了改动。

图1 HMM示意图

HMM全名为Hidden Markov Model，其示意图如上图所示，是一个概率图。观测变量，顾名思义，就是我们观测到的量，比如语音识别里就是我们听到的声音信号；状态变量就是隐藏的特征，在语音识别里，可以是发音单元Phoneme，甚至是更小的单元Tri-phone，不管是什么，这必须要是一个离散的可枚举的集合。当状态变量变成连续变量的时候，如果连续变量是线性的，典型的代表就是Kalman Filter，如果是连续变量是非线性的，典型的代表就是Particle Filter。本文只讲HMM。

同一个时刻，从状态变量转变为观测变量，服从某个分布，一般是混合高斯分布（GMM）。

状态变量从前一个时刻转变为下一个时刻的，也服从某个分布，一般也是GMM。

每个时刻的观测变量之间，必须不是独立同分布的。

2 符号说明

令整个序列共有T个time step。

状态序列为$S=[s_1,s_2,...,s_t,...,s_{T-1},s_T]$，$s_t$是可枚举的离散变量，值域为 { q 1 , q 2 , . . . , q N } \{q_1, q_2, ..., q_N\} {q1​,q2​,...,qN​}，$N$表示有$N$种状态。

观测序列为$O=[o_1,o_2,...,o_t,...,o_{T-1},o_T]$，$o_t$可以是连续变量。

${\pi }_i$为初始时刻的状态概率，即${\pi }_i=P(s_1=q_i)$。

$A$为状态转移矩阵$[a_{ij}{]}_{N\times N}$，矩阵中$a_{ij}=P(s_t=q_j|s_{t-1}=q_i)$，表示任意两个相邻时间点时间状态从$q_i$转变为$q_j$的概率。

$b_j(o_t)$为发射概率，表示从状态$q_j$变成观测值$o_t$的概率，即$b_j(o_t)=P(o_t|s_t=q_j)$。

$\lambda =(\pi ,a,b)$表示模型中所有的可学习的参数。

有了符号的图1就变成了

图2 HMM示意图（有符号）

3 两点假设

（1）齐次Markov假设
$t+1$时刻的状态只和$t$时刻的状态有关。
$$P(s_{t+1}|s_1,s_2,...,s_t,o_1,o_2,...,o_t)=P(s_{t+1}|s_t)\text{\hspace{1em}(3-1)}$$

（2）观测独立假设
$t$时刻的观测变量只和$t$时刻的状态变量有关。
$$P(o_t|s_1,s_2,...,s_t,o_1,o_2,...,o_{t-1})=P(o_t|s_t)\text{\hspace{1em}(3-2)}$$

这两个假设在后面的推导中有着极其重要的作用。

4 Evaluation

Evaluation要做的事情是，给定所有的模型参数，即$\lambda$，之后，求得到某个观测序列$O=[o_1,o_2,...,o_T]$的概率，记作$P(O|\lambda )$。注意，现在我们是知道模型所有的参数的，只是做一个inference的过程。

我们先来看下直接求解的情况是如何的。我们把这个条件概率稍微变一下，把状态变量给引进来

$$P(O|\lambda )=\sum _{allS}P(S,O|\lambda )=\sum _{allS}P(O|S,\lambda )P(S|\lambda )\text{\hspace{1em}(4-1)}$$

这个没问题吧，我们把所有可能的$S$序列都考虑进来了，这是个全概率。

然后我们把$P(S|\lambda )$展开来看下，$\lambda$只表示已知所有模型参数，可写可不写

$$\begin{array}{rl}P(S|\lambda )& =P(s_1,s_2,...,s_T|\lambda )\\ & =P(s_T|s_1,s_2,...,s_{T-1},\lambda )P(s_1,s_2,...,s_{T-1},\lambda )\\ & 利用齐次Markov假设(3-1)，\lambda 可写可不写\\ & =P(s_T|s_{T-1})P(s_1,s_2,...,s_{T-1})\\ & 继续拆拆拆\\ & =P(s_T|s_{T-1})P(s_{T-1}|s_{T-2})...P(s_2|s_1)P(s_1)\\ & 除了最后一项都是状态转移矩阵里的\\ & =\prod _{t=1}^{T-1}{a}_{s_t,s_{t+1}}{\pi }_{s_1}\end{array}\text{\hspace{1em}(4-2)}$$

接着我们再把$P(O|S,\lambda )$展开看下

$$\begin{array}{rl}P(O|S,\lambda )& =P(o_1,o_2,...,o_T|s_1,s_2,...,s_T,\lambda )\\ & =P(o_T|o_1,o_2,...,o_{T-1},s_1,s_2,...,s_T,\lambda )P(o_1,o_2,...,o_{T-1}|s_1,s_2,...,s_T,\lambda )\\ & 利用观测独立假设(3-2)\\ & =P(o_T|s_T)P(o_1,o_2,...,o_{T-1}|s_1,s_2,...,s_T,\lambda )\\ & 继续拆拆拆\\ & =P(o_T|s_T)P(o_{T-1}|s_{T-1})...P(o_1|s_1)\\ & 用发射概率函数\\ & =\prod _{t=1}^T{b}_{s_t}(o_t)\end{array}\text{\hspace{1em}(4-3)}$$

将(4-2)和(4-3)带入(4-1)可得

$$\begin{array}{rl}P(O|\lambda )& =\sum _{allS}\prod _{t=1}^T{b}_{s_t}(o_t)\prod _{t=1}^{T-1}{a}_{s_t,s_{t+1}}{\pi }_{s_1}\\ & 把allS展开\\ & =\sum _{s_1=q_1}^{q_N}\sum _{s_2=q_1}^{q_N}...\sum _{s_T=q_1}^{q_N}\prod _{t=1}^T{b}_{s_t}(o_t)\prod _{t=1}^{T-1}{a}_{s_t,s_{t+1}}{\pi }_{s_1}\end{array}\text{\hspace{1em}(4-4)}$$

这个的复杂度是$O(N^T)$的，计算量随着序列的变长而爆炸式指数增长，没法用的。

于是，就有人提出了前向和后向算法来降低计算成本。

4.1 前向算法（forward algorithm）

如下图3所示，前向算法考虑的就是橘黄色方框内变量的联合概率，记作

$${\alpha }_t(q_i)=P(o_1,o_2,...,o_t,s_t=q_i|\lambda )\text{\hspace{1em}(4-5)}$$

若要问为什么要是$(4-5)$，那我还真答不上来，这是一种设计，如果有其他的设计应该也可以。

图3 前向算法示意图

我们来看看${\alpha }_T(q_i)$是怎么样的

$${\alpha }_T(q_i)=P(o_1,o_2,...,o_T,s_T=q_i|\lambda )=P(O,s_T=q_i|\lambda )\text{\hspace{1em}(4-6)}$$

$s_T$的状态是可枚举的，我们遍历所有的$s_T$的可能性，然后求个和，就有了

$$\sum _{i=1}^N{\alpha }_T(q_i)=\sum _{i=1}^{N}P(O,s_T=q_i|\lambda )=P(O|\lambda )\text{\hspace{1em}(4-7)}$$

这不，$P(O|\lambda )$出现了。

接下来我们的工作就是，看看怎么求这个${\alpha }_t(q_i)$。${\alpha }_1(q_i)$我们是知道的

$${\alpha }_1(q_i)=P(o_1,s_1=q_i|\lambda )=P(o_1|s_1=q_i,\lambda )P(s_1=q_i|\lambda )$$

看出来了吗？一个是我们的发射概率，一个是我们的初始概率，所以有

$${\alpha }_1(q_i)=b_i(o_1){\pi }_i\text{\hspace{1em}(4-8)}$$

既然知道了${\alpha }_1(q_i)$，那我们如果还能知道${\alpha }_t(q_i)$到${\alpha }_{t+1}(q_i)$的递推公式，这个问题不久解决了吗？我们来试试看！$\lambda$写不写都无所谓，我们心里知道就好，下面就不写了。

$$\begin{array}{rl}{\alpha }_{t+1}(q_i)& =P(o_1,o_2,...,o_{t+1},s_{t+1}=q_i)\\ & 用全概率凑个s_t出来试试\\ & =\sum _{j=1}^{N}P(o_1,o_2,...,o_{t+1},s_t=q_j,s_{t+1}=q_i)\\ & 提o_{t+1}\\ & =\sum _{j=1}^{N}P(o_{t+1}|o_1,o_2,...o_t,s_t=q_j,s_{t+1}=q_i)P(o_1,o_2,...,o_t,s_t=q_j,s_{t+1}=q_i)\\ & 利用观测独立假设(3-2)\\ & =\sum _{j=1}^{N}P(o_{t+1}|s_{t+1}=q_i)P(o_1,o_2,...,o_t,s_t=q_j,s_{t+1}=q_i)\\ & 提后项的s_{t+1}\\ & =\sum _{j=1}^{N}P(o_{t+1}|s_{t+1}=q_i)P(s_{t+1}=q_i|o_1,o_2,...,o_t,s_t=q_j)P(o_1,o_2,...,o_t,s_t=q_j)\\ & 利用齐次Markov假设(3-1)\\ & =\sum _{j=1}^{N}P(o_{t+1}|s_{t+1}=q_i)P(s_{t+1}=q_i|s_t=q_j)P(o_1,o_2,...,o_t,s_t=q_j)\end{array}$$

发现了吗？这三项分别是发射概率，状态转移概率和${\alpha }_t(q_j)$。

于是我们就得到了递推式

$${\alpha }_{t+1}(q_i)=\sum _{j=1}^N{b}_i(o_{t+1})a_{ji}{\alpha }_t(q_j)\text{\hspace{1em}(4-9)}$$
结合$(4-8)$和$(4-9)$我们可以得到所有状态下的${\alpha }_T(q_i)$，$(4-7)$得解，此时的复杂度为$O((TN{)}^2)$。

4.2 后向算法（backward algorithm）

如下图4所示，后向算法考虑的就是青绿色方框内变量的联合概率，记作

$${\beta }_t(q_i)=P(o_{t+1},...,o_{T-1},o_T|s_t=q_i,\lambda )\text{\hspace{1em}(4-10)}$$

这也是一种设计，和前向的互补。注意看与$(4-5)$的区别，后向的推导相比前向要绕一点。

图4 后向算法示意图

我们来看下${\beta }_1(q_i)$是怎么样的

$${\beta }_1(q_i)=P(o_2,...,o_{T-1},o_T|s_1=q_i|\lambda )\text{\hspace{1em}(4-11)}$$

我们来看下这个${\beta }_1(q_i)$和我们要求的$P(O|\lambda )$有什么关系

$$\begin{array}{rl}P(O|\lambda )& =P(o_1,o_2,...,o_T|\lambda )\\ & 省略\lambda ，引入s_1\\ & =\sum _{i=1}^{N}P(o_1,o_2,...,o_T,s_1=q_i)\\ & 把s_1当成条件\\ & =\sum _{i=1}^{N}P(o_1,o_2,...,o_T|s_1=q_i)P(s_1=q_i)\\ & 拆出o_1，注意后向为初始概率\\ & =\sum _{i=1}^{N}P(o_1|o_2,...,o_T,s_1=q_i)P(o_2,...,o_T|s_1=q_i){\pi }_i\\ & 利用观测独立假设(3-2)\\ & =\sum _{i=1}^{N}P(o_1|s_1=q_i)P(o_2,...,o_T|s_1=q_i){\pi }_i\\ & 代入(4-11)和发射概率\\ & =\sum _{i=1}^N{b}_i(o_1){\beta }_1(q_i){\pi }_i\end{array}\text{\hspace{1em}(4-12)}$$

这样以来，$P(O|\lambda )$和${\beta }_1(q_i)$的一个关系就找到了，要做的就是来求一下这个${\beta }_1(q_i)$了。

我们令

$${\beta }_T(q_i)=1\text{\hspace{1em}(4-13)}$$

然后再来算一下${\beta }_t(q_i)$和${\beta }_{t+1}(q_j)$递推关系，$\lambda$我就直接省略了。

$$\begin{array}{rl}{\beta }_t(q_i)& =P(o_{t+1},...,o_{T-1},o_T|s_t=q_i)\\ & 利用全概率引入s_{t+1}\\ & =\sum _{j=1}^{N}P(o_{t+1},...,o_T,s_{t+1}=q_j|s_t=q_i)\\ & 把s_{t+1}引到条件当中去\\ & =\sum _{j=1}^{N}P(o_{t+1},...,o_T|s_{t+1}=q_j,s_t=q_i)P(s_{t+1}=q_j|s_t=q_i)\\ & 前项中的o_{t+1},...,o_T只和s_{t+1}=q_j有关，这个可证，但这里不证\\ & 后项为状态转移概率\\ & =\sum _{j=1}^{N}P(o_{t+1},...,o_T|s_{t+1}=q_j)a_{ij}\\ & 把o_{t+1}拿出来\\ & =\sum _{j=1}^{N}P(o_{t+1}|o_{t+2},...,o_T,s_{t+1}=q_j)P(o_{t+2},...,o_T|s_{t+1}=q_j)a_{ij}\\ & 利用观测独立假设(3-2)\\ & =\sum _{j=1}^{N}P(o_{t+1}|s_{t+1}=q_j){\beta }_{t+1}(q_j)a_{ij}\\ & 前项为发射概率\\ & =\sum _{j=1}^N{b}_j(o_{t+1}){\beta }_{t+1}(q_j)a_{ij}\end{array}\text{\hspace{1em}(4-14)}$$

结合$(4-13)$和$(4-14)$，我们就可以求得${\beta }_1(q_i)$，也就可以求得$P(O|\lambda )$。

值得注意的是，在任意时刻$t$，我们结合前向和后向算法，都可以有

$$P(O|\lambda )=\sum _{i=1}^N{\alpha }_t(q_i){\beta }_t(q_i)\text{\hspace{1em}(4-15)}$$

这个这里简单说下，这个会依赖于一个不靠谱的假设，就是$o_1,...,o_t$和$o_{t+1},...,o_T$是不相关的。

$$\begin{array}{rl}P(O|\lambda )& =\sum _{i=1}^{N}P(O,s_t=q_i|\lambda )\\ & =\sum _{i=1}^{N}P(o_1,...,o_t,o_{t+1},..,o_T,s_t=q_i|\lambda )\\ & =\sum _{i=1}^{N}P(o_1,...,o_t,s_t=q_i|\lambda )P(o_{t+1},..,o_T|o_1,...,o_t,s_t=q_i)\\ & 前项为{\alpha }_t(q_i)，后向依赖于假设忽略o_1,...,o_t就是{\beta }_t(q_i)\\ & =\sum _{i=1}^N{\alpha }_t(q_i){\beta }_t(q_i)\end{array}$$

虽然这个假设不靠谱，但是为了简化计算，都是这么做的。

5 Learning

Learning要做的一件事情就是，在给定了观测序列之后，找到得到该观测序列概率最大的那组参数$\lambda$，即

$${\lambda }_{MLE}=arg\underset{\lambda }{\max }P(O|\lambda )\text{\hspace{1em}(5-1)}$$

这里的MLE就是Max Likelyhood Estimation。

讲道理，如果能够把导数的表达式求出来的话，这个${\lambda }_{MLE}$一下子就出来了，但这里的$P(O|\lambda )$一般都是混合高斯函数，没法直接求导，所以就需要用到EM算法了。

这篇不讲EM算法是什么，我们直接使用EM算法，想知道EM算法是什么的话，推荐看徐亦达老师的EM算法讲解。一句话概括就是，本来没法直接求导的，我们添加了一个隐变量，变成去求另一个方程的导数，这个方程求导相对简单，但是没法一步到位，需要不断迭代，逐渐逼近局部最优点。

EM算法的迭代公式为

$${\theta }^{(t+1)}=arg\underset{\theta }{\max }{\int }_{z}logP(x,z|\theta )P(z|x,{\theta }^{(t)})dz\text{\hspace{1em}(5-2)}$$

这里的$\theta$就是我们的模型参数$\lambda$，$x$就是我们的观测变量$O$，$z$就是我们的状态变量$S$，我们这里的$S$是离散的，积分变累加。我们来改写一下$(5-2)$就有

$${\lambda }^{(t+1)}=arg\underset{\lambda }{\max }\sum _{allS}logP(O,S|\lambda )P(S|O,{\lambda }^{(t)})\text{\hspace{1em}(5-3)}$$

我们这里再对$P(S|O,{\lambda }^{(t)})$做一个小的改动

$$P(S|O,{\lambda }^{(t)})=\frac{P(S,O|{\lambda }^{(t)})}{P(O|{\lambda }^{(t)})}$$

这里的${\lambda }^{(t)}$是一个常数，$O$又是和$\lambda$无关的，所以$P(O|{\lambda }^{(t)})$这项是个常数，可以忽略。故把$(5-3)$改为

$${\lambda }^{(t+1)}=arg\underset{\lambda }{\max }\sum _{allS}logP(O,S|\lambda )P(O,S|{\lambda }^{(t)})\text{\hspace{1em}(5-4)}$$

实际操作时，就是不断迭代$(5-4)$。但看到这里，这个argmax还是不会求呀。这个其实也挺复杂的，下面会以求初始概率参数$\pi$为例，简单说明下，其他的就不求了，太复杂了，吃不消。

我们定义

$$Q(\lambda ,{\lambda }^{(t)})=\sum _{allS}logP(O,S|\lambda )P(O,S|{\lambda }^{(t)})\text{\hspace{1em}(5-5)}$$

我们把式$(4-4)$代进来看下

$$\begin{array}{rl}Q(\lambda ,{\lambda }^{(t)})& =\sum _{s_1=q_1}^{q_N}...\sum _{s_T=q_1}^{q_N}log(\prod _{t=1}^T{b}_{s_t}(o_t)\prod _{t=1}^{T-1}{a}_{s_t,s_{t+1}}{\pi }_{s_1})P(O,S|{\lambda }^{(t)})\\ & =\sum _{s_1=q_1}^{q_N}...\sum _{s_T=q_1}^{q_N}(log{\pi }_{s_1}+\sum _{t=1}^{T}log{b}_{s_t}(o_t)+\sum _{t=1}^{T-1}log{a}_{s_t,s_{t+1}})P(O,S|{\lambda }^{(t)})\end{array}$$

好，我们令第$t$次迭代过程中，初始概率的参数为${\pi }^{(t)}$，那么

$$\begin{array}{rl}{\pi }^{(t+1)}& =arg\underset{\pi }{\max }Q(\lambda ,{\lambda }^{(t)})\\ & 过滤掉和\pi 无关的变量\\ & =arg\underset{\pi }{\max }\sum _{s_1=q_1}^{q_N}...\sum _{s_T=q_1}^{q_N}log{\pi }_{s_1}P(O,s_1,...,s_T|{\lambda }^{(t)})\\ & 和s_1无关的状态变量用全概率去掉\\ & =arg\underset{\pi }{\max }\sum _{s_1=q_1}^{q_N}log{\pi }_{s_1}P(O,s_1|{\lambda }^{(t)})\end{array}\text{\hspace{1em}(5-6)}$$

这样一来，我们就方便去求导了。不过这里别忘了有一个约束条件，就是
$$s.t.\sum _{s_1=q_1}^{q_N}{\pi }_{s_1}=1\text{\hspace{1em}(5-7)}$$

有约束求极值，拉格朗日乘子法来也。令

$$L({\pi }_{s_1},\eta )=\sum _{s_1=q_1}^{q_N}(log{\pi }_{s_1}P(O,s_1|{\lambda }^{(t)}))+\eta (\sum _{s_1=q_1}^{q_N}{\pi }_{s_1}-1)\text{\hspace{1em}(5-8)}$$

对$(5-8)$进行求偏导，有
$$\frac{\partial L}{\partial {\pi }_{s_1}}=\frac{1}{{\pi }_{s_1}}P(O,s_1|{\lambda }^{(t)}))+\eta \text{\hspace{1em}(5-9)}$$

令偏导等于0，有

$$P(O,s_1|{\lambda }^{(t)}))+{\pi }_{s_1}^{(t+1)}\eta =0\text{\hspace{1em}(5-10)}$$

对所有的状态变量求和，有

$$\sum _{s_1=q_1}^{q_N}(P(O,s_1|{\lambda }^{(t)})+{\pi }_{s_1}^{(t+1)}\eta )=0$$

故有

$$P(O|{\lambda }^{(t)})+\eta =0$$

即
$$\eta =-P(O|{\lambda }^{(t)})\text{\hspace{1em}(5-11)}$$

将$(5-11)$代入到$(5-12)$有

$${\pi }_{s_1}^{(t+1)}=\frac{P(O,s_1|{\lambda }^{(t)})}{P(O|{\lambda }^{(t)})}\text{\hspace{1em}(5-12)}$$

终于求出来了，其他的参数用类似的方法求即可，不过会复杂一些。

6 Decoding

Decoding要解决的问题是

$$\hat{S}=arg\underset{S}{\max }P(S|O,\lambda )\text{\hspace{1em}(6-1)}$$

翻译过来就是，给定了模型参数，对应的状态变量序列最优可能是哪一组。

由于$P(O|\lambda )$是已经观测到的变量，我们也可以认为$(6-1)$等价于

$$\hat{S}=arg\underset{S}{\max }P(S|O,\lambda )P(O|\lambda )=arg\underset{S}{\max }P(S,O|\lambda )\text{\hspace{1em}(6-2)}$$

我们画个图看下

图5 decoding示意图

看图一下子就明白了，我们的每time step的状态变量都有$N$个状态，我们在每个time step选择一个状态变量，形成一条路径，使得经过整条路径的联合概率最大。

这里一共有$N^T$条路径，如果把每条路径的概率都算一遍，再取找概率最大的那条，时间复杂度就太高了。因此，我们用动态规划的思路去求解这个问题，也叫做Viterbi algorithm。

我们令

$${\delta }_t(q_i)=\underset{s_1,...,s_{t-1}}{\max }P(o_1,...,o_t,s_1,...,s_{t-1},s_t=q_i|\lambda )\text{\hspace{1em}(6-3)}$$

翻译一下就是，当$t$时刻的状态取$q_i$时，使得到$t$时刻为止的联合概率最大的状态路径$[s_1,...,s_{t-1}]$为${\delta }_t(q_i)$。

我们来看下${\delta }_{t+1}(q_j)$时刻和${\delta }_t(q_i)$的关系

$$\begin{array}{rl}{\delta }_{t+1}(q_j)& =\underset{s_1,...,s_t}{\max }P(o_1,...,o_{t+1},s_1,...,s_t,s_{t+1}=q_j|\lambda )\\ & 遍历t时刻所有的{\delta }_t(q_i)\\ & =\underset{1\le i\le N}{\max }{\delta }_t(q_i)a_{ij}{b}_j(o_{t+1})\end{array}\text{\hspace{1em}(6-4)}$$

这就是递推式了。用这个递推式，把所有的${\delta }_t(q_i)$算出来。

同时我们也要定义一个记录到每个time step的各个状态的最优的前一个状态

$${\psi }_t(j)=arg\underset{i}{\max }{\delta }_t(q_i)a_{ij}\text{\hspace{1em}(6-5)}$$

这个$(6-5)$是用来回溯路径的，这个学过动态规划的就再熟悉不过了。

最终我们会在最后一个time step找到是的概率最大的那个状态，记作

$$q_T^{\ast }=arg\underset{i}{\max }({\delta }_T(q_i))\text{\hspace{1em}(6-6)}$$

然后用$(6-5)$不断回溯就行了。

参考资料

[1] 机器学习-白板推导系列(十四)-隐马尔可夫模型HMM