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

《第6章隐马尔科夫链模型2.doc》由会员分享，可在线阅读，更多相关《第6章隐马尔科夫链模型2.doc（9页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、如有侵权，请联系网站删除，仅供学习与交流第三节第四节第五节第六节第七节第八节第九节 第6章隐马尔科夫链模型2【精品文档】第 246 页第十节 后验解码问题 前面我们介绍了在给定一个符号序列时，在不同的可能状态序列中找出概率（可能性）最大的状态序列（或路径）。在实际情况中可能会碰到一个问题：不同的路径其概率相差不大，这为我们选择哪一条路径增加了一定的困难。但我们知道一条状态序列，它是由各个位置的状态组成的，比如图6-11中，从第1个位置到第45个位置均由状态F组成，紧接后面的21个是L状态，依此类推。因此克服上述困难我们可以每个位置各个不同状态的概率，这里面我们必须分清一点：某一条序列对应于某个

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

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

4、的路径（状态序列）很多，而且随着符号序列的长度的增加，姿态态序列的个数N呈指数函数的方式往上增长，以CpG岛为例：当序列长度为3时，其可能的状态序列个数为，当序列长度为4时，则有。因此如何计算（6-11）则需要一定的技巧。借鉴Viterbi算法，我们也不妨来个沿着符号序列逐步计算，这种逐步的方式有两种：一种是从符号序列的左端向右端（DNA序列的5端到3端，蛋白质序列的N端到C端），相应的算法叫前向算法（forward algorithm）；另一种是从右端向左端，相应的算法为后向算法(backward algorithm)。我们首先介绍前向算法。一、前向算法(Forward Algorithm)

5、 在介绍前向算法之前，我们先看一下式（6-7）的第一个公式即： （6-13）从中我们可发现它取的是每个位置最大值概率，因为对符号序列，其每个位置的概率等其各个状态概率之和，因此，整个符号序列的概率是每个位置概率之积，沿着符号序列的方向，我们有： (6-14)整个算法可写成如下：初始化： (6-15a)迭代过程： (6-15b)终止： (6-15c)我们可将前向算法以图示方式表示：图6-17 前向算法的基本框架图图6-18 前向算法的上体计算图。这里我们仅选状态1的计算过程，目的是使图形显得简洁明了，因为其它状的算法与此是相同的。计算例子：符号序列“3151”。我们仍以Viterbi算法中的例子

6、即符号序列“3151”来说明：计算第一个点：初始条件如Viterbi算法所示，同时设定。首先，我们计算第一个位置即，我们有：当时，我们有：当，我们有：当时，我们有：在终止阶段，我们有：与Viterbi算法相似，前向算法也存着“因计算机浮点数的限制导致后面的结果失真”，因此同样的需要用对数运算来代替相应的乘法运算，有文献曾试图通过如下的转换：克服这一点。但我们认为这样做其实没有实质性的意义，因为中间一个“exp”的运算，如果log(p)小到一定程度，exp(log(p)仍旧受到浮点数的影响。因此，我们并不认为这是一个好的处理方式，甚至认为这种做法除了增加运算量外，没有其它的实际意义。有人曾引进一

7、个标度因子，我们认为这种方法相比较理想，而且，我们也自编了相应的程序，发现它的确可以克服计算机浮点数有限的限制。为此，我们这里详细介绍这个方法。该方法的基本原理是将每一步（或每一时刻）进行归一化计算，具体为：每一步的前向概率进行如下转换： (6-16a)即有 (6-16b)通过这样的转换，我们可以得到式（6-15b）的迭代公式为： (6-17)通过这样的转换，要求转换后的符合如下条件： (6-18)根据公式（6-18），我们可以得到式（6-17）中的归一化因子的表达式为： (6-19)应用这种方法计算前向概率与后向概率显然可以在基本上克服计算机浮点数带来的不足，因为计算得到的其值肯定在状态数倒

8、数值附近浮动，与序列长度相比，值显然要小得多，如在蛋白质二级结构中有，在CpG岛中有。应用这种校正方法有一个很显著的优点就是计算序列的概率值即，它与归一化因子的关系为： (6-20a)因此很自然地就可以将它转化为相应的对数即： (6-20b)这里提前说明一下：有关文献在介绍这种归一化因子法是，将我们刚才前面介绍的前向算法的归一化因子也用在我们接下来要介绍的后向算法中。我们通过深入的研究发现：在后向算法中直接应用前向算法中的归一化因子不是很妥当，在计算后验解码问题中“计算机浮点数有限”的问题仍未得到有效的解决。为此，我们在后向算法中重新计算归一化因子，具体地将在后向算法中介绍。虽然这种改变是微小

9、的，但它在后验解码问题及后面我们要介绍的HMM模型的概率参数即的求解算法之一的Bauch-Welch算法中会发现其具有非常独特的优越性。二、后向算法（Backward Algorithm）后向算法基本上与前向算法基本类似，只不过计算的方向不同：从符号序列的右边向左边运算（对DNA序列是从3端到5端；对蛋白质序列则是从C端到N端）。下面我们直接将其算法叙述如下：初始化： 对所有状态，有 (6-20a)迭代过程： (6-20b)终止： (6-20c)计算例子：符号序列“3151”我们仍以符号序列“3151”来说明后向算法的计算过程：计算第一个点：设定，根据初始条件的设置我们有：首先，我们计算倒数第

10、二个位置即，此时，我们有：其次，计算倒数第三个位置即时，我们有：最后，我们计算倒数第四个位置也就是第一个位置即时，我们有：这样我们计算得到的有： 我们将后向算法的结果与前向算法的结果相比，它们基本一致，其差异来自于计算过程中的舍入误差。 我们在前一小节中提到“如果将前向算法的归一化因子用于后向算法，则在后验解码问题及概率参数的求解仍未摆脱计算机浮点数有限的束缚”，并说明可以重新对后向算法各迭代步骤中的概率进行归一化，具体的方法很简单，其归一化因子为： (6-21)这里的代表后向算法中的归一化因子，转化后的后向概率即与其实际的概率即的关系为： (6-22)这样我们就有： (6-23)三 后验解码

11、问题的算法本节一开头我们就提到过：计算后验解码问题需要前向算法与后向算法。现在我们已经介绍完它们，自然地，接下来就要计算后验解码问题即计算：根据条件概率公式： (6-24)我们有： (6-25)而上式右边的分子则为：于是式（6-20）可转化为： (6-26)当然可以用前向算法或后向算法得到。从式（6-25）我们可以看出的具体意义就是对某个符号序列，当第个位置的符号为状态时的概率。在CpG岛问题中，如果为“+”，则说明这个位置的碱基在CpG岛区域内；如果为“-”，则说明这个位置的碱基不在CpG岛区域内。现在我们根据式（6-26）来说明针对后向算法我们自己提出的归一化因子在后验解码中的意义。将实际

12、的前向概率与后向概率即与分别应用归一化后的概率即与来表示，即将式(6-16b)，（6-20a）及（6-23）代入式（6-26），我们有： （6-27）由于一般会在“1”附近，因此上式中它们的乘积部分显然不会太小或也不会太大，即可以将它们的乘积控制在计算机浮点数允许的范围内。自然就会摆脱了计算机浮点数有限的束缚。按老办法，我们现在仍以前面的例子即符号序列“来说明“后验解码”的计算方式。后验解码的计算例子显然我们有，根据我们前面计算得到的一系列与，我们很容易就可以计算出各个位置中各个状态的概率即。当时，我们有：当时，我们有：当时，我们有：当时，我们有：以上我们仅选了前4个点的计算，为了使读者对全部

13、数据序列有一个比较直观的认识，我们应用C语言编制了一个计算机程序，将应用前向算法和后向算法计算得到的各位置上的概率即列于图6-15，同时将也将各个位置上对应的各状态概率即总结在图6-16中。同时，应用前向算法与后向算法计算得到该序列的概率均为： 读者如果对隐马尔科夫链方法感兴趣，可以根据我们提供的例子手算几个位置，这样可以加深对隐马尔科法的印象。而且，正如我们前面所说的，一旦我们熟悉了具体的算法，则可以很“自然”将它们应用到我们熟悉的其的领域中去，显然对扩大隐马尔科夫链方法的应用以及生物信息学新方法的建立是很有帮助的。此外，根据图6-20，我们显然也可以据此判断每个位置是“F态”还是“L态”。

14、具体的方法是如果位置的某个状态概率大于另一个状态概率，则可以确定张三用的是这个状态所代表的骰子，具体为： 如果，则该位置的状态为F，反之则为L我们将上面的结果图示如下： -3.00 -5.32 -7.58 -9.76 -11.86 -12.27 -14.57 -16.81 -17.35 -19.68 -21.96 -22.55 -23.29 -25.67 -28.03 -30.35 -32.60 -34.79 -36.88 -38.90 -40.85 -42.76 -44.64 -46.50 -46.74 -49.00 -51.20 -53.30 -53.72 -56.02 -58.26 -6

15、0.41 -62.49 -64.48 -66.42 -68.32 -68.58 -70.85 -73.05 -75.16 -77.19 -79.16 -81.07 -82.95 -84.82 -85.06 -87.32 -89.52 -90.02 -92.33 -2.48 -4.27 -6.07 -7.89 -9.72 -11.55 -13.34 -15.16 -16.98 -18.75 -20.56 -22.37 -24.13 -25.76 -27.49 -29.27 -31.08 -32.90 -34.73 -36.56 -38.39 -40.23 -42.06 -43.90 -45.73

16、 -47.54 -49.36 -51.18 -53.01 -54.81 -56.62 -58.44 -60.27 -62.10 -63.93 -65.77 -67.60 -69.41 -71.23 -73.05 -74.88 -76.72 -78.55 -80.38 -82.22 -84.05 -85.86 -87.68 -89.50 -91.29-516.29 -514.19 -512.00 -509.75 -507.43 -506.68 -504.42 -502.10 -501.35 -499.08 -496.76 -496.01 -495.36 -493.52 -491.68 -489.

17、83 -487.97 -486.10 -484.20 -482.27 -480.27 -478.20 -476.04 -473.80 -473.11 -471.08 -468.97 -466.77 -466.11 -464.25 -462.36 -460.44 -458.47 -456.43 -454.31 -452.10 -451.43 -449.52 -447.55 -445.52 -443.41 -441.21 -438.94 -436.61 -434.25 -433.47 -431.11 -428.73 -427.94 -425.55-515.35 -513.52 -511.70 -5

18、09.89 -508.11 -506.38 -504.57 -502.79 -501.06 -499.25 -497.47 -495.74 -493.93 -492.10 -490.26 -488.43 -486.59 -484.76 -482.92 -481.09 -479.26 -477.43 -475.61 -473.80 -472.00 -470.17 -468.35 -466.53 -464.72 -462.89 -461.05 -459.22 -457.39 -455.56 -453.73 -451.92 -450.11 -448.28 -446.45 -444.62 -442.7

19、9 -440.97 -439.17 -437.39 -435.67 -434.08 -432.34 -430.71 -429.37 -427.78图6-19 前50个位置应用前向算法和后向算法的概率：显然，如果将换成，换成则“掷骰子问题”即刻变成蛋白质二级结构预测问题了。结合前面最大路径的求解，我们马上可以得到蛋白质二级结构预测的两个方法。足见将数学方法彻底搞懂的好处。图6-20 300个位置上各种状态（F与L）的概率四、后验解码问题的应用后验解码问题的应用之一就是我们前面刚提到的选取某个位置概率最大的状态，这些状态同样也可以组成一条路径即： （6-28）我们将式（6-18）与最大几率状态序列

20、（路径）相比： （6-28）可以发现：后验解码问题将各个位置的状态概率分开来考虑，其路径是由每个位置概率最大的状态组成，但由此路径与符号序列组成的联合概率不一定最大，而最大几率状态路径是指整条序列与相应的路径的联合概率最大，但具体到局部某个位置就不一定是最大的了。很明显，当几条路径与某个符号序列的联系概率差不多大时，选择哪条路径就会变成一件较为困难的事。这个时候，后验解码问题的优势就显示出来了。因此最好的办法是将二者联合起来考虑，这样就会形成“优势互补”的状态了。后验解码的另一个应用就是根据各个位置中各个可能状态的概率判断这个位置来源于哪个状态。其实这个问题与刚才我们讲的具有很多相似的问题。在

21、这里具体的做法是给一个阈值函数： (6-29)这里我是状态集中的某个子集，比如在蛋白质二级预测问题中，状态集，则其子集可以为：等，这样以为例，当或时，而C由于不属于这个子解，因此有。应用这个阈值函数进行如下计算： (6-30)这样我们就可以确定位置i的状态是否来自某个子状态集了。 到目前为止，我们已经详细介绍了隐马尔科夫链的两个问题，也就是说我们回答了隐马氏链的两个问题，归结起来是一个问题：如何求解其状态序列。虽然隐马氏链是一条高度抽象的链，但从我们前面所举的例子来看，它又是非常具体实际的：是许多实际对象抽象出来的一个数学模型。接下来我们要面临的是隐马氏链的第三个问题：在许多实际情况中，隐马氏模型中的状态转换概率与发射概率（统称为）往往是未知的，因此在应用它用于解决实际问题是首先得知道。这就是我们通常所说的隐马氏模型的参数估计问题，也就是如何求解隐马氏链三要素中的发射概率与状态转移概率即。