序列比对（17）第二部分的小结

最新推荐文章于 2021-01-14 15:28:14 发布

生信了（公众号同名）

最新推荐文章于 2021-01-14 15:28:14 发布

阅读量330

点赞数

分类专栏： # 序列算法文章标签：序列比对 HMM

本文链接：https://blog.csdn.net/biocity/article/details/99597791

版权

序列算法专栏收录该内容

32 篇文章 16 订阅

订阅专栏

序列比对的系列文章第二部分主要介绍了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)$ 如何计算？就是
$\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}$

终止条件是：
$\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}$

终止条件是：
$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}$

并且有：
$\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{aligned} \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{aligned}$

五、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$ 或者迭代次数超过最大迭代次数 $M a x I t e r$ ，即：
$\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{aligned} \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{aligned}$

$\begin{aligned} \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{aligned}$

如果用前向算法和后向算法中的符号表示的话，我们可以推导得到：
$\begin{aligned} \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{aligned}$

$\begin{aligned} \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{aligned}$