概率图模型(1)--隐马尔科夫模型(2)

HMM模型参数求解概述

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

第一种情况较为简单，就是我们已知 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=\left[a_{ij}\right],其中a_{ij}=\frac{A_{ij}}{\sum _{s=1}^N{A}_{is}}$$

假设样本隐藏状态为 $q_i$ 且观测状态为 $v_k$ 的频率计数是 $B_{jk}$ ,那么观测状态概率矩阵为：

$$B=[b_j(k)],其中b_j(k)=\frac{B_{jk}}{\sum _{s=1}^M{B}_{js}}$$

假设所有样本中初始隐藏状态为 $q_i$ 的频率计数为 $C(i)$ ,那么初始概率分布为：

$$\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​)} 是已知的，此时我们能不能求出合适的HMM模型参数呢？

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

鲍姆-韦尔奇(Baum-Welch)算法原理

鲍姆-韦尔奇算法原理既然使用的就是EM算法的原理，那么我们需要在E步求出联合分布 $P(O,I|\lambda )$ 基于条件概率 $P(I|O,\overline{\lambda })$ 的期望，其中 $\overline{\lambda }$ 为当前的模型参数，然后再M步最大化这个期望，得到更新的模型参数 $\lambda$ 。接着不停的进行EM迭代，直到模型参数的值收敛为止。

首先来看看E步，当前模型参数为 $\overline{\lambda }$ , 联合分布 $P(O,I|\lambda )$ 基于条件概率 $P(I|O,\overline{\lambda })$ 的期望表达式为：

$$L(\lambda ,\overline{\lambda })=\sum _{I}P(I|O,\overline{\lambda })logP(O,I|\lambda )$$

在M步，我们极大化上式，然后得到更新后的模型参数如下：

$$\overline{\lambda }=arg\underset{\lambda }{\max }\sum _{I}P(I|O,\overline{\lambda })logP(O,I|\lambda )$$

通过不断的E步和M步的迭代，直到 $\overline{\lambda }$ 收敛。下面我们来看看鲍姆-韦尔奇算法的推导过程。

鲍姆-韦尔奇算法的推导

我们的训练数据为 { ( 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​)} ，其中任意一个观测序列 O d = { o 1 ( d ) , o 2 ( d ) , . . . o T ( d ) } O_d = \{o_1^{(d)}, o_2^{(d)}, ... o_T^{(d)}\} Od​={o1(d)​,o2(d)​,...oT(d)​} ,其对应的未知的隐藏状态序列表示为： I d = { i 1 ( d ) , i 2 ( d ) , . . . i T ( d ) } I_d = \{i_1^{(d)}, i_2^{(d)}, ... i_T^{(d)}\} Id​={i1(d)​,i2(d)​,...iT(d)​}

首先看鲍姆-韦尔奇算法的E步，我们需要先计算联合分布 $P(O,I|\lambda )$ 的表达式如下：

$$P(O,I|\lambda )=\prod _{d=1}^D{\pi }_{i_1^{(d)}}{b}_{i_1^{(d)}}(o_1^{(d)})a_{i_1^{(d)}{i}_2^{(d)}}{b}_{i_2^{(d)}}(o_2^{(d)})...a_{i_{T-1}^{(d)}{i}_T^{(d)}}{b}_{i_T^{(d)}}(o_T^{(d)})$$

我们的E步得到的期望表达式为：

$$L(\lambda ,\overline{\lambda })=\sum _{I}P(I|O,\overline{\lambda })logP(O,I|\lambda )$$

在M步我们要极大化上式。由于 $P(I|O,\overline{\lambda })=P(I,O|\overline{\lambda })/P(O|\overline{\lambda })$ ,而 $P(O|\overline{\lambda })$ 是常数，因此我们要极大化的式子等价于：

$$\overline{\lambda }=arg\underset{\lambda }{\max }\sum _{I}P(O,I|\overline{\lambda })logP(O,I|\lambda )$$

我们将上面 $P(O,I|\lambda )$ 的表达式带入我们的极大化式子，得到的表达式如下：

$$\overline{\lambda }=arg\underset{\lambda }{\max }\sum _{d=1}^D\sum _{I}P(O,I|\overline{\lambda })(log{\pi }_{i_1}+\sum _{t=1}^{T-1}log{a}_{i_t,i_{t+1}}+\sum _{t=1}^{T}log{b}_{i_t}(o_t))$$

我们的隐藏模型参数 $\lambda =(A,B,\pi )$ ,因此下面我们只需要对上式分别对 $A,B,\pi$ 求导即可得到我们更新的模型参数 $\overline{\lambda }$

首先我们看看对模型参数 $\pi$ 的求导。由于 $\pi$ 只在上式中括号里的第一部分出现，因此我们对于 $\pi$ 的极大化式子为：

$$\overline{{\pi }_i}=arg\underset{{\pi }_{i_1}}{\max }\sum _{d=1}^D\sum _{I}P(O,I|\overline{\lambda })log{\pi }_{i_1}=arg\underset{{\pi }_i}{\max }\sum _{d=1}^D\sum _{i=1}^{N}P(O,i_1^{(d)}=i|\overline{\lambda })log{\pi }_i$$

由于 ${\pi }_i$ 还满足 $\sum _{i=1}^N{\pi }_i=1$ ，因此根据拉格朗日子乘法，我们得到 ${\pi }_i$ 要极大化的拉格朗日函数为：

$$arg\underset{{\pi }_i}{\max }\sum _{d=1}^D\sum _{i=1}^{N}P(O,i_1^{(d)}=i|\overline{\lambda })log{\pi }_i+\gamma (\sum _{i=1}^N{\pi }_i-1)$$

其中， $\gamma$ 为拉格朗日系数。上式对 ${\pi }_i$ 求偏导数并令结果为0， 我们得到：

$$\sum _{d=1}^{D}P(O,i_1^{(d)}=i|\overline{\lambda })+\gamma {\pi }_i=0$$

令 i 分别等于从1到 N，从上式可以得到 N 个式子，对这 N 个式子求和可得：

$$\sum _{d=1}^{D}P(O|\overline{\lambda })+\gamma =0$$

从上两式消去 $\gamma$ ,得到 ${\pi }_i$ 的表达式为：

$${\pi }_i=\frac{\sum _{d=1}^{D}P(O,i_1^{(d)}=i|\overline{\lambda })}{\sum _{d=1}^{D}P(O|\overline{\lambda })}=\frac{\sum _{d=1}^{D}P(O,i_1^{(d)}=i|\overline{\lambda })}{DP(O|\overline{\lambda })}=\frac{\sum _{d=1}^{D}P(i_1^{(d)}=i|O,\overline{\lambda })}{D}=\frac{\sum _{d=1}^{D}P(i_1^{(d)}=i|O^{(d)},\overline{\lambda })}{D}$$

利用我们在前向概率的定义可得：

$$P(i_1^{(d)}=i|O^{(d)},\overline{\lambda })={\gamma }_1^{(d)}(i)$$

因此最终我们在M步 ${\pi }_i$ 的迭代公式为：

$${\pi }_i=\frac{\sum _{d=1}^D{\gamma }_1^{(d)}(i)}{D}$$

现在我们来看看 A 的迭代公式求法。方法和$\pi$的类似。由于 A 只在最大化函数式中括号里的第二部分出现，而这部分式子可以整理为：

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

由于 $a_{ij}$ 还满足 $\sum _{j=1}^N{a}_{ij}=1$ 。和求解 ${\pi }_i$ 类似，我们可以用拉格朗日子乘法并对 $a_{ij}$ 求导，并令结果为0，可以得到 $a_{ij}$ 的迭代表达式为：

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

利用前向概率的定义和 ${\xi }_t(i,j)$ 的定义可得们在M步 $a_{ij}$ 的迭代公式为：

$$a_{ij}=\frac{\sum _{d=1}^D\sum _{t=1}^{T-1}{\xi }_t^{(d)}(i,j)}{\sum _{d=1}^D\sum _{t=1}^{T-1}{\gamma }_t^{(d)}(i)}$$

现在我们来看看 B 的迭代公式求法。方法和$\pi$ 的类似。由于 B 只在最大化函数式中括号里的第三部分出现，而这部分式子可以整理为：

$$\sum _{d=1}^D\sum _I\sum _{t=1}^{T}P(O,I|\overline{\lambda })log{b}_{i_t}(o_t)=\sum _{d=1}^D\sum _{j=1}^N\sum _{t=1}^{T}P(O,i_t^{(d)}=j|\overline{\lambda })log{b}_j(o_t)$$

由于 $b_j(o_t)$ 还满足 $\sum _{k=1}^M{b}_j(o_t=v_k)=1$ 。和求解 ${\pi }_i$ 类似，我们可以用拉格朗日子乘法并对 $b_j(k)$ 求导，并令结果为0，得到 $b_j(k)$ 的迭代表达式为：

$$b_j(k)=\frac{\sum _{d=1}^D\sum _{t=1}^{T}P(O,i_t^{(d)}=j|\overline{\lambda })I(o_t^{(d)}=v_k)}{\sum _{d=1}^D\sum _{t=1}^{T}P(O,i_t^{(d)}=j|\overline{\lambda })}$$

其中 $I(o_t^{(d)}=v_k)$ 当且仅当 $o_t^{(d)}=v_k$ 时为1，否则为0. 利用前向概率的定义可得 $b_j(o_t)$ 的最终表达式为：

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

有了 ${\pi }_i,a_{ij},b_j(k)$ 的迭代公式，我们就可以迭代求解HMM模型参数了。

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

这里我们概括总结下鲍姆-韦尔奇算法的流程。

输入： D 个观测序列样本 { ( O 1 ) , ( O 2 ) , . . . ( O D ) } \{(O_1), (O_2), ...(O_D)\} {(O1​),(O2​),...(OD​)}

输出：HMM模型参数

1)随机初始化所有的 ${\pi }_i,a_{ij},b_j(k)$

2) 对于每个样本 $d=1,2,...D$，用前向后向算法计算${\gamma }_t^{(d)}(i)，{\xi }_t^{(d)}(i,j),t=1,2...T$

3) 更新模型参数：

$${\pi }_i=\frac{\sum _{d=1}^D{\gamma }_1^{(d)}(i)}{D}$$

$$a_{ij}=\frac{\sum _{d=1}^D\sum _{t=1}^{T-1}{\xi }_t^{(d)}(i,j)}{\sum _{d=1}^D\sum _{t=1}^{T-1}{\gamma }_t^{(d)}(i)}$$

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

4. 如果 ${\pi }_i,a_{ij},b_j(k)$ 的值已经收敛，则算法结束，否则回到第2）步继续迭代。

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(I^{\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∗​} 。要这样近似求解不难，利用定义：在给定模型λλ和观测序列 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解码的方法。

维特比算法

维特比算法是一个通用的解码算法，是基于动态规划的求序列最短路径的方法。

既然是动态规划算法，那么就需要找到合适的局部状态，以及局部状态的递推公式。在HMM中，维特比算法定义了两个局部状态用于递推。

第一个局部状态是在时刻 t 隐藏状态为i所有可能的状态转移路径 $i_1,i_2,...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)$ 的定义可以得到 $\delta$ 的递推表达式：

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

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

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

维特比算法流程总结

现在我们来总结下维特比算法的流程：

输入：HMM模型 $\lambda =(A,B,\pi )$，观测序列$O=(o_1,o_2,...o_T)$

输出：最有可能的隐藏状态序列 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)=\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]b_i(0_t),i=1,2...N$$

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

3. 计算时刻 T 最大的 ${\delta }_T(i)$ ,即为最可能隐藏状态序列出现的概率。计算时刻 T最大的 ${\Psi }_t(i)$ ,即为时刻 T 最可能的隐藏状态。

$$P\ast =\underset{1\le j\le N}{\max }{\delta }_T(i)$$

$$i_T^{\ast }=arg\underset{1\le j\le N}{\max }[{\delta }_T(i)]$$

4. 利用局部状态 ${\Psi }_t(i)$ 开始回溯。对于 $t=T-1,T-2,...,1$：

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

最终得到最有可能的隐藏状态序列 I ∗ = { i 1 ∗ , i 2 ∗ , . . . i T ∗ } I^*= \{i_1^*,i_2^*,...i_T^*\} I∗={i1∗​,i2∗​,...iT∗​}

HMM模型维特比算法总结

如果了解文本挖掘的分词中的维特比算法，就会发现这两篇个维特比算法稍有不同。主要原因是在中文分词时，我们没有观察状态和隐藏状态的区别，只有一种状态。但是维特比算法的核心是定义动态规划的局部状态与局部递推公式，这一点在中文分词维特比算法和HMM的维特比算法是相同的，也是维特比算法的精华所在。

维特比算法也是寻找序列最短路径的一个通用方法，和dijkstra算法有些类似，但是dijkstra算法并没有使用动态规划，而是贪心算法。同时维特比算法仅仅局限于求序列最短路径，而dijkstra算法是通用的求最短路径的方法。