HMM学习最佳范例七：前向-后向算法

前向-后向算法(Forward-backward algorithm)

根据观察序列生成隐马尔科夫模型(Generating a HMM from a sequence of obersvations)

　　与HMM模型相关的“有用”的问题是评估（前向算法）和解码（维特比算法）——它们一个被用来测量一个模型的相对适用性，另一个被用来推测模型隐藏的部分在做什么（“到底发生了”什么）。可以看出它们都依赖于隐马尔科夫模型（HMM）参数这一先验知识——状态转移矩阵，混淆（观察）矩阵，以及pi向量（初始化概率向量）。
　　然而，在许多实际问题的情况下这些参数都不能直接计算的，而要需要进行估计——这就是隐马尔科夫模型中的学习问题。前向-后向算法就可以以一个观察序列为基础来进行这样的估计，而这个观察序列来自于一个给定的集合，它所代表的是一个隐马尔科夫模型中的一个已知的隐藏集合。
　　一个例子可能是一个庞大的语音处理数据库，其底层的语音可能由一个马尔可夫过程基于已知的音素建模的，而其可以观察的部分可能由可识别的状态（可能通过一些矢量数据表示）建模的，但是没有（直接）的方式来获取隐马尔科夫模型（HMM）参数。
　　前向-后向算法并非特别难以理解，但自然地比前向算法和维特比算法更复杂。由于这个原因，这里就不详细讲解前向-后向算法了（任何有关HMM模型的参考文献都会提供这方面的资料的）。
　　总之，前向-后向算法首先对于隐马尔科夫模型的参数进行一个初始的估计（这很可能是完全错误的），然后通过对于给定的数据评估这些参数的的价值并减少它们所引起的错误来重新修订这些HMM参数。从这个意义上讲，它是以一种梯度下降的形式寻找一种错误测度的最小值。
　　之所以称其为前向-后向算法，主要是因为对于网格中的每一个状态，它既计算到达此状态的“前向”概率（给定当前模型的近似估计），又计算生成此模型最终状态的“后向”概率（给定当前模型的近似估计）。 这些都可以通过利用递归进行有利地计算，就像我们已经看到的。可以通过利用近似的HMM模型参数来提高这些中间概率进行调整，而这些调整又形成了前向-后向算法迭代的基础。

注：关于前向-后向算法，原文只讲了这么多，后继我将按自己的理解补充一些内容。

要理解前向-后向算法，首先需要了解两个算法：后向算法和EM算法。后向算法是必须的，因为前向-后向算法就是利用了前向算法与后向算法中的变量因子，其得名也因于此；而EM算法不是必须的，不过由于前向-后向算法是EM算法的一个特例，因此了解一下EM算法也是有好处的，说实话，对于EM算法，我也是云里雾里的。好了，废话少说，我们先谈谈后向算法。

1、后向算法（Backward algorithm）
　　其实如果理解了前向算法，后向算法也是比较好理解的，这里首先重新定义一下前向算法中的局部概率at(i)，称其为前向变量，这也是为前向-后向算法做点准备：
　　　
　　相似地，我们也可以定义一个后向变量Bt(i)（同样可以理解为一个局部概率）：
　　　
　　后向变量(局部概率)表示的是已知隐马尔科夫模型及t时刻位于隐藏状态Si这一事实，从t+1时刻到终止时刻的局部观察序列的概率。同样与前向算法相似，我们可以从后向前（故称之为后向算法）递归地计算后向变量：
　　1）初始化，令t=T时刻所有状态的后向变量为1：
　　　　　
　　2）归纳，递归计算每个时间点，t=T-1,T-2,…,1时的后向变量：
　　
　　这样就可以计算每个时间点上所有的隐藏状态所对应的后向变量，如果需要利用后向算法计算观察序列的概率，只需将t=1时刻的后向变量（局部概率）相加即可。下图显示的是t+1时刻与t时刻的后向变量之间的关系：
　　　
　　上述主要参考自HMM经典论文《A tutorial on Hidden Markov Models and selected applications in speech recognition》。下面我们给出利用后向算法计算观察序列概率的程序示例，这个程序仍然来自于UMDHMM。

后向算法程序示例如下（在backward.c中):

void Backward(HMM *phmm, int T, int *O, double **beta, double *pprob)
{
　　int i, j; /* state indices */
　　int t; /* time index */
　　double sum;

　　/* 1. Initialization */
　　for (i = 1; i <= phmm->N; i++)
　　　　beta[T][i] = 1.0;

　　/* 2. Induction */
　　for (t = T - 1; t >= 1; t--)
　　{
　　　　for (i = 1; i <= phmm->N; i++)
　　　　{
　　　　　　sum = 0.0;
　　　　　　for (j = 1; j <= phmm->N; j++)
　　　　　　　　sum += phmm->A[i][j] *
　　　　　　　　　　　　　　(phmm->B[j][O[t+1]])*beta[t+1][j];
　　　　　　beta[t][i] = sum;
　　　　}
　　}

　　/* 3. Termination */
　　*pprob = 0.0;
　　for (i = 1; i <= phmm->N; i++)
　　　　*pprob += beta[1][i];
}

　　好了，后向算法就到此为止了，下一节我们粗略的谈谈EM算法。

前向-后向算法是Baum于1972年提出来的，又称之为Baum-Welch算法，虽然它是EM(Expectation-Maximization)算法的一个特例，但EM算法却是于1977年提出的。那么为什么说前向-后向算法是EM算法的一个特例呢？这里有两点需要说明一下。
　　第一，1977年A. P. Dempster、N. M. Laird、 D. B. Rubin在其论文“Maximum Likelihood from Incomplete Data via the EM Algorithm”中首次提出了EM算法的概念，但是他们也在论文的介绍中提到了在此之前就有一些学者利用了EM算法的思想解决了一些特殊问题，其中就包括了Baum在70年代初期的相关工作，只是这类方法没有被总结而已，他们的工作就是对这类解决问题的方法在更高的层次上定义了一个完整的EM算法框架。
　　第二，对于前向-后向算法与EM算法的关系，此后在许多与HMM或EM相关的论文里都被提及，其中贾里尼克（Jelinek）老先生在1997所著的书“Statistical Methods for Speech Recognition”中对于前向-后向算法与EM算法的关系进行了完整的描述，读者有兴趣的话可以找来这本书读读。
　　关于EM算法的讲解，网上有很多，这里我就不献丑了，直接拿目前搜索“EM算法”在Google排名第一的文章做了参考，希望读者不要拍砖：

　　EM 算法是 Dempster，Laind，Rubin 于 1977 年提出的求参数极大似然估计的一种方法，它可以从非完整数据集中对参数进行 MLE 估计，是一种非常简单实用的学习算法。这种方法可以广泛地应用于处理缺损数据，截尾数据，带有讨厌数据等所谓的不完全数据(incomplete data)。
　　假定集合Z = (X,Y)由观测数据 X 和未观测数据Y 组成，Z = (X,Y)和 X 分别称为完整数据和不完整数据。假设Z的联合概率密度被参数化地定义为P(X，Y|Θ)，其中Θ 表示要被估计的参数。Θ 的最大似然估计是求不完整数据的对数似然函数L(X;Θ)的最大值而得到的：
　　　L(Θ; X )= log p(X |Θ) = ∫log p(X ,Y |Θ)dY ；(1)
　　EM算法包括两个步骤：由E步和M步组成，它是通过迭代地最大化完整数据的对数似然函数Lc( X;Θ )的期望来最大化不完整数据的对数似然函数，其中：
　　　Lc(X;Θ) =log p(X，Y |Θ) ； (2)
　　假设在算法第t次迭代后Θ 获得的估计记为Θ(t ) ，则在（t+1）次迭代时，
　　E-步：计算完整数据的对数似然函数的期望，记为：
　　　Q(Θ |Θ (t) ) = E{Lc(Θ;Z)|X;Θ(t) }； (3)
　　M-步：通过最大化Q(Θ |Θ(t) ) 来获得新的Θ 。
　　通过交替使用这两个步骤，EM算法逐步改进模型的参数，使参数和训练样本的似然概率逐渐增大，最后终止于一个极大点。
　　直观地理解EM算法，它也可被看作为一个逐次逼近算法：事先并不知道模型的参数，可以随机的选择一套参数或者事先粗略地给定某个初始参数λ0 ，确定出对应于这组参数的最可能的状态，计算每个训练样本的可能结果的概率，在当前的状态下再由样本对参数修正，重新估计参数λ ，并在新的参数下重新确定模型的状态，这样，通过多次的迭代，循环直至某个收敛条件满足为止，就可以使得模型的参数逐渐逼近真实参数。
　　EM算法的主要目的是提供一个简单的迭代算法计算后验密度函数，它的最大优点是简单和稳定，但容易陷入局部最优。
　　参考原文见：http://49805085.spaces.live.com/Blog/cns!145C40DDDB4C6E5!670.entry

　　注意上面那段粗体字，读者如果觉得EM算法不好理解的话，就记住这段粗体字的意思吧！
　　有了后向算法和EM算法的预备知识，下一节我们就正式的谈一谈前向-后向算法。

隐马尔科夫模型（HMM）的三个基本问题中，第三个HMM参数学习的问题是最难的，因为对于给定的观察序列O，没有任何一种方法可以精确地找到一组最优的隐马尔科夫模型参数（A、B、pi）使P(O| )最大。因而，学者们退而求其次，不能使P(O| )全局最优，就寻求使其局部最优（最大化）的解决方法，而前向-后向算法（又称之为Baum-Welch算法）就成了隐马尔科夫模型学习问题的一种替代（近似）解决方法。
　　我们首先定义两个变量。 给定观察序列O及隐马尔科夫模型，定义t时刻位于隐藏状态Si的概率变量为：
　　　　　　　　
　　回顾一下 第二节中关于前向变量at(i)及后向变量Bt(i)的定义，我们可以很容易地将上式用前向、后向变量表示为：
　　　
　　其中分母的作用是确保：
　　 给定观察序列O及隐马尔科夫模型，定义t时刻位于隐藏状态Si及t+1时刻位于隐藏状态Sj的概率变量为：
　　　　
　　该变量在网格中所代表的关系如下图所示：
　
　　同样，该变量也可以由前向、后向变量表示：
　　　
　　而上述定义的两个变量间也存在着如下关系：
　　　　　　　　　　　　
　　如果对于时间轴t上的所有 相加，我们可以得到一个总和，它可以被解释为从其他隐藏状态访问Si的期望值（网格中的所有时间的期望），或者，如果我们求和时不包括时间轴上的t=T时刻，那么它可以被解释为从隐藏状态Si出发的状态转移期望值。相似地，如果对 在时间轴t上求和（从t=1到t=T-1），那么该和可以被解释为从状态Si到状态Sj的状态转移期望值。即：
　　　
　　　
上一节我们定义了两个变量及相应的期望值，本节我们利用这两个变量及其期望值来重新估计隐马尔科夫模型（HMM）的参数pi，A及B：

　　如果我们定义当前的HMM模型为，那么可以利用该模型计算上面三个式子的右端；我们再定义重新估计的HMM模型为，那么上面三个式子的左端就是重估的HMM模型参数。Baum及他的同事在70年代证明了因此如果我们迭代地的计算上面三个式子，由此不断地重新估计HMM的参数，那么在多次迭代后可以得到的HMM模型的一个最大似然估计。不过需要注意的是，前向-后向算法所得的这个结果（最大似然估计）是一个局部最优解。
　　关于前向-后向算法和EM算法的具体关系的解释，大家可以参考HMM经典论文《A tutorial on Hidden Markov Models and selected applications in speech recognition》，这里就不详述了。下面我们给出UMDHMM中的前向-后向算法示例，这个算法比较复杂，这里只取主函数，其中所引用的函数大家如果感兴趣的话可以自行参考UMDHMM。

前向-后向算法程序示例如下（在baum.c中):

void BaumWelch(HMM *phmm, int T, int *O, double **alpha, double **beta, double **gamma, int *pniter, double *plogprobinit, double *plogprobfinal)
{
　　int i, j, k;
　　int t, l = 0;

　　double logprobf, logprobb, threshold;
　　double numeratorA, denominatorA;
　　double numeratorB, denominatorB;

　　double ***xi, *scale;
　　double delta, deltaprev, logprobprev;

　　deltaprev = 10e-70;

　　xi = AllocXi(T, phmm->N);
　　scale = dvector(1, T);

　　ForwardWithScale(phmm, T, O, alpha, scale, &logprobf);
　　*plogprobinit = logprobf; /* log P(O |intial model) */
　　BackwardWithScale(phmm, T, O, beta, scale, &logprobb);
　　ComputeGamma(phmm, T, alpha, beta, gamma);
　　ComputeXi(phmm, T, O, alpha, beta, xi);
　　logprobprev = logprobf;

　　do
　　{

　　　　/* reestimate frequency of state i in time t=1 */
　　　　for (i = 1; i <= phmm->N; i++)
　　　　　　phmm->pi[i] = .001 + .999*gamma[1][i];

　　　　/* reestimate transition matrix and symbol prob in
　　　　　　　　each state */
　　　　for (i = 1; i <= phmm->N; i++)
　　　　{
　　　　　　denominatorA = 0.0;
　　　　　　for (t = 1; t <= T - 1; t++) 　　　　　　　　denominatorA += gamma[t][i];　　　　　　for (j = 1; j <= phmm->N; j++)
　　　　　　{
　　　　　　　　numeratorA = 0.0;
　　　　　　　　for (t = 1; t <= T - 1; t++) 　　　　　　　　　　numeratorA += xi[t][i][j];　　　　　　　　phmm->A[i][j] = .001 +
　　　　　　　　　　　　　　　　　.999*numeratorA/denominatorA;
　　　　　　}

　　　　　　denominatorB = denominatorA + gamma[T][i];
　　　　　　for (k = 1; k <= phmm->M; k++)
　　　　　　{
　　　　　　　　numeratorB = 0.0;
　　　　　　　　for (t = 1; t <= T; t++) 　　　　　　　　{　　　　　　　　　　if (O[t] == k) 　　　　　　　　　　　　numeratorB += gamma[t][i];　　　　　　　　}　　　　　　　　phmm->B[i][k] = .001 +
　　　　　　　　　　　　　　　　　.999*numeratorB/denominatorB;
　　　　　　}
　　　　}

　　　　ForwardWithScale(phmm, T, O, alpha, scale, &logprobf);
　　　　BackwardWithScale(phmm, T, O, beta, scale, &logprobb);
　　　　ComputeGamma(phmm, T, alpha, beta, gamma);
　　　　ComputeXi(phmm, T, O, alpha, beta, xi);

　　　　/* compute difference between log probability of
　　　　　　two iterations */
　　　　delta = logprobf - logprobprev;
　　　　logprobprev = logprobf;
　　　　l++;

　　}
　　while (delta > DELTA); /* if log probability does not
　　　　　　　　　　　　　　change much, exit */

　　*pniter = l;
　　*plogprobfinal = logprobf; /* log P(O|estimated model) */
　　FreeXi(xi, T, phmm->N);
　　free_dvector(scale, 1, T);
}

　　
　　前向-后向算法就到此为止了。

未完待续：总结

注：原创文章，转载请注明出处“我爱自然语言处理”：www.52nlp.cn

本文链接地址：http://www.52nlp.cn/hmm-learn-best-practices-seven-forward-backward-algorithm-5