HMM Algorithm Theory

0. 简述

    本文主要参考"A Tutorial On Hidden Markov Models"，这个文章可以在http://vision.ai.uiuc.edu/dugad/上面找到。主要是通过阅读这篇技术报告来了解一些HMM的学习算法，即分段K-Means和向前向后算法（也称为BaumWelch算法）。

    前面两篇文章介绍过HMM的基础了，这里就更多的直接上公式了。那两篇文章是：隐马尔可夫模型-HMM-简述-1-原理-示例，隐马尔可夫模型-HMM-简述-2-评估-解码-学习。

1. 关于符号说明   

    这里简单说明一下，这篇文章将隐藏状态简称为状态，将观察状态简称为符号。

    N：状态的种类，M：符号的种类，T：符号序列长度（这里假设观察样本长度都是有一样长的）。

    i_t：t时刻的状态。

    {V}：符号集合。

    Π：状态的初始概率向量。

    a_i_j：P(t+1时刻处于状态j | t时刻处于状态i)，这是状态之间的转移，是相邻时间的，t时刻到t+1时刻。

    b_j(k)：P(t时刻的符号是k | t时刻处于状态j)，这是状态到符号的条件概率，是同一时间的，都是t时刻。

    O_t：t时刻的符号。

    λ=(A,B,Π)：HMM模型，或者说是HMM的参数。

    原文中的符号说明：    

2. 评估

    评估问题：已知HMM参数，计算符号序列的概率，即P(O | λ)。

    首先介绍一种直接的计算方式，计算该符号序列对应每一种状态序列的概率，然后将这些概率值求和，即为这个符号序列的概率。 

    I = i_1,i_2,...,i_T表示任意一个状态序列，O = O_1,O_2,...,O_T表示任意一个符号序列。

    P(O,I | λ) = P(O | λ,I) * P(I | λ)

    P(O | λ,I) = b_i_1(O_1) * b_i_2(O_2) * ··· * b_i_T(O_T)

    P(I | λ) = Π_i_1 * a_i_1_i_2 * ··· * a_i_T-1_i_T

    P(O,I | λ)对所有的I求和，即为P(O | λ)

    这样计算复杂度较高，对于每个状态I，计算时间为2T，状态I序列的数量为N^T，所以时间复杂度为2T * N^T，是T的指数函数。原文中的公式推导：

    下面给出向前算法，来解决这个问题。首先定义向前因子，α_t(i)=P(O_1,O_2,...,O_t, i_t=i | λ)，就是P(O, I|λ)从时刻1到时刻t的概率。这个α_t(i)是可以动态规划计算的。

    α_1(i) =  Π_i * b_i(O_1), 1<=i<=N

    α_t+1(j) =  Σ_i (α_t(i) * a_i_j) * b_j(O_t+1)

    P(O|λ) = Σ_i α_T(i)

  这里公式与前面公式中下标的区别值得注意，前面公式遍历每一种状态序列，i_1表示状态序列i的第一个状态，向前算法以及后面的向后算法，只是单纯针对状态来说的，因此1表示第一个状态。

    这样求α_1到α_T，一共T次，每次计算需要2*N^2次乘法，复杂度为2T * N^2。

    原文的公式推导：

    

    下面给出向后算法，来解决这个问题。定义向后因子，β_t(i) = P (O_t+1,O_t+2, ...O_T | i_t=i, λ)。递推公式如下：

    β_T(i) = 1, 1<=i<=N。

    β_t(i) = Σ_j (a_i_j) * b_j(O_t+1) * β_t+1(j)，从状态i转移到状态j，然后从状态j确定符号O_t+1，然后就是递归

    P(O|λ) = Σ_i  Π_i * b_i(O_1) * β_1(i) 

    向后算法的复杂度与向前算法一样，也是2T*N^2。

    

3. 解码

    解码问题：已知HMM参数，求产生某个符号序列概率最大的状态序列。即计算概率最大的那个P(O,I | λ)。这个问题实际上，可以使用前面的直接方法，向前算法呢和向后算法分别解决。直接方法就是算每个P(O,I | λ)，然后取最大值。向前算法和向后算法就是将公式中的求和变成求最大值即可。

    向前算法的变化：   

    α_1(i) =  Π_i * b_i(O_1), 1<=i<=N

    α_t+1(j) =  max{ (α_t(i) * a_i_j)  } * b_j(O_t+1), 1<=i<=N

    max{ P(O, I|λ) } = max{ α_T(i) }, 1<=i<=N

    最大状态序列= {i_1, i_2, ..., i_T}， 其中α_t(i_t) = max{a_t(j)}, 1<=j<=N

    向后算法类似，这里忽略了。解码常用的算法是维特比算法，即Viterbi Alogrithm，其实就是向前算法的变化，把概率变为-ln(概率)，这样通过乘法求最大值，变成了通过加法求最小值，计算量大大减小了。

    原文公式推导：   

    

4. 学习

    学习问题：根据符号序列和与其相关的状态集合，计算是符号序列概率最大的HMM参数。

    有两个方法，一个是分段KMeans（The Segmental K-means Algorithm），另一个是向前向后算法（也称为Baum-Welch算法）。

    The Segmental K-means Algorithm：通过最大化P(O,I | λ)，来调整λ=(A,B,Π)，其中I根据解码方法得到的最优状态序列。

    The Baum-Welch re-estimation Formulas：通过增加P(O | λ)到一个最大值，来调整λ=(A,B,Π)。P(O | λ)的计算需要针对涉及所有的状态序列I，通过将P(O,I | λ)求和得到。因此这个算法关注的是所有的状态序列，而segmental K-means关注的是最优的状态序列。

4.1 分割K-means算法 (The Segmental K-means Algorithm)

    假设λ为当前参数，λ'为新计算的参数，I为λ对应计算出的最优状态序列，I'为λ'计算出的最优状态序列，如果P(O,I|λ) < P(O,I'|λ')，那么继续迭代，令λ=λ'。训练这个模型，我们需要一些符号序列。假设有w个符号，每个符号序列为O=O1,O2,...,OT，长度都为T。每个符号假设是一个D(>=1)的向量。这个算法由以下几个步骤组成：

   第一步，随机选择N个符号(D维的向量)，然后将wT个符号的每一个符号，分到与其欧氏距离最小的类别中。这样产生了N个类，每个类称为一个状态(从1到N)。这个初始类的划分不会决定最终HMM的参数，但是会决定得到最终HMM参数需要迭代的次数。在我们的例子中，我们均匀的选取N个符号序列，使得这N个符号序列之间相间的符号序列个数相同，并且从每个符号序列中选取一个符号，即在第1个序列中选择第1个符号，在2个序列中选择第2个符号，...，在第N个符号序列中选取第N个符号。具体如何选择可以根据个人意愿进行选择。

    第二步，计算初始概率和转移概率，估计Π_i和a_i_j。Π_i=在类别i中O1的个数 / O1的个数。这个很好理解，比如，初始有w个符号序列，那么就有w个O1，分母就是w，分子就是这w个O1有多少个被分到类别i中。a_i_j=O_t在类别i中且O_t+1在类别j中的个数 / O_t在类别i中的个数。

     

    第三步，计算每个状态（指每个类别）的均值和方差。   

    

    第四步，估计观察状态到符号的概率。通过计算符号在每个状态（分类）中的分布来估计。这里假定符号的分布服从高斯分布。
    

    第五步，根据估计的HMM参数，计算每个观察序列的最优状态序列I'。如果I'与I不同，那么对应不同的Oi会被重新分配，比如I'与I在第3个元素上不同，那么O3就会重新分配到I'第三个元素对应的分类上去。

    第六步，如果在第二步中有重新分配的情况，那么跳转到第二步，进行迭代。

    分段K-means算法，在一些观察密度函数和我们假设的高斯密度函数上，收敛到“状态最优似然函数”。

4.2 向前向后算法（Baum-Welch 算法）

  这个方法首先假设有一个初始的HMM模型，然后通过下文介绍的方法来改进HMM模型使得P(O|λ)最大化。一个初始的HMM模型可以使用任何的方法构造，但是我们可以使用分段K-means算法的前五步来构造，这样可以得到一个相对合理的HMM的估计。不管怎样，我们经假设一个初始的HMM是已知的。这个算法通过调整参数λ，来使得P(O | λ)最大化。这个优化标准被称为最大似然标标准。P(O | λ)被称为似然函数。

    下面首先说明一些计算变量：
     γ_t(i)表示t时刻处于状态i的概率，分子刚好是向前因子与向后因子的乘积。
    

   ξ_t(i,j)表示t时刻在状态i，t+1时刻在状态j的概率。

    参数估计的公式为：   

    

    其中，Π_i的公式感觉应该是Π_i= γ_1(i)，即第一个时刻是状态i的概率。（不知道作者是不是写错了，感觉是写错了）
    a_i_j的分子是从i状态转移到j状态的联合概率，分母是i状态的概率，相除刚好是从i状态到j状态的转移概率。
    b_j(k)的分子是在j状态且符号为k的联合概率，分母是在j状态的概率，相除刚好是已知j状态，输出符号为k的条件概率。

    这样计算分两个步骤：第一，根据λ，计算λ'；第二，如果P(O | λ')>P(O | λ)，那么更新λ为λ'迭代计算，否则停止计算。

    Baum-Welch算法需要很多乘法，因此数值可能容易溢出，这需要一些可能的处理技术。而Segmental K-means算法使用Viterbi算法来产生I，因此溢出问题不严重。Segmental K-means的优势在于，有时我们比较关注最优状态序列产生的观察序列，而对所有状态序列产生的观察序列并不关心。此外Segmental K-means的计算量很小。