第6章隐马尔科夫链模型2

239第三节后验解码问题前面我们介绍了在给定一个符号序列时，在不同的可能状态序列中找出概率（可能性）最大的状态序列（或路径）。在实际情况中可能会碰到一个问题：不同的路径其概率相差不大，这为我们选择哪一条路径增加了一定的困难。但我们知道一条状态序列，它是由各个位置的状态组成的，比如图6-11中，从第1个位置到第45个位置均由状态F组成，紧接后面的21个是L状态，依此类推。因此克服上述困难我们可以每个位置各个不同状态的概率，这里面我们必须分清一点：某一条序列对应于某个状态序列的概率含义与状态序列中某个位置上具体某个状态的概率是不同的，我们将它们的区分总结于如图6-14：图6-14某个DNA序列对应的不同可能的CpG岛分布（状态序列）的概率,xP与各个位置上某个状态的概率之间的区别（“+”代表CpG岛区域；“-”代表非CpG岛区域）。我们在上节的“Viterbi”算法中重点介绍的是如何计算图6-14中右边部分，并将其中最大的概率对应的路径*找出来。而我们这一节要介绍的是如何计算图6-14中最下面的概率，这就是通常我们所说的“后验解码问题”：计算概率xkPi|为多少？为了计算这个概率，首先得介绍所给定序列对应的概率，由于不同的状态序列可以产生同一条序列，因此，我们有图6-15的描述：图6-15隐马尔克夫链中符号序列的概率与路径的关系（这里假设符号序列集只有两条序列1x与2x，对应的状态序列（路径）中也只有两个路径1与2。240图6-15中的1xP与2xP就是要求的某个符号序列的概率，由于此时路径1与2已知，因此它与符号序列的联合概率表达式可写为：LiiLiiLiiLiiiiiiiiiiiiiiaxeaxPaxeaxPaxeaxPaxeaxP12022120121102111011212222111111212221111111,,,,（6-11）公式（6-11）罗列了一大堆符号，这对学生物学的人来说要比较快地接受它们可能会要多花点时间，为使我们有一个清晰的直观印象，我们这里以两条DNA序列及相应的CpG岛分布图来说明。图6-16DNA序列中CpG岛分布的隐马尔科夫链模型中在公式（6-6）中的示意图241如果我们以x代表Mxxx,,,21中的某个序列，相应的状态序列有N,,,21，则有：NiixPxP1,(6-12)在实际情况中，可能的路径（状态序列）很多，而且随着符号序列的长度的增加，姿态态序列的个数N呈指数函数的方式往上增长，以CpG岛为例：当序列长度为3时，其可能的状态序列个数为823N，当序列长度为4时，则有1624N。因此如何计算（6-11）则需要一定的技巧。借鉴Viterbi算法，我们也不妨来个沿着符号序列逐步计算，这种逐步的方式有两种：一种是从符号序列的左端向右端（DNA序列的5’端到3’端，蛋白质序列的N端到C端），相应的算法叫前向算法（forwardalgorithm）；另一种是从右端向左端，相应的算法为后向算法(backwardalgorithm)。我们首先介绍前向算法。一、前向算法(ForwardAlgorithm)在介绍前向算法之前，我们先看一下式（6-7）的第一个公式即：klkkillaivxeiv1max（6-13）从中我们可发现它取的是每个位置最大值概率，因为对符号序列x，其每个位置的概率等其各个状态概率之和，因此，整个符号序列x的概率xP是每个位置概率之积，沿着符号序列的方向，我们有：kklkillaifxeif11(6-14)整个算法可写成如下：初始化0i：000,100kffk(6-15a)迭代过程Li,,2,1：klkkillaifxeif1(6-15b)终止：kkkaLfxP0(6-15c)我们可将前向算法以图示方式表示：图6-17前向算法的基本框架图242图6-18前向算法的上体计算图。这里我们仅选状态1的计算过程，目的是使图形显得简洁明了，因为其它状的算法与此是相同的。计算例子：符号序列“3151”。我们仍以Viterbi算法中的例子即符号序列“3151”来说明：计算第一个点：初始条件如Viterbi算法所示，同时设定5.0,5.0,,FbeginLbeginaa。首先，我们计算第一个位置即31x，我们有：08333.05.016103105.05.011.0031,1,1FbeginbeginFFLbeginbeginLLafxefafxef当2i时，12x，我们有：2232104027.195.008333.01.005.0611112109167.405.008333.09.005.01.01112FFFLFLFFFLFLLLLLafafxefafafxef当3i，53x，我们有：24332334233103028.295.0104027.11.0109167.4612253101264.505.0104027.19.0109167.41.02253FFFLFLFFFLFLLLLLafafxefafafxef当4i时，14x，我们有：43445344107315.395.0103028.21.0101264.5613314107652.505.0103028.29.0101264.51013314FFFLFLFFFLFLLLLLafafxefafafxef在终止阶段，我们有：445,,101540.25.0107315.35.0107652.544endFFendLLafafxP与Viterbi算法相似，前向算法也存着“因计算机浮点数的限制导致后面的结果失真”，因此同样的需要用对数运算来代替相应的乘法运算，有文献曾试图通过如下的转换：qpqplogexplogexploglog克服这一点。但我们认为这样做其实没有实质性的意义，因为中间一个“exp”的运算，如果log(p)小到一定程度，exp(log(p))仍旧受到浮点数的影响。因此，我们并不认为这是一个好的处理方式，甚至认为这种做法除了增加运算量外，没有其它的实际意义。有人曾引进一个标度因子，我们认为这种方法相比较理想，而且，我们也自编了相应的程序，发现它的确可以克服计算机浮点数有限的限制。为此，我们这里详细介绍这个方法。该方法的基本原理是将每一步（或每一时刻）进行归一化计算，具体为：每一步的前向概率进行如下转换：ijjllsifif1~(6-16a)即有ijjllsifif1~(6-16b)通过这样的转换，我们可以得到式（6-15b）的迭代公式为：kklkililaifxesif~11~11(6-17)通过这样的转换，要求转换后的ifl~符合如下条件：2441~llif(6-18)根据公式（6-18），我们可以得到式（6-17）中的归一化因子is的表达式为：lkklkiliaifxes~11(6-19)应用这种方法计算前向概率与后向概率显然可以在基本上克服计算机浮点数带来的不足，因为计算得到的ifl~其值肯定在状态数q倒数值q1附近浮动，与序列长度相比，q值显然要小得多，如在蛋白质二级结构中有3q，在CpG岛中有2q。应用这种校正方法有一个很显著的优点就是计算序列x的概率值即xP，它与归一化因子的关系为：LjjsxP1(6-20a)因此很自然地就可以将它转化为相应的对数即：jjsxPlog(6-20b)这里提前说明一下：有关文献在介绍这种归一化因子法是，将我们刚才前面介绍的前向算法的归一化因子也用在我们接下来要介绍的后向算法中。我们通过深入的研究发现：在后向算法中直接应用前向算法中的归一化因子is不是很妥当，在计算后验解码问题中“计算机浮点数有限”的问题仍未得到有效的解决。为此，我们在后向算法中重新计算归一化因子，具体地将在后向算法中介绍。虽然这种改变是微小的，但它在后验解码问题及后面我们要介绍的HMM模型的概率参数即的求解算法之一的Bauch-Welch算法中会发现其具有非常独特的优越性。二、后向算法（BackwardAlgorithm）后向算法基本上与前向算法基本类似，只不过计算的方向不同：从符号序列的右边向左边运算（对DNA序列是从3’端到5’端；对蛋白质序列则是从C端到N端）。下面我们直接将其算法叙述如下：初始化Li：对所有状态k，有0kkaLb(6-20a)迭代过程1,,1Li：llilklkibxeaib11(6-20b)终止：llllbxeaiP110(6-20c)计算例子：符号序列“3151”我们仍以符号序列“3151”来说明后向算法的计算过程：计算第一个点：设定5.000FLaa，根据初始条件的设置我们有：2455.044FLbb首先，我们计算倒数第二个位置即53x，此时，141xxi我们有：244244101667.8615.095.01015.005.014143103333.5615.01.01015.09.014143xebaxebabxebaxebabFFFFLLFLFFFLFLLLLL其次，计算倒数第三个位置即2i时，12x，531xxi我们有：2223332233103197.161101667.895.0101103333.505.053532101611.661101667.81.0101103333.59.053532xebaxebabxebaxebabFFFFLLFLFFFLFLLLLL最后，我们计算倒数第四个位置也就是第一个位置即1i时，1,321xx，我们有：3232242322101203.261103197.195.0101101611.605.012121107445.761103197.11.0101101611.69.012121xebaxebabxebaxebabFFFFLLFLFFFLFLLLLL这样我们计算得到的XP有：4341,1,101541.2101203.2615.0107445.71.05.011FFFbeginLLLbeginbxeabxeaxP我们将后向算法的结果与前向算法的结果相比，它们基本一致，其差异来自于计算过程中的舍入误差。我们在前一小节中提到“如果将前向算法的归一化因子用于后向算法，则在后验解码问题及概率参数的求解仍未摆脱计算机浮点数有限的束缚”，并说明可以重新对后向算法各迭代步骤中的概率ibk进行归一化，具体的方法很简单，其归一化因子bis为：kllilklbiibxeas1~1(6-21)246这里的bis代表后向算法中的归一化因子，转化后