隐马尔科夫模型的三个基本问题

前面简单介绍了隐马尔科夫模型，下面继续看看HMM的三个基本问题。

问题一、评估观察序列概率

问题描述：假设我们已知HMM模型的参数$\lambda =(A,B,\Pi )$，其中A是隐藏状态转移概率的矩阵，B是观测状态生成概率的矩阵， $\Pi$是隐藏状态的初始概率分布。同时给定一个观测序列 O = { o 1 , o 2 , . . . o T } O =\{o_1,o_2,...o_T\} O={o1​,o2​,...oT​}，现在我们要求这个观测序列?在模型?下出现的条件概率$P(O|\lambda )$。

1. 暴力求解法

我们可以列举出所有可能出现的长度为T的隐藏序列 I = { i 1 , i 2 , . . . , i T } I = \{i_1,i_2,...,i_T\} I={i1​,i2​,...,iT​}，分别求出这些隐藏序列与观测序列O的联合概率分布$P(O,I|\lambda )$，这样我们就可以很容易的求出边缘分布$P(O|\lambda )$了，也就是我们要求的观测序列$O$在模型$\lambda$下出现的条件概率$P(O|\lambda )$。

(1) 首先，任意一个隐藏序列 I = { i 1 , i 2 , . . . , i T } I = \{i_1,i_2,...,i_T\} I={i1​,i2​,...,iT​}出现的概率是$P(I|\lambda )$

(2) 对于固定的状态序列 $I$，我们要求的观察序列 $O$ 出现的概率是：$P(O|I,\lambda )$， 则$I$和$O$联合出现的概率是：$P(O,I|\lambda )=P(I|\lambda )P(O|I,\lambda )$

(3) 然后求边缘概率分布，即可得到观测序列 $O$ 在模型 $\lambda$下出现的条件概率：$P(O|\lambda )=\sum _{I}P(O,I|\lambda )$

虽然上述方法有效，但是如果我们的隐藏状态数?非常多的那就麻烦了，此时我们预测状态有$N^T$种组合，算法的时间复杂度是$O(T{N}^T)$阶的。因此对于一些隐藏状态数极少的模型，我们可以用暴力求解法来得到观测序列出现的概率，但是如果隐藏状态多，则上述算法太耗时，我们需要寻找其他简洁的算法。

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

2. 前向算法

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

在前向算法中，通过定义“前向概率”来定义动态规划的这个局部状态。什么是前向概率呢, 其实定义很简单：假定t时刻隐藏状态为$q_i$，且观测状态序列为$o_1,o_2,...o_t$ 的概率就是前向概率，并记为：${\alpha }_t(i)=P(o_1,o_2,...o_t,h_t=q_i|\lambda )$。

递推关系式如下：

$${\alpha }_{t+1}(j)=\sum _{i=1}^N{\alpha }_t(i)a_{ij}{b}_j(o_{t+1})$$

我们的动态规划从时刻1开始，到时刻T结束，由于${\alpha }_T(i)$表示在时刻 T 观测序列为$o_1,o_2,...o_T$并且时刻T隐藏状态$q_i$的概率。因此我们只要将所有隐藏状态对应的概率相加，即$\sum _{i=1}^N{\alpha }_T(i)$就得到了在时刻T观测序列为$o_1,o_2,...o_T$的概率，也就是我们要求的这个观测序列$O$在模型 $\lambda$ 下出现的条件概率$P(O|\lambda )$。

总结下前向算法的全部过程：

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

(1) 计算时刻1的各个隐藏状态前向概率：${\alpha }_1(i)={\pi }_i{b}_i(o_1),i=1,2,...N$，${\pi }_i$是时刻时刻1隐藏状态为i的概率，$b_i(o_1)$是观测状态为$o_1$的概率，${\alpha }_1(i)$是时刻1状态为i的前向概率。其实就是i取1,2,...,N，要计算所有状态的前向概率。

(2) 递推时刻2,3,…T时刻各个隐藏状态的前向概率：${\alpha }_{t+1}(i)=\sum _{i=1}^N{\alpha }_t(i)a_{ij}{b}_j(o_{t+1}),i=1,2,...N$。这里j也取1,2,...,N，也要计算所有状态的前向概率。

(3) 计算最终结果：$P(O|\lambda )=\sum _{i=1}^N{\alpha }_T(i)$，这里相加的是T时刻各隐藏状态的前向概率。

从递推公式可以看出，我们的算法时间复杂度是$O(T{N}^2)$，比暴力解法的时间复杂度$O(T{N}^T)$少了几个数量级。

3. 前向算法实例

以隐马尔科夫模型简介中的盒子与球的例子来显示前向概率的计算。

我们的观察集合是:V={红，白}，M=2，我们的状态集合是:Q ={盒子1，盒子2，盒子3}， N=3。而观察序列和状态序列的长度为3。

初始状态分布为：$\Pi =(0.2,0.4,0.4{)}^T$

状态转移概率分布矩阵为：$A=\left(\begin{array}{ccc}0.5& 0.2& 0.3\\ 0.3& 0.5& 0.2\\ 0.2& 0.3& 0.5\end{array}\right)$

观测状态概率矩阵为：$B=\left(\begin{array}{cc}0.5& 0.5\\ 0.4& 0.6\\ 0.7& 0.3\end{array}\right)$

球的颜色的观测序列:O={红，白，红}

按照我们上一节的前向算法。首先计算时刻1三个状态的前向概率：

(1) 时刻1是红色球:

隐藏状态是盒子1的概率为：${\alpha }_1(1)={\pi }_1{b}_1(o_1)=0.2\times 0.5=0.1$

隐藏状态是盒子2的概率为：${\alpha }_1(2)={\pi }_2{b}_2(o_1)=0.4\times 0.4=0.16$

隐藏状态是盒子3的概率为：${\alpha }_1(3)={\pi }_3{b}_3(o_1)=0.4\times 0.7=0.28$

(2) 现在我们可以开始递推了，首先递推时刻2三个状态的前向概率

时刻2是白色球:

隐藏状态是盒子1的概率为：${\alpha }_2(1)=[\sum _{i=1}^3{\alpha }_1(i)a_{i1}]b_1(o_2)=[0.1\ast 0.5+0.16\ast 0.3+0.28\ast 0.2]\times 0.5=0.077$

隐藏状态是盒子2的概率为：${\alpha }_2(2)=[\sum _{i=1}^3{\alpha }_1(i)a_{i2}]b_2(o_2)=[0.1\ast 0.2+0.16\ast 0.5+0.28\ast 0.3]\times 0.6=0.1104$

隐藏状态是盒子3的概率为：${\alpha }_2(3)=[\sum _{i=1}^3{\alpha }_1(i)a_{i3}]b_3(o_2)=[0.1\ast 0.3+0.16\ast 0.2+0.28\ast 0.5]\times 0.3=0.0606$

(3) 继续递推，现在我们递推时刻3三个状态的前向概率

时刻3是红色球:

隐藏状态是盒子1的概率为：${\alpha }_3(1)=[\sum _{i=1}^3{\alpha }_2(i)a_{i1}]b_1(o_3)=[0.077\ast 0.5+0.1104\ast 0.3+0.0606\ast 0.2]\times 0.5=0.04187$

隐藏状态是盒子2的概率为：${\alpha }_3(2)=[\sum _{i=1}^3{\alpha }_2(i)a_{i2}]b_2(o_3)=[0.077\ast 0.2+0.1104\ast 0.5+0.0606\ast 0.3]\times 0.4=0.03551$

隐藏状态是盒子3的概率为：${\alpha }_3(3)=[\sum _{i=1}^3{\alpha }_3(i)a_{i3}]b_3(o_3)=[0.077\ast 0.3+0.1104\ast 0.2+0.0606\ast 0.5]\times 0.7=0.05284$

(4) 最终我们求出观测序列:O={红，白，红}的概率为：

$$P(O|\lambda )=\sum _{i=1}^3{\alpha }_3(i)=0.13022$$

4. 后向算法

后向算法和前向算法非常类似，都是用的动态规划，唯一的区别是选择的局部状态不同，后向算法用的是“后向概率”，那么后向概率是如何定义的呢？

定义时刻t时隐藏状态为$q_i$, 从时刻t+1到最后时刻T的观测状态的序列为$o_{t+1},o_{t+2},...o_T$的概率为后向概率。记为：

$${\beta }_t(i)=P(o_{t+1},o_{t+2},...o_T|h_t=q_i,\lambda )$$

后向概率的动态规划递推公式和前向概率是相反的:

• 现在我们假设我们已经找到了在时刻t+1时各个隐藏状态的后向概率。${\beta }_{t+1}(j)$，现在我们需要递推出时刻t时各个隐藏状态的后向概率。
• 如下图，我们可以计算出观测状态的序列为$o_{t+2},o_{t+3},...o_T$， 且t时刻隐藏状态为$q_i$, t+1时刻隐藏状态为$q_j$的概率为$a_{ij}{\beta }_{t+1}(j)$。
• 接着可以得到观测状态的序列为$o_{t+1},o_{t+2},...o_T$， t时隐藏状态为$q_i$, 时刻t+1隐藏状态为$q_j$的概率为$a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)$, 则把下面所有线对应的概率加起来，我们可以得到观测状态的序列为$o_{t+1},o_{t+2},...o_T$， t时隐藏状态为$q_i$的概率为$\sum _{j=1}^N{a}_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)$，这个概率即为时刻t的后向概率。

这样我们得到了后向概率的递推关系式如下：

$${\beta }_t(i)=\sum _{j=1}^N{a}_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)$$

现在我们总结下后向算法的流程,注意下和前向算法的相同点和不同点：

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

(1) 初始化时刻T的各个隐藏状态后向概率：${\beta }_T(i)=1,i=1,2,...N$。注意这里所有状态的后向概率为1。

(2) 递推时刻T−1,T−2,…1时刻的后向概率：${\beta }_t(i)=\sum _{j=1}^N{a}_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j),i=1,2,...N$

(3) 计算最终结果：$P(O|\lambda )=\sum _{i=1}^N{\pi }_i{b}_i(o_1){\beta }_1(i)$
　　　　
此时我们的算法时间复杂度仍然是$O(T{N}^2)$。

5. HMM常用概率的计算

利用前向概率和后向概率，我们可以计算出HMM中单个状态和两个状态的概率公式。

(1) 给定模型$\lambda$和观测序列$O$,在时刻t处于状态$q_i$的概率记为:

$${\gamma }_t(i)=P(h_t=q_i|O,\lambda )=\frac{P(h_t=q_i,O|\lambda )}{P(O|\lambda )}$$

利用前向概率和后向概率的定义可知：$$P(h_t=q_i,O|\lambda )={\alpha }_t(i){\beta }_t(i)$$

于是我们得到：

$${\gamma }_t(i)=\frac{{\alpha }_t(i){\beta }_t(i)}{\sum _{j=1}^N{\alpha }_t(j){\beta }_t(j)}$$

(2) 给定模型$\lambda$和观测序列$O$，在时刻t处于状态$q_i$，且时刻t+1处于状态$q_j$的概率记为:

$${\xi }_t(i,j)=P(h_t=q_i,h_{t+1}=q_j|O,\lambda )=\frac{P(h_t=q_i,h_{t+1}=q_j,O|\lambda )}{P(O|\lambda )}$$

而$P(h_t=q_i,h_{t+1}=q_j,O|\lambda )$可以由前向后向概率来表示为:

$$P(h_t=q_i,h_{t+1}=q_j,O|\lambda )={\alpha }_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)$$

从而最终我们得到${\xi }_t(i,j)$的表达式如下：

$${\xi }_t(i,j)=\frac{{\alpha }_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)}{\sum _{r=1}^N\sum _{s=1}^N{\alpha }_t(r)a_{rs}{b}_s(o_{t+1}){\beta }_{t+1}(s)}$$

(3) 将${\gamma }_t(i)$和${\xi }_t(i,j)$在各个时刻t求和，可以得到：

• 在观测序列O下状态 i 出现的期望值: $\sum _{t=1}^T{\gamma }_t(i)$
• 在观测序列O下由状态 i 转移的期望值: $\sum _{t=1}^{T-1}{\gamma }_t(i)$
• 在观测序列O下由状态 i 转移到状态 j 的期望值: $\sum _{t=1}^{T-1}{\xi }_t(i,j)$

上面这些常用的概率值在求解HMM问题二，即求解HMM模型参数的时候需要用到。

问题二、最可能隐藏状态序列求解

问题描述：在HMM模型的解码问题中，给定模型$\lambda =(A,B,\Pi )$和观测序列 O = { o 1 , o 2 , . . . o T } O =\{o_1,o_2,...o_T\} O={o1​,o2​,...oT​}，求给定观测序列O条件下，最可能出现的对应的状态序列 I ∗ = { i 1 ∗ , i 2 ∗ , . . . i T ∗ } I^*= \{i_1^*,i_2^*,...i_T^*\} I∗={i1∗​,i2∗​,...iT∗​},即$P(T^{\ast }|O)$要最大化。

一个可能的近似解法是求出观测序列O在每个时刻t最可能的隐藏状态$i_t^{\ast }$然后得到一个近似的隐藏状态序列 I ∗ = { i 1 ∗ , i 2 ∗ , . . . i T ∗ } I^*= \{i_1^*,i_2^*,...i_T^*\} I∗={i1∗​,i2∗​,...iT∗​}。要这样近似求解不难，利用前一个问题的前向后向算法评估观察序列概率中第五节(HMM常用概率)的定义：在给定模型$\lambda$和观测序列O时，在时刻t处于状态$q_i$的概率是${\gamma }_t(i)$，这个概率可以通过HMM的前向算法与后向算法计算。这样我们有：

$$i_t^{\ast }=arg\underset{1\le i\le N}{\max }[{\gamma }_t(i)],t=1,2,...T$$

近似算法很简单，但是却不能保证预测的状态序列是整体是最可能的状态序列，因为预测的状态序列中某些相邻的隐藏状态可能存在转移概率为0的情况。

而维特比算法可以将HMM的状态序列作为一个整体来考虑，避免近似算法的问题，维特比算法是一个通用的求序列最短路径的动态规划算法，也可以用于很多其他问题。下面我们来看看维特比算法进行HMM解码的方法。

1. 维特比算法概述

维特比算法是一个通用的解码算法，是基于动态规划的求序列最短路径的方法。既然是动态规划算法，那么就需要找到合适的局部状态，以及局部状态的递推公式。在HMM中，维特比算法定义了两个局部状态用于递推。
　　　　
(1) 第一个局部状态是在时刻t隐藏状态为i所有可能的状态转移路径$i_1,i_2,...i_t$（$i\to i_t$）中的概率最大值。记为${\delta }_t(i)$:
$${\delta }_t(i)=\underset{i_1,i_2,...i_{t-1}}{\max }P(i_t=i,i_1,i_2,...i_{t-1},o_t,o_{t-1},...o_1|\lambda ),i=1,2,...N$$

由${\delta }_t(i)$的定义可以得到?的递推表达式：

$$\begin{array}{ll}{\delta }_{t+1}(i)& ={\max }_{i_1,i_2,...i_t}P(i_{t+1}=i,i_1,i_2,...i_t,o_{t+1},o_t,...o_1|\lambda )\\ & ={\max }_{1\le j\le N}[{\delta }_t(j)a_{ji}]b_i(o_{t+1})\end{array}$$

(2)第二个局部状态由第一个局部状态递推得到。我们定义在时刻t隐藏状态为i的所有单个状态转移路径$(i_1,i_2,...,i_{t-1},i)$中概率最大的转移路径中第t−1个时刻的隐藏状态为$i_{t-1}={\Psi }_t(i_t)$,其递推表达式可以表示为：

$${\Psi }_t(i_t)=arg\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]$$

有了这两个局部状态，我们就可以先从时刻0一直递推到时刻T，然后再利用${\Psi }_t(i_t)$记录的前一个最可能的状态节点回溯，直到找到最优的隐藏状态序列。

2. 维特比算法的流程

输入：HMM模型$\lambda =(A,B,\Pi )$，观测序列 O = { o 1 , o 2 , . . . o T } O =\{o_1,o_2,...o_T\} O={o1​,o2​,...oT​}
输出：最有可能的隐藏状态序列 I ∗ = { i 1 ∗ , i 2 ∗ , . . . i T ∗ } I^*= \{i_1^*,i_2^*,...i_T^*\} I∗={i1∗​,i2∗​,...iT∗​}

(1) 初始化局部状态：

${\delta }_1(i)={\pi }_i{b}_i(o_1),i=1,2...N$
${\Psi }_1(i)=0,i=1,2...N$

(2) 进行动态规划递推时刻t=2,3,…,T时刻的局部状态：

${\delta }_t(i)={\max }_{1\le j\le N}[{\delta }_{t-1}(j)a_{ji}]b_i(0_t),i=1,2...N$
${\Psi }_t(i)=arg{\max }_{1\le j\le N}[{\delta }_{t-1}(j)a_{ji}],i=1,2...N$

(3) 如果t=1,2,…,T所有时刻的两个局部状态${\delta }_t(i)$和${\Psi }_t(i)$都已经计算完成

那就要从时刻T开始，计算该时刻最大的${\delta }_T(i)$,即为最可能隐藏状态序列出现的概率。计算时刻T最大的${\Psi }_t(i)$,即为时刻T最可能的隐藏状态$i_T^{\ast }$。

${\delta }_T(i)$中的最大概率：$P^{\ast }={\max }_{1\le j\le N}{\delta }_T(i)$
对应了T时刻最可能的隐藏状态：$i_T^{\ast }=arg{\max }_{1\le j\le N}[{\delta }_T(i)]$

(4) 利用局部状态$\Psi (i)$开始回溯。对于t=T−1,T−2,…,1：

$i_t^{\ast }={\Psi }_{t+1}(i_{t+1}^{\ast })$

也就是将t+1时刻的隐藏状态$i_{t+1}^{\ast }$代入到上面这个等式，从而求得t时刻的最可能隐藏状态$i_t^{\ast }$，最终得到最有可能的隐藏状态序列 I ∗ = { i 1 ∗ , i 2 ∗ , . . . i T ∗ } I^*= \{i_1^*,i_2^*,...i_T^*\} I∗={i1∗​,i2∗​,...iT∗​}。

3. HMM维特比算法求解实例

仍使用上面的盒子与球的例子来看看HMM维特比算法求解。

按照我们上一节的维特比算法:

(1) 首先需要得到三个隐藏状态在时刻1时对应的各自两个局部状态${\delta }_1(i)$和${\Psi }_1(i)$，此时观测状态为1（红球）：

${\delta }_1(1)={\pi }_1{b}_1(o_1)=0.2\times 0.5=0.1$
${\delta }_1(2)={\pi }_2{b}_2(o_1)=0.4\times 0.4=0.16$
${\delta }_1(3)={\pi }_3{b}_3(o_1)=0.4\times 0.7=0.28$
${\Psi }_1(1)={\Psi }_1(2)={\Psi }_1(3)=0$

(2) 现在开始递推三个隐藏状态在时刻2时对应的各自两个局部状态，此时观测状态为2（白球）：

${\delta }_2(1)={\max }_{1\le j\le 3}[{\delta }_1(j)a_{j1}]b_1(o_2)={\max }_{1\le j\le 3}[0.1\times 0.5,0.16\times 0.3,0.28\times 0.2]\times 0.5=0.028$
${\Psi }_2(1)=arg{\max }_{1\le j\le 3}[{\delta }_1(j)a_{j1}]=arg{\max }_{1\le j\le 3}[0.1\times 0.5,0.16\times 0.3,0.28\times 0.2]=3$

${\delta }_2(2)={\max }_{1\le j\le 3}[{\delta }_1(j)a_{j2}]b_2(o_2)={\max }_{1\le j\le 3}[0.1\times 0.2,0.16\times 0.5,0.28\times 0.3]\times 0.6=0.0504$
${\Psi }_2(3)=arg{\max }_{1\le j\le 3}[{\delta }_1(j)a_{j2}]=arg{\max }_{1\le j\le 3}[0.1\times 0.2,0.16\times 0.5,0.28\times 0.3]=3$

${\delta }_2(3)={\max }_{1\le j\le 3}[{\delta }_1(j)a_{j3}]b_3(o_2)={\max }_{1\le j\le 3}[0.1\times 0.3,0.16\times 0.2,0.28\times 0.5]\times 0.3=0.042$
${\Psi }_2(3)=arg{\max }_{1\le j\le 3}[{\delta }_1(j)a_{j3}]=arg{\max }_{1\le j\le 3}[0.1\times 0.3,0.16\times 0.2,0.28\times 0.5]=3$

(3) 继续递推三个隐藏状态在时刻3时对应的各自两个局部状态，此时观测状态为1（红球）：

${\delta }_3(1)={\max }_{1\le j\le 3}[{\delta }_2(j)a_{j1}]b_1(o_3)={\max }_{1\le j\le 3}[0.028\times 0.5,0.0504\times 0.3,0.042\times 0.2]\times 0.5=0.00756$
${\Psi }_3(1)=2$

${\delta }_3(2)={\max }_{1\le j\le 3}[{\delta }_2(j)a_{j2}]b_2(o_3)={\max }_{1\le j\le 3}[0.028\times 0.2,0.0504\times 0.5,0.042\times 0.3]\times 0.4=0.01008$
${\Psi }_3(2)=2$

${\delta }_3(3)={\max }_{1\le j\le 3}[{\delta }_2(j)a_{j3}]b_3(o_3)={\max }_{1\le j\le 3}[0.028\times 0.3,0.0504\times 0.2,0.042\times 0.5]\times 0.7=0.0147$
${\Psi }_3(3)=3$

(4) 此时已经到最后的时刻，我们开始准备回溯。
此时最大概率为$P^{\ast }={\delta }_3(3)$,从而得到$i_3^{\ast }=3$，也就是时刻3最可能的隐藏状态是盒子3。将$i_3^{\ast }=3$代入$i_2^{\ast }={\Psi }_3(i_3^{\ast })$，所以$i_2^{\ast }=3$。同样将$i_2^{\ast }=3$代入$i_1^{\ast }={\Psi }_2(i_2^{\ast })$，所以$i_1^{\ast }=3$。从而得到最终的最可能的隐藏状态序列为：(3,3,3)。

问题三、模型参数学习问题

HMM模型参数求解根据已知的条件可以分为两种情况。

1. 第一种情况

这种情况较为简单，就是我们已知D个长度为T的观测序列和对应的隐藏状态序列，即 { ( O 1 , I 1 ) , ( O 2 , I 2 ) , . . . ( O D , I D ) } \{(O_1, I_1), (O_2, I_2), ...(O_D, I_D)\} {(O1​,I1​),(O2​,I2​),...(OD​,ID​)}是已知的，此时我们可以很容易的用最大似然来求解模型参数。

• 假设样本从隐藏状态$q_i$转移到$q_j$的频率计数是$A_{ij}$,那么状态转移矩阵A求得为：
  $A=\left[a_{ij}\right],其中a_{ij}=\frac{A_{ij}}{\sum _{s=1}^N{A}_{is}}$

• 假设样本隐藏状态为$q_j$且观测状态为$v_k$的频率计数是$B_{jk}$,那么观测状态概率矩阵B为：
  $B=[b_j(k)],其中b_j(k)=\frac{B_{jk}}{\sum _{s=1}^M{B}_{js}}$

• 假设所有样本中初始隐藏状态为$q_i$的频率计数为C(i),那么初始概率分布Π为：
  $\Pi =\pi (i)=\frac{C(i)}{\sum _{s=1}^{N}C(s)}$

可见第一种情况下求解模型还是很简单的，其实就是基于统计计算模型参数。但是在很多时候，我们无法得到HMM样本观察序列对应的隐藏序列，只有D个长度为T的观测序列，即 { ( O 1 ) , ( O 2 ) , . . . ( O D ) } \{(O_1), (O_2), ...(O_D)\} {(O1​),(O2​),...(OD​)}是已知的，那就属于第二种情况了。

2. 第二种情况

也是我们本文要讨论的重点。它的解法最常用的是鲍姆-韦尔奇（Baum-Welch）算法，其实就是基于EM算法的求解，只不过鲍姆-韦尔奇算法出现的时代，EM算法还没有被抽象出来，所以我们本文还是说鲍姆-韦尔奇算法。

假设给定训练数据中包含S个长度为T的观测序列$O=(o_1,o_2,...,o_T)$，而没有对应的隐藏状态序列，目标是学习隐马尔可夫模型$\lambda =(A,B,\pi )$的参数。我们将观测序列数据看作观测数据$O$，状态序列数据看作不可观测的隐数据$I$，那么隐马尔可夫模型事实上是一个含有隐变量的概率模型 :
$$P(O|\lambda )=\sum _{I}P(O|I,\lambda )P(I|\lambda )\text{\hspace{1em}(3.2.1)}$$

(1) 确定完全数据的对数似然函数
所有观测数据写成$O=(o_1,o_2,...,o_T)$，所有隐数据写成$I=(i_1,i_2,...,i_T)$，完全数据是$(O,I)=(o_1,o_2,...,o_T,i_1,i_2,...,i_T)$。完全数据的对数似然函数是$\log P(O,I|\lambda )$（似然函数还是条件概率的形式，但是目的是求参数）。

(2) EM算法的E步：写出Q函数

首先Q函数的定义：$Q(\lambda ,\overline{\lambda })=E_I[logP(O,I|\lambda )|O,\overline{\lambda }]$。

Q函数是完全数据的对数似然函数关于给定模型参数和观测变量的前提下对隐变量的条件概率分布的期望。按照Q函数的定义在当前问题的具体形式如下：
$$Q(\lambda ,\overline{\lambda })=\sum _{I}logP(O,I|\lambda )P(I|O,\overline{\lambda })\text{\hspace{1em}(3.2.2)}$$

这里仅仅取$P(O,I|\lambda )$的对数，$P(I|O,\overline{\lambda })$是在对数的外面。其中，$\overline{\lambda }$是隐马尔可夫模型参数的当前估计值，$\lambda$是要极大化的隐马尔可夫模型参数。

并且我们注意到Q函数的后部分，根据条件联合分布公式可得 ：$P(I|O,\overline{\lambda })=\frac{P(O,I|\overline{\lambda })}{P(O|\overline{\lambda })}$。而这个$1/P(O|\overline{\lambda })$其实是一个常数因子，可以省略，因此得到$Q(\lambda ,\overline{\lambda })$为：
$$Q(\lambda ,\overline{\lambda })=\sum _I\log P(O,I|\lambda )P(O,I|\overline{\lambda })\text{\hspace{1em}(3.2.3)}$$

又因为$P(O,I|\lambda )={\pi }_{i_1}{b}_{i_1}(o_1)a_{i_1{i}_2}{b}_{i_2}(o_2)\cdots a_{i_{T-1}{i}_T}{b}_{i_T}(o_T)$，所以

$$\begin{array}{ll}Q(\lambda ,\overline{\lambda })=& {\sum }_{I}log{\pi }_{i_1}P(O,I|,\overline{\lambda })\\ & +{\sum }_I\left({\sum }_{t=1}^{T-1}log{a}_{i_t{i}_{t+1}}\right)P(O,I|,\overline{\lambda })\\ & +{\sum }_I\left({\sum }_{t=1}^{T}log{b}_{i_t}(o_t)\right)P(O,I|,\overline{\lambda })\end{array}\text{\hspace{1em}(3.2.4)}$$

此时我们看到待估计的参数刚好分别出现在三个项中，所以只需对各个项分别极大化。

(3) EM算法的M步：极大化Q函数（公式3.2.4）

第一项：${\sum }_{I}log{\pi }_{i_1}P(O,I|,\overline{\lambda })$

=${\sum }_I\log {\pi }_{i_1}P(O,I|\overline{\lambda })={\sum }_{j=1}^N\log {\pi }_{j}P(O,i_1=j|\overline{\lambda })$

注意到${\pi }_j$满足约束条件${\sum }_{j=1}^N{\pi }_j=1$，利用拉格朗日乘子法，写出拉格朗日函数：
$f_1={\sum }_{j=1}^N\log {\pi }_{j}P(O,i_1=j|\overline{\lambda })+\theta ({\sum }_{j=1}^N{\pi }_j-1)$

对其求偏导数并令结果为0：
$\frac{\partial f_1}{{\partial }_{{\pi }_j}}=\frac{\partial }{{\partial }_{{\pi }_j}}[{\sum }_{j=1}^N\log {\pi }_{j}P(O,i_1=j|\overline{\lambda })+\theta ({\sum }_{j=1}^N{\pi }_j-1)]=0$

注意这里的${\pi }_j$其实是个向量，所以后半部分$\frac{\partial }{\partial {\pi }_j}\gamma ({\sum }_{j=1}^N{\pi }_j-1)=\theta {\pi }_j$，于是得到：
$P(O,i_1=j|\overline{\lambda })+\theta {\pi }_j=0$

对j求和得到拉格朗日因子：$\theta =-P(O|\overline{\lambda })$

于是得：
${\pi }_j=\frac{P(O,i_1=j|\overline{\lambda })}{P(O|\overline{\lambda })}$

根据上面常用概率的关于$\gamma$的定义$:{\gamma }_t(i)=P(h_t=q_i|O,\lambda )=\frac{P(h_t=q_i,O|\lambda )}{P(O|\lambda )}$
可得初始状态概率为：
$${\pi }_j={\gamma }_1(j)\text{\hspace{1em}(3.2.5)}$$

第二项：${\sum }_I\left({\sum }_{t=1}^{T-1}log{a}_{i_t{i}_{t+1}}\right)P(O,I|,\overline{\lambda })$

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

类似第一项，应用具有约束条件${\sum }_{j=1}^N{a}_{ij}=1$的拉格朗日乘子法可以求出：
$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 })}$

同样根据${\gamma }_t(j)$定义以及${\xi }_t(i,j)$的定义${\xi }_t(i,j)=P(h_t=q_i,h_{t+1}=q_j|O,\lambda )=\frac{P(h_t=q_i,h_{t+1}=q_j,O|\lambda )}{P(O|\lambda )}$可得状态转移概率为：
$$a_{ij}=\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{{\sum }_{t=1}^T{\gamma }_t(j)}\text{\hspace{1em}(3.2.6)}$$

第三项：${\sum }_I\left({\sum }_{t=1}^{T}log{b}_{i_t}(o_t)\right)P(O,I|,\overline{\lambda })$

=${\sum }_{j=1}^N{\sum }_{t=1}^T\log b_j(o_t)P(O,i_t=j|\overline{\lambda })$

同样用拉格朗日乘子法，约束条件是${\sum }_{k=1}^M{b}_j(k)=1$。注意，只有在$o_t=v_k$时$b_j(o_t)$对$b_j(k)$的偏导数才不为0，以$I(o_t=v_k)$表示，求得：
$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 })}$

同样根据${\gamma }_t(j)$定义可得观测状态概率为：
$$b_j(k)=\frac{{\sum }_{t=1,o_t=v_k}^T{\gamma }_t(j)}{{\sum }_{t=1}^T{\gamma }_t(j)}\text{\hspace{1em}(3.2.7)}$$

3. 鲍姆-韦尔奇算法流程总结

输入：观测数据$O=(o_1,o_2,...,o_T)$
输出：隐马尔可夫模型参数

(1) 初始化
对n=0，选取$a_{ij}^{(0)},b_j(k{)}^{(0)},{\pi }_i^{(0)}$，得到模型${\lambda }^{(0)}=(A^{(0)},B^{(0)},{\pi }^{(0)})$

(2) 递推，对n=1,2,…
$$a_{ij}^{(n+1)}=\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{{\sum }_{t=1}^{T-1}{\gamma }_t(i)}$$

$$b_j^{(n+1)}(k)=\frac{{\sum }_{t=1,o_t=v_k}^T{\gamma }_t(j)}{{\sum }_{t=1}^T{\gamma }_t(j)}$$

$${\pi }_i^{(n+1)}={\gamma }_1(i)$$

(3) 模型收敛（满足条件）则终止递推，得到模型参数${\lambda }^{(n+1)}=(A^{(n+1)},B^{(n+1)},{\pi }^{(n+1)})$

---

参考：

• 隐马尔可夫模型之Baum-Welch算法详解
• HMM的Baum-Welch算法和Viterbi算法公式推导细节
• 隐马尔科夫模型HMM（三）鲍姆-韦尔奇算法求解HMM参数
• 《统计学习方法（李航）- 第10章 隐马尔科可夫模型》

---

THE END.