隐马尔可夫模型（HMM）推导详解

1、隐马尔可夫模型基本概念

隐马尔可夫（Hidden Markov Model, HMM）是马尔可夫链(Markov chain)的扩展，同样适合于对有限状态演进的随机过程建模。在马尔可夫链中，状态是可观察得到的，然而在某些情况下，状态隐藏在观测事件下，状态与观测变量不再满足一一对应的关系，此时观测变量是状态的概率分布函数。举例来说，医生可以通过患者的行为（发热、打喷嚏等）判断病人的身体状况，但无法不通过观察直接得到。此时，发热、打喷嚏等行为即观测事件，身体状况即病人此时的状态。HMM就是用来描述此类的马尔可夫链。

2、隐马尔可夫模型基本参数

HMM各个时期会处于各种状态，设$Q$是所有可能状态的集合；每个状态可以产生各种观测事件，设$V$是所有可能的观测事件的集合。$Q$，$V$的定义如下：
Q = { q 1 , q 2 , . . . q N } , V = { v 1 , v 2 , . . . v M } Q=\{q_1, q_2, ... q_N\}, V=\{v_1, v_2, ... v_M\} Q={q1​,q2​,...qN​},V={v1​,v2​,...vM​}
其中$N$是可能的状态数，$M$是可能的观测数。还是以上面的医生看病为例，这里的$Q$表示{感冒，健康}，而$V$表示{发热、打喷嚏、无异常}。

然后定义$S$是长度为$T$的模型状态序列，定义$O$为对应的观测序列。
S = { s 1 , s 2 , . . . s T } , O = { o 1 , o 2 , . . . o T } S=\{s_1, s_2, ... s_T\}, O=\{o_1, o_2, ... o_T\} S={s1​,s2​,...sT​},O={o1​,o2​,...oT​}
这两个序列都是整数列表，其中每个整数对应的是状态集合$Q$和观测事件集合$V$中的索引。比如：$q_{s_2}$表示模型在时刻2时的状态，$v_{o_5}$表示模型在时刻5时的观测。

除此之外，还有三个分布被定义来描述一个隐马尔可夫模型：

• $\pi$：初始概率分布，表示在初始时刻各个状态出现的概率：
  $$\begin{gathered} \pi =[{\pi }_1,{\pi }_2,...,{\pi }_N{]}^T \\ {\pi }_i=P(s_1=i),i=1,2,...,N \end{gathered}$$
• $A$：状态转移分布，表示上一时刻到下一时刻不同状态之间转换的概率：
  $$\begin{gathered} A=[a_{ij}{]}_{N\times N} \\ {a}_{ij}=P(s_{t+1}=j|s_t=i) \end{gathered}$$
• $B$：观测概率分布，表示某个状态下各种观测事件出现的概率：
  $$\begin{gathered} B=[b_{ij}{]}_{N\times M} \\ {b}_{ij}=P(o_t=j|s_t=i) \end{gathered}$$

隐马尔可夫模型由以上三个分布决定，因此可以用一个三元符号表示：
$$\lambda =(A,B,\pi )$$

3、隐马尔可夫模型基本假设

HMM 是表示随机过程 { s 1 , s 2 , . . . , s T , o 1 , o 2 , . . . , o T } \{s_1, s_2, ..., s_T, o_1, o_2, ..., o_T\} {s1​,s2​,...,sT​,o1​,o2​,...,oT​}的联合概率分布。HMM 基于两条假设使得问题得到简化：

• 齐次马尔可夫性假设。任意时刻的状态只依赖于前一时刻的状态，与其它时刻的状态及观测无关：
  $$P(s_t|s_{t-1},o_{t-1},...,s_1,o_1)=P(s_t|s_{t-1})$$
• 观测独立性假设。任意时刻的观测只依赖于此刻的状态，与其它无关：
  $$P(o_t|s_t,o_{t-1},...,s_1,o_1)=P(o_t|s_t)$$

4、隐马尔可夫模型基本问题

经典HMM主要解决三个基本问题。

4.1 评估问题

给定模型$\lambda$和观测序列$O$，计算在该模型下观测序列$O$出现的概率$P(O|\lambda )$。

4.2 学习问题

已知观测序列$O$，估计模型$\lambda$的参数，使得在该模型下观测到这个序列的概率$P(O|\lambda )$最大。

4.3 预测问题

已知模型$\lambda$，给定观测序列$O$，求与之对应的状态序列$S$，使得概率$P(S|O,\lambda )$最大。

5、隐马尔可夫模型基本问题解法

5.1 评估问题

对于给定的模型参数$\lambda =(A,B,\pi )$，计算观测序列$O=(o_1,o_2,\dots ,o_T)$出现的概率，最简单的就是把所有可能的状态序列的概率都计算出来，然后选择最大的那个，也就是俗称的穷举法。但是这个方法计算复杂度是极大的，高达$O(T{N}^T)$阶，所以不可行。为解决这一问题，两种计算可行的算法被提出。

5.1.1 前向算法

首先对前向概率进行定义：
$${\alpha }_t(i)=P(o_1,o_2,...,o_t,s_t=i|\lambda )$$
也就是说，${\alpha }_t(i)$表示在时刻$t$时，观测序列为$(o_1,o_2,\dots ,o_t)$，且此时状态为$q_i$的概率。前向算法基本流程如下：

• 定义初值。即$t=1$时在观察到$o_1$的情况下所有状态发生的概率，一共有$N$个（因为有$N$种可能的状态）：
  $${\alpha }_1(i)=P(o_1,s_1=i|\lambda )={\pi }_i{b}_{i{o}_1},i=1,2,...,N$$

• 递推。对$t=2,3,...,T$，计算所有状态发生的概率，每个时间点下同样都要计算N个概率：
  $$\begin{array}{rl}{\alpha }_t(i)& =P(o_1,o_2,...,o_t,s_t=i|\lambda )\\ & =\sum _{j=1}^{N}P(o_1,o_2,...,o_t,s_{t-1}=j,s_t=i|\lambda )\\ & =\sum _{j=1}^{N}P(o_t|o_1,o_2,...,o_{t-1},s_{t-1}=j,s_t=i,\lambda )P(o_1,o_2,...,o_{t-1},s_{t-1}=j,s_t=i|\lambda )\\ & =\sum _{j=1}^{N}P(o_t|s_t=i,\lambda )P(s_t=i|o_1,o_2,...,o_{t-1},s_{t-1}=j,\lambda )P(o_1,o_2,...,o_{t-1},s_{t-1}=j|\lambda )\\ & =P(o_t|s_t=i,\lambda )\sum _{j=1}^{N}P(s_t=i|s_{t-1}=j,\lambda )P(o_1,o_2,...,o_{t-1},s_{t-1}=j|\lambda )\\ & =b_{i{o}_t}\sum _{j=1}^N{a}_{ji}{\alpha }_{t-1}(j),i=1,2,...,N\end{array}$$
  其中：

  • 第一行表示${\alpha }_t(i)$的定义；
  • 第二行表示考虑了$t-1$时刻所有状态的概率；
  • 第三行中
    • 第二项是在该模型参数下取得观察序列 { o 1 , o 2 , . . . , o t − 1 } \{o_1, o_2, ..., o_{t-1}\} {o1​,o2​,...,ot−1​}以及第$t-1$时刻状态为$j$和第$t$时刻状态为$i$的概率
    • 第一项是在已知观察序列 { o 1 , o 2 , . . . , o t − 1 } \{o_1, o_2, ..., o_{t-1}\} {o1​,o2​,...,ot−1​}以及第$t-1$时刻状态为$j$和第$t$时刻状态为$i$的情况下取得$t$时刻观察值为$o_t$的概率；
  • 第四行中
    • 第三项是在该模型参数下取得观察序列 { o 1 , o 2 , . . . , o t − 1 } \{o_1, o_2, ..., o_{t-1}\} {o1​,o2​,...,ot−1​}以及第$t-1$时刻状态为$j$的概率
    • 第二项是在已知观察序列 { o 1 , o 2 , . . . , o t − 1 } \{o_1, o_2, ..., o_{t-1}\} {o1​,o2​,...,ot−1​}以及第$t-1$时刻状态为$j$的情况下取得$t$时刻状态为$i$的概率
    • 第一项是在已知$t$时刻状态为$i$的情况下取得$t$时刻观察值为$o_t$的概率；
  • 第五行中的操作是将第四行中与变量$j$无关的项提到求和操作外面。

• 计算最后的概率。
  $$P(O|\lambda )=\sum _{i=1}^{N}P(o_1,o_2,...,o_T,s_T=i|\lambda )=\sum _{j=1}^N{\alpha }_T(i)$$

5.1.2 后向算法

首先对后向概率进行定义：
$${\beta }_t(i)=P(o_{t+1},o_{t+2},...,o_T|s_t=i,\lambda )$$
也就是说，${\beta }_t(i)$表示已知时刻$t$时状态为$q_i$的条件下，之后观测序列为$(o_{t+1},o_{t+2},...,o_T)$的概率。后向算法基本流程如下：

• 定义初值。由于是从后往前推，所以要定义${\beta }_T(i)$的值：
  $${\beta }_T(i)=1,i=1,2,...,N$$

• 递推。对$t=T-1,T-2,...,1$，计算：
  $$\begin{array}{rl}{\beta }_t(i)& =P(o_{t+1},o_{t+2},...,o_T|s_t=i,\lambda )\\ & =\sum _{j=1}^{N}P(o_{t+1},o_{t+2},...,o_T,s_{t+1}=j|s_t=i,\lambda )\\ & =\sum _{j=1}^{N}P(o_{t+1}|o_{t+2},...,o_T,s_{t+1}=j,s_t=i,\lambda )P(o_{t+2},...,o_T,s_{t+1}=j|s_t=i,\lambda )\\ & =\sum _{j=1}^{N}P(o_{t+1}|s_{t+1}=j,\lambda )P(o_{t+2},...,o_T|s_{t+1}=j,s_t=i,\lambda )P(s_{t+1}=j|s_t=i,\lambda )\\ & =\sum _{j=1}^{N}P(o_{t+1}|s_{t+1}=j,\lambda )P(o_{t+2},...,o_T|s_{t+1}=j,\lambda )P(s_{t+1}=j|s_t=i,\lambda )\\ & =\sum _{j=1}^N{a}_{ij}{b}_{j{o}_{t+1}}{\beta }_{t+1}(j),i=1,2,...,N\end{array}$$
  其中：

  • 第一行表示${\beta }_t(i)$的定义；
  • 第二行表示考虑了$t+1$时刻所有状态的概率；
  • 第三行中
    • 第二项是在该模型参数下已知第$t$时刻状态为$i$的情况下取得观察序列 { o t + 2 , . . . , o T } \{o_{t+2}, ..., o_T\} {ot+2​,...,oT​}以及第$t+1$时刻状态为$j$的概率
    • 第一项是在已知第$t$时刻状态为$i$、第$t+1$时刻状态为$j$以及观察序列 { o t + 2 , . . . , o T } \{o_{t+2}, ..., o_T\} {ot+2​,...,oT​}的情况下取得$t+1$时刻观察值为$o_{t+1}$的概率；
  • 第四行中
    • 第三项是在已知第$t$时刻状态为$i$的情况下取得第$t+1$时刻状态为$j$的概率
    • 第二项是在已知第$t$时刻状态为$i$以及第$t+1$时刻状态为$j$的情况下取得观察序列 { o t + 2 , . . . , o T } \{o_{t+2}, ..., o_T\} {ot+2​,...,oT​}的概率
    • 第一项是根据观测独立性假设对第三行第一项的简化表示；
  • 第五行相对于第四行只有第二项有改变，这一项将$s_t=i$这个条件去掉了，这是因为，根据观测独立性假设，第$t$时刻的状态对第$t+2$及以后时刻的观察值不产生影响。

• 计算最后的概率。
  $$P(O|\lambda )=\sum _{i=1}^N{\pi }_i{b}_{i{o}_1}{\beta }_1(i)$$

5.2 学习问题

在学习问题中，只有观测序列是已知的，我们需要根据这些观测值来估计概率模型的参数$\lambda$。另外，一般会有不止一个观测序列被用于该训练过程，即会进行多次实验得到多个观测序列。这种问题常使用Baum-Welch算法（实际上就是EM算法的应用）。

Q函数

因为只有观测序列，其未知的状态序列就可以看做隐变量。在EM算法推导详解一文中，我们知道求解带有隐函数的概率模型参数需要使用EM算法，而此算法中最重要的一步就是构造Q函数，然后通过Q函数不断更新模型参数$\lambda$。

假设用于训练的观测序列集合为：
O = { O 1 , O 2 , . . . , O ∣ O ∣ } \Bbb{O}=\{O_1, O_2, ..., O_{|\Bbb{O}|}\} O={O1​,O2​,...,O∣O∣​}
其中任意观测序列$O_i$的长度都为$T$，共有$|\mathbb{O}|$个观测序列被用于模型训练过程。

假定隐变量状态序列的所有可能性集合为：
S = { S = ( s 1 , s 2 , . . . , s T ) ∣ s t = 1 , 2 , . . . , N , t = 1 , 2 , . . . , T } \Bbb{S}=\{S=(s_1, s_2, ..., s_T)|s_t=1, 2, ..., N, t=1, 2, ..., T\} S={S=(s1​,s2​,...,sT​)∣st​=1,2,...,N,t=1,2,...,T}
因此，根据EM算法推导详解一文中推导出的Q函数，我们可以定义第$k$次迭代的Q函数：
$$Q(\lambda ,{\lambda }^k)=\sum _{\mathbb{O}}\sum _{\mathbb{S}}P(S|O,{\lambda }^k)logP(S,O|\lambda )$$
代入待优化参数，然后进行分解：
$$\begin{array}{rl}Q(\lambda ,{\lambda }^k)& =\sum _{\mathbb{O}}\sum _{\mathbb{S}}P(S|O,{\lambda }^k)log({\pi }_{s_1}{b}_{s_1{o}_1}{a}_{s_1{s}_2}{b}_{s_2{o}_2}\cdot \cdot \cdot a_{s_{T-1}{s}_T}{b}_{s_T{o}_T})\\ & =\sum _{\mathbb{O}}\sum _{\mathbb{S}}P(S|O,{\lambda }^k)(log{\pi }_{s_1}+\sum _{t=1}^{T-1}log{a}_{s_t{s}_{t+1}}+\sum _{t=1}^{T}log{b}_{s_t{o}_t})\end{array}$$
可以看到待优化式子分解为三个不同待优化变量的求和，因此可以将它们分开分别进行优化。

迭代

使用EM算法迭代计算HMM参数的步骤如下：

• 1、初始化${\lambda }^0=({\pi }^0,A^0,B^0)$；
• 2、第$k$次迭代，在${\lambda }^{k-1}$的基础上对$\pi$，$a$和$b$进行优化，得到${\lambda }^k$；
• 3、判断终止条件，满足即完成迭代，否则回到第二步执行下一次迭代。

5.3 预测问题

预测问题最常用的方法是维特比算法，它实际上就是动态规划，在递推的过程中保持步步最优，从而最终达到全局最优的目的。

定义在$t$时刻以内，且$t$时刻状态为$q_i$的所有状态序列的出现概率的最大值为：
$${\delta }_t(i)=\underset{s_1,...,s_{t-1}}{\max }P(s_t=i,s_{t-1},...,s_1,o_t,...,o_1|\lambda )$$
由定义可以获得$\delta$的递推公式：
$$\begin{array}{rl}{\delta }_{t+1}(i)& =\underset{s_1,...,s_t}{\max }P(s_{t+1}=i,s_t,...,s_1,o_{t+1},...,o_1|\lambda )\\ & =\underset{s_t}{\max }\underset{s_1,...,s_{t-1}}{\max }P(s_{t+1}=i,s_t,...,s_1,o_{t+1},...,o_1|\lambda )\\ & =\underset{1\le j\le N}{\max }\underset{s_1,...,s_{t-1}}{\max }P(s_{t+1}=i,s_t=j,s_{t-1},...,s_1,o_{t+1},...,o_1|\lambda )\\ & =\underset{1\le j\le N}{\max }\underset{s_1,...,s_{t-1}}{\max }[P(o_{t+1}|s_{t+1}=i,\lambda )P(s_{t+1}=i,s_t=j,s_{t-1},...,s_1,o_t,...,o_1|\lambda )]\\ & =P(o_{t+1}|s_{t+1}=i,\lambda )\underset{1\le j\le N}{\max }\underset{s_1,...,s_{t-1}}{\max }[P(s_{t+1}=i|s_t=j,\lambda )P(s_t=j,s_{t-1},...,s_1,o_t,...,o_1|\lambda )]\\ & =P(o_{t+1}|s_{t+1}=i,\lambda )\underset{1\le j\le N}{\max }[P(s_{t+1}=i|s_t=j,\lambda )\underset{s_1,...,s_{t-1}}{\max }P(s_t=j,s_{t-1},...,s_1,o_t,...,o_1|\lambda )]\\ & =b_{i{o}_{t+1}}\underset{1\le j\le N}{\max }{a}_{ji}{\delta }_t(j)\end{array}$$
以上算的是概率，并没有获取每个时刻下的状态编号。所以再定义${\Psi }_t(i)$函数来回溯获取状态序列：
$${\Psi }_t(i)=\underset{1\le j\le N}{\mathrm{argmax}}[{\delta }_{t-1}(j)a_{ji}],i=1,2,...,N$$
维特比算法用已知观测序列预测状态序列过程如下：

• 1、初始化。直接计算${\delta }_1(i)$；
• 2、对$t=2,3,...,T$，使用递推公式计算${\delta }_t(i)$；
• 3、使用${\Psi }_t(i)$函数回溯，获取最大概率状态序列。

参考资料

[1] 隐马尔可夫模型——概念与算法详解
[2] EM 算法推导详解