隐马尔科夫链模型.docx

后验解码问题
前面我们介绍了在给定一个符号序列时，在不同的可能状态序列中找出概率（可能性）最大的状态序列（或路径）。在实际情况中可能会碰到一个问题：不同的路径其概率相差不大，这为我们选择哪一条路径增加了一定的困难。但我们知道一条状态序列，它是由各个位置的状态组成的，比如图6-11中，从第1个位置到第45个位置均由状态F组成，紧接后面的21个是L状态，依此类推。因此克服上述困难我们可以每个位置各个不同状态的概率，这里面我们必须分清一点：某一条序列对应于某个状态序列的概率含义与状态序列中某个位置上具体某个状态的概率是不同的，我们将它们的区分总结于如图6-14：

图6-14某个DNA序列对应的不同可能的CpG岛分布（状态序列）的概率与各个位置上某个状态的概率之间的区别（“+”代表CpG岛区域；“-”代表非CpG岛区域）。
我们在上节的“Viterbi”算法中重点介绍的是如何计算图6-14中右边部分，并将其中最大的概率对应的路径找出来。而我们这一节要介绍的是如何计算图6-14中最下面的概率，这就是通常我们所说的“后验解码问题”：计算概率为多少？
为了计算这个概率，首先得介绍所给定序列对应的概率，由于不同的状态序列可以产生同一条序列，因此，我们有图6-15的描述：

图6-15隐马尔克夫链中符号序列的概率与路径的关系（这里假设符号序列集只有两条序列与，对应的状态序列（路径）中也只有两个路径与。
图6-15中的与就是要求的某个符号序列的概率，由于此时路径与已知，因此它与符号序列的联合概率表达式可写为：
（6-11）
公式（6-11）罗列了一大堆符号，这对学生物学的人来说要比较快地接受它们可能会要多花点时间，为使我们有一个清晰的直观印象，我们这里以两条DNA序列及相应的CpG岛分布图来说明。

图6-16DNA序列中CpG岛分布的隐马尔科夫链模型中在公式（6-6）中的示意图
如果我们以代表中的某个序列，相应的状态序列有，则有：
(6-12)
在实际情况中，可能的路径（状态序列）很多，而且随着符号序列的长度的增加，姿态态序列的个数N呈指数函数的方式往上增长，以CpG岛为例：当序列长度为3时，其可能的状态序列个数为，当序列长度为4时，则有。因此如何计算（6-11）则需要一定的技巧。借鉴Viterbi算法，我们也不妨来个沿着符号序列逐步计算，这种逐步的方式有两种：一种是从符号序列的左端向右端（DNA序列的5’端到3’端，蛋白质序列的N端到C端），相应的算法叫前向算法（forwardalgorithm）；另一种是从右端向左端，相应的算法为后向算法(backwardalgorithm)。我们首先介绍前向算法。
一、前向算法(ForwardAlgorithm)
在介绍前向算法之前，我们先看一下式（6-7）的第一个公式即：
（6-13）
从中我们可发现它取的是每个位置最大值概率，因为对符号序列，其每个位置的概率等其各个状态概率之和，因此，整个符号序列的概率是每个位置概率之积，沿着符号序列的方向，我们有：
(6-14)
整个算法可写成如下：
初始化：(6-15a)
迭代过程：(6-15b)
终止：(6-15c)
我们可将前向算法以图示方式表示：

图6-17前向算法的基本框架图

图6-18前向算法的上体计算图。这里我们仅选状态1的计算过程，目的是使图形显得简洁明了，因为其它状的算法与此是相同的。

计算例子：符号序列“3151”。
我们仍以Viterbi算法中的例子即符号序列“3151”来说明：
计算第一个点：
初始条件如Viterbi算法所示，同时设定。
首先，我们计算第一个位置即，我们有：

当时，，我们有：
当，，我们有：
当时，，我们有：
在终止阶段，我们有：

与Viterbi算法相似，前向算法也存着“因计算机浮点数的限制导致后面的结果失真”，因此同样的需要用对数运算来代替相应的乘法运算，有文献曾试图通过如下的转换：

克服这一点。但我们认为这样做其实没有实质性的意义，因为中间一个“exp”的运算，如果log(p)小到一定程度，exp(log(p))仍旧受到浮点数的影响。因此，我们并不认为这是一个好的处理方式，甚至认为这种做法除了增加运算量外，没有其它的实际意义。有人曾引进一个标度因子，我们认为这种方法相比较理想，而且，我们也自编了相应的程序，发现它的确可以克服计算机浮点数有限的限制。为此，我们这里详细介绍这个方法。
该方法的基本原理是将每一步（或每一时刻）进行归一化计算，具体为：
每一步的前向概率进行如下转换：
(6-16a)
即有
(6-16b)
通过这样的转换，我们可以得到式（6-15b）的迭代公式为：
(6-17)
通过这样的转换，要求转换后的符合如下条件：
(6-18)
根据公式（6-18），我们可以得到式（6-17）中的归一