美文网首页
序列比对（十七）——第二部分的小结

序列比对（十七）——第二部分的小结

作者: 生信了 | 来源:发表于2019-11-14 11:32 被阅读0次

原创：hxj7
序列比对的系列文章第二部分主要介绍了HMM（隐马尔科夫模型），包含了八篇文章：
《序列比对（九）从掷骰子说起HMM》
《序列比对（十）viterbi算法求解最可能路径》
《序列比对（11）计算符号序列的全概率》
《序列比对（12）：计算后验概率》
《序列比对（13）后验解码》
《序列比对（14）viterbi算法和后验解码的比较》
《序列比对（15）EM算法以及Baum-Welch算法的推导》
《序列比对（16）Baum-Welch算法估算HMM参数》

HMM模型最关键的一点就是在一个状态序列中，某一步状态的概率只与上一步的状态有关。也正是因为与前面一步状态有关，所以HMM模型天然地适用动态规划算法。

这几篇文章都是以掷骰子为例。

一、如何随机生成状态序列和符号序列：

假定要随机生成一个状态序列\pi和一个符号序列x，那么：
P(\pi_i|\pi_1,\pi_2,...,\pi_{i-1}) = P(\pi_i|\pi_{i-1}) \tag{1.1}
P(x_i|x_1,x_2,...,x_{i-1},\pi_1,\pi_2,...,\pi_i) = P(x_i|\pi_{i}) \tag{1.2}

给出了转移概率和发射概率的符号定义：
a_{kl} = P(\pi_i=l|\pi_{i-1}=k) \tag{1.3}
e_k(b) = P(x_i=b|\pi_i=k) \tag{1.4}

若是给初始状态0和初始符号x_0建模，则有：
e_0(b) = \begin{cases} 0, & \text{if $b \neq x_0$} \\ 1, & \text{if $b = x_0$} \end{cases} \tag{1.5}

若是给结束状态0建模，我们假定：
a_{k0} = 1, \ \ \ \text{for each possible $k \neq 0$} \tag{1.6}

二、Viterbi算法求解最可能路径

一条路径就是一个状态序列。在符号序列已知的情况下，如何猜出来最有可能的路径\pi^*。即

\pi^* = \underset{\pi}{\mathrm{argmax}} \, P(x, \pi) \tag{2.1}

这个可以用Viterbi算法求解。Viterbi算法是一种动态规划算法。假设v_k(i)是以状态k以及符号x_i结束的所有序列中概率的最大值。那么
Viterbi算法的初始条件是：
v_k(0) = \begin{cases} 0, & \text{if $k \neq 0$} \\ 1, & \text{if $k = 0$} \end{cases} \tag{2.2}

迭代公式是：
v_l(i+1) = e_l(x_{i+1}) \underset{k}{\max} \big[ v_k(i)a_{kl} \big] \\ \mathrm{ptr}_{i+1}(l) = \underset{k}{\mathrm{argmax}} \big[ v_k(i)a_{kl} \big] \tag{2.3}

终止条件是：
P(x,\pi^*) = \underset{k}{\max} \, v_k(L) \\ \pi_L^* = \underset{k}{\mathrm{argmax}} \, v_k(L) \tag{2.4}

回溯公式是：
\pi_i^* = \mathrm{ptr}_{i+1}(\pi_{i+1}^*) \tag{2.5}

三、前向算法和后向算法计算符号序列的全概率

假设符号序列x已知，而路径\pi未知，那么符号序列的全概率P(x)如何计算？就是
P(x) = \sum_{\pi} P(x, \pi) \tag{3.1}

具体计算仍然使用动态规划迭代计算，按照序列迭代的方向可分为前向算法和后向算法。

所谓前向算法，就是从序列头部开始迭代，直至序列尾部。其初始条件是：
f_k(i) = P(x_1,x_2,...,x_i,\pi_i=k) \\ f_k(0) = \begin{cases} 0, & \text{if $k \neq 0$} \\ 1, & \text{if $k = 0$} \end{cases} \tag{3.2}

迭代公式是：
f_l(i+1) = e_l(x_{i+1}) \sum_k f_k(i)a_{kl} \tag{3.3}

终止条件是：
P(x) = \sum_k f_k(L) \tag{3.4}

所谓后向算法，就是从序列尾部开始迭代，直至序列头部。其初始条件是：
b_k(i) = P(x_{i+1},x_{i+2},...,x_L|\pi_i=k) \\ b_k(L) = 1, \ \ \ \ \ \text{for each possible $k > 0$}\tag{3.5}

迭代公式是：
b_k(i) = \sum_l a_{kl}e_l(x_{i+1})b_l(i+1) \tag{3.6}

终止条件是：
P(x) = b_0(0) \tag{3.7}

在实际的计算中，为了解决多个概率相乘可能导致的下溢问题，要对概率进行适当的缩放。我们选择了缩放因子这一方法。我们首先定义缩放后的概率：
\tilde{f}_l(i) = \frac{f_l(i)}{\displaystyle \prod_{j=1}^i s_j} \\ \tilde{b}_l(i) = \frac{b_l(i)}{\displaystyle \prod_{j=i}^L s_j} \tag{3.8}

缩放因子s_i的计算方法是：
s_{i+1} = \sum_l e_l(x_{i+1}) \sum_k \tilde{f}_k(i)a_{kl}, \ \ \ \text{when $\sum_l \tilde{f}_l(i) = 1$} \tag{3.9}

并且有：
P(x) = \prod_{j=1}^L s_j \tag{3.10}

用缩放后的概率表示前向算法：
\tilde{f}_k(0) = \begin{cases} 0, & \text{if $k \neq 0$} \\ 1, & \text{if $k = 0$} \end{cases} \\ s_{i+1} = \sum_l e_l(x_{i+1}) \sum_k \tilde{f}_k(i)a_{kl} \\ \tilde{f}_l(i+1) = \frac{1}{s_{i+1}} e_l(x_{i+1}) \sum_k \tilde{f}_l(i) a_{kl} \tag{3.11}

用缩放后的概率表示后向算法：
\tilde{b}_k(L) = \frac{1}{s_L}\\ \tilde{b}_k(i) = \frac{1}{s_i} \sum_l a_{kl} \tilde{b}_l(i+1)e_l(x_{i+1}) \tag{3.12}

四、后验概率的计算与后验解码

当符号序列已知而路径未知时，我们除了关心最可能路径外，也对某一位置状态的概率P(\pi_i=k|x)感兴趣。很明显，这是一个后验概率，通过推导，可以通过下面的公式计算：
P(\pi_i=k|x) = \frac{f_k(i)b_k(i)}{P(x)} \tag{4.1}

如果利用缩放后的概率计算，上述公式变为：
P(\pi_i=k|x) = \tilde{f}_k(i) \tilde{b}_k(i) s_i \tag{4.2}

所谓后验解码，就是对每一个位置，计算得到该位置最有可能的状态。即
\begin{align} \hat{\pi}_i & = \underset{k}{\mathrm{argmax}} \, P(\pi_i=k|x) \\ & = \underset{k}{\mathrm{argmax}} \, \tilde{f}_k(i) \tilde{b}_k(i) \tag{4.3} \end{align}

五、Baum-Welch算法估算HMM模型参数

EM算法是一种在有缺失数据的情况下利用最大似然法或最大对数似然法估算概率模型参数的迭代方法。
\hat{\theta} = \underset{\theta}{\mathrm{argmax}} \, P(x|\theta) \\ \hat{\theta} = \underset{\theta}{\mathrm{argmax}} \, \log P(x|\theta) \tag{5.1}

假设缺失数据是y，那么EM算法的步骤就是：
E步骤（Expectation）要计算Q函数：
Q(\theta|\theta^t) = \sum_y P(y|x,\theta^t) \log P(x,y|\theta) \tag{5.2}

M步骤（Maximization）根据步骤t的结果计算得到步骤t+1的参数：
\theta^{t+1} = \underset{\theta}{\mathrm{argmax}} \, Q(\theta|\theta^t) \tag{5.3}

终止条件就是当（对数）似然值变化小于一个阈值T或者迭代次数超过最大迭代次数MaxIter，即：
\begin{aligned} & P(x|\theta^{t+1}) - P(x|\theta^t) < T \\ & \text{\textbf{or}} \ \ \ \log P(x|\theta^{t+1}) - \log P(x|\theta^t) < T \\ & \text{\textbf{or}} \ \ \ t > MaxIter \tag{5.4} \end{aligned}

Baum-Welch算法是EM算法的一个特例。当符号序列已知而路径未知时，可以用Baum-Welch算法来估算HMM模型的参数。这时，缺失数据就是路径\pi。那么，Baum-Welch算法的Q函数变为：
Q(\theta|\theta^t) = \sum_{\pi} P(\pi|x,\theta^t) \log P(x,\pi|\theta) \tag{5.5}

HMM模型的参数主要是转移概率a_{kl}以及发射概率e_k(b)，可以证明，当Baum-Welch算法采用如下迭代方式计算参数时，可以让对数似然取到局部最大值：
a_{kl}^{t+1} = \frac{A_{kl}^t + r_{kl}}{\displaystyle \sum_{l'} \Big[A_{kl'}^t + r_{kl'} \Big]} \\ e^{t+1}_k(b) = \frac{E^{t+1}_k(b) + r_k(b)}{\displaystyle \sum_{b'} \Big[E^{t+1}_k(b') + r_k(b') \Big]} \tag{5.6}

其中，t表示迭代的次数，A_{kl}是a_{kl}在所有训练序列中发生的期望次数，E_k(b)是e_k(b)在所有训练序列中发生的期望次数。即
\begin{align} \displaystyle A_{kl} & = \sum_{j} \sum_{\pi} P(\pi^j|x^j,\theta) A_{kl}(\pi^j) \\ & = \sum_{j} \sum_i P(\pi_i^j=k, \pi_{i+1}^j=l|x^j,\theta) \tag{5.7} \end{align}

\begin{align} \displaystyle E_k(b) & = \sum_j \sum_\pi P(\pi^j|x^j, \theta) E_k(b, \pi^j) \\ & = \sum_{j} \sum_i P(\pi_i^j=k, x_i^j=b|x^j,\theta) \\ & = \sum_{j} \sum_{\{i|x_i^j=b\}} P(\pi_i^j=k|x^j,\theta) \tag{5.8} \end{align}

如果用前向算法和后向算法中的符号表示的话，我们可以推导得到：
\begin{align} \displaystyle A_{kl} & = \sum_j \frac {1} {P(x^j|\theta)} \sum_i f^j_k(i) b^j_l(i+1) a_{kl} e_l(x^j_{i+1}) \\ & = \sum_j \sum_i \tilde{f}^j_k(i) a_{kl} e_l(x^j_{i+1}) \tilde{b}^j_l(i+1) \tag{5.9} \end{align}

\begin{align} \displaystyle E_k(b) & = \sum_{j} \frac {1} {P(x^j|\theta)} \sum_{\{i|x_i^j=b\}} f^j_k(i) b_k^j(i)\\ & = \sum_{j} \sum_{\{i|x_i^j=b\}} \tilde{f}^j_k(i) \tilde{b}^j_k(i) s^j_i \tag{5.10} \end{align}

伪计数之和反映了我们?认为先验知识的可靠程度。如果伪计数之和较大，则我们认为先验概率比较可靠，需要用更多的数据?取改变它。反之不太可靠。最开始的概率参数可以用随机数?。?只要保证相应的概率和为1即可。

至此，序列比对系列文章已经介绍了“第一部分：二序列比对的常见问题和相应算法”以及“第二部分：HMM模型及相关算法简介”。在学习“将HMM模型应用到序列比对（或序列发现）”之前，我们接下来准备先学习一下多序列比对的常见问题和相关算法。

（公众号：生信了）

相关文章

网友评论

      本文标题：序列比对（十七）——第二部分的小结

      本文链接：https://www.haomeiwen.com/subject/elxrictx.html

      延伸阅读

      深度阅读

      • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

      栏目导航

      热点阅读

      关于我们|服务条款|联系我们|序列比对（十七）——第二部分的小结|投稿指南|网站地图|RSS订阅|排版工具|手机版

      提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

      Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

      备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

      本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

      所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！