序列比对(十五)——EM算法以及Baum-Welch算法的推导

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

序列比对(十五)——EM算法以及Baum-Welch算法的推导
EM算法推导
序列比对(十六)——Baum-Welch算法估算HMM参数
3️⃣ 多序列比对(2):方法
EM算法
EM 算法
全局比对与局部比对
序列比对(三)——局部联配Smith-Waterman算法
算法文献阅读14：Waterman's new algorith
生物信息学

原创：hxj7

本文介绍了 EM算法和 Baum-Welch算法的推导过程。

一般地，估算概率模型参数可以用最大似然法，即找到
$\hat{\theta} = \mathrm{argmax}_\theta \, P(x|\theta)$

有时为了简便运算，也可以用最大对数似然代替最大似然，即以下公式
$\hat{\theta} = \mathrm{argmax}_\theta \, \mathrm{log} \, P(x|\theta)$

很多时候上述公式没有解析解或者解析解运算过于复杂，这个时候可以用迭代的方法求解。预先设置终止条件，满足该条件时即停止迭代。

EM算法

EM算法就是这样一种迭代算法，该算法在有缺失值存在的情况下估算概率模型的参数（或参数集，下文统称参数）。其大致过程如下：
E步骤：计算Q函数。
M步骤：相较于 $\theta$ ，最大化 $Q(\theta|\theta^t)$ 。

其具体推导过程如下：
最终目的是找到最大对数似然对应的参数。
$\hat{\theta} = \mathrm{argmax}_\theta \, \mathrm{log} \, P(x|\theta) \tag{1}$

假设 $y$ 为缺失数据，我们首先可以得到：
$\mathrm{log} \, P(x|\theta) = \mathrm{log} \, P(x,y|\theta) - \mathrm{log} \, P(y|x,\theta) \tag{2}$

公式(2)很容易推导，因为：
$P(x,y|\theta) = P(x|\theta)P(y|x,\theta) \tag{2.1}$
所以
$P(x|\theta) = \frac{P(x,y|\theta)}{P(y|x,\theta)} \tag{2.2}$
等式两边取 $\mathrm{log} \,$ 可以得到公式(2)。

在求取最大对数似然的迭代过程中，假设步骤 $t$ 中得到了参数 $\theta^t$ ，对应的对数似然是 $\mathrm{log}\, P(x|\theta^t)$ 。那么步骤 $t+1$ 中的对数似然应该不小于 $\mathrm{log}\, P(x|\theta^t)$ 。即
$\mathrm{log}\, P(x|\theta^{t+1}) - \mathrm{log}\, P(x|\theta^t) \geq 0 \tag{3}$

为了得到 $\theta^{t+1}$ ，我们首先变换得到：
$\begin{align} \displaystyle \mathrm{log}\, P(x|\theta)= &\sum_y P(y|x,\theta^t)\mathrm{log}\, P(x,y|\theta)\\ &-\sum_y P(y|x,\theta^t)\mathrm{log}\, P(y|x,\theta) \tag{4} \end{align}$

公式（4）可以这样推导得到：
将公式（2）等式两边乘上 $P(y|x,\theta^t)$ ，再对所有的 $y$ 求和，得到：
$\begin{align} \displaystyle \sum_y P(y|x,\theta^t)\mathrm{log} \, P(x|\theta) = & \sum_y P(y|x,\theta^t)\mathrm{log} \, P(x,y|\theta)\\ & -\sum_y P(y|x,\theta^t)\mathrm{log} \, P(y|x,\theta) \tag{4.1} \end{align}$
很容易看出等式左边：
$\begin{align} \displaystyle \sum_y P(y|x,\theta^t)\mathrm{log} \, P(x|\theta) & = \mathrm{log} \, P(x|\theta)\sum_y P(y|x,\theta^t) \\ & = \mathrm{log} \, P(x|\theta) \tag{4.2} \end{align}$
由公式（4.1）和（4.2）可以得到公式（4）。

如果令
$Q(\theta|\theta^t) = \sum_y P(y|x,\theta^t)\mathrm{log} \, P(x,y|\theta) \tag{5}$

那么我们可以得到：
$\begin{align} \displaystyle & \mathrm{log} \, P(x|\theta) - \mathrm{log} \, P(x|\theta^t)\\ & = Q(\theta|\theta^t) - Q(\theta^t|\theta^t) + \sum_y P(y|x,\theta^t)\mathrm{log} \, \frac{P(y|x,\theta^t)}{P(y|x,\theta)} \tag{6} \end{align}$

公式（6）的推导也很容易：
$\begin{align} &\mathrm{log} \, P(x|\theta) - \mathrm{log} \, P(x|\theta^t)\\ &=\bigg[\sum_y P(y|x,\theta^t)\mathrm{log} \, P(x,y|\theta) - \sum_y P(y|x,\theta^t)\mathrm{log} \, P(y|x,\theta) \bigg]\\ &\ \ \ \ \ - \bigg[\sum_y P(y|x,\theta^t)\mathrm{log} \, P(x,y|\theta^t) - \sum_y P(y|x,\theta^t)\mathrm{log} \, P(y|x,\theta^t) \bigg]\\ &=\bigg[Q(\theta|\theta^t) - \sum_y P(y|x,\theta^t)\mathrm{log} \, P(y|x,\theta)\bigg] \\ &\ \ \ \ \ - \bigg[Q(\theta^t|\theta^t) - \sum_y P(y|x,\theta^t)\mathrm{log} \, P(y|x,\theta^t)\bigg] \\ &=\bigg[Q(\theta|\theta^t) - Q(\theta^t|\theta^t)\bigg]\\ &\ \ \ \ \ +\bigg[\sum_y P(y|x,\theta^t)\mathrm{log} \, P(y|x,\theta^t) - \sum_y P(y|x,\theta^t)\mathrm{log} \, P(y|x,\theta)\bigg]\\ &=Q(\theta|\theta^t) - Q(\theta^t|\theta^t) + \sum_y P(y|x,\theta^t)\mathrm{log} \, \frac{P(y|x,\theta^t)}{P(y|x,\theta)} \tag{6.1} \end{align}$

其中 $\sum_yP(y|x,\theta^t)\mathrm{log} \, \frac{P(y|x,\theta^t)}{P(y|x,\theta)}$ 一项是 $P(y|x,\theta^t)$ 对 $P(y|x,\theta)$ 的相对熵，所以不小于0。只要 $Q(\theta|\theta^t)$ 不小于 $Q(\theta^t|\theta^t)$ ，那么就可以保证公式（3）成立。那么可以取
$\theta^{t+1}=\mathrm{argmax}_\theta \, Q(\theta|\theta^t) \tag{7}$

我们定义概率分布 $P(x)$ 相对于概率分布 $Q(x)$ 的相对熵为：
$\displaystyle H(P||Q) = \sum_i\ P(x_i) \mathrm{log} \, \frac{P(x_i)}{Q(x_i)} \tag{7.1}$
那么 $H(P||Q) \geq 0$ ，证明如下：
我们知道：
$\mathrm{log} \, x \leq x - 1, \qquad when \quad x > 0 \tag{7.2}$
那么：
$-\mathrm{log} \, x \geq 1 - x, \qquad when \quad x > 0 \tag{7.3}$
所以：
$\begin{align} \mathrm{log} \, \frac{P(x_i)}{Q(x_i)} &= -\mathrm{log} \, \frac{Q(x_i)}{P(x_i)} \\ &\geq 1 - \frac{Q(x_i)}{P(x_i)} \tag{7.4} \end{align}$
因而：
$\begin{align} \displaystyle H(P||Q) &= \sum_i\ P(x_i) \mathrm{log} \, \frac{P(x_i)}{Q(x_i)}\\ &\geq \sum_i\ P(x_i) \bigg(1 - \frac{Q(x_i)}{P(x_i)}\bigg)\\ &=\sum_i\ P(x_i)-\sum_i\ Q(x_i)\\ &=1-1\\ &=0 \tag{7.5} \end{align}$
等号成立当且仅当对所有的 $i$ ， $P(x_i)=Q(x_i)$ 都成立。

Baum-Welch算法

Baum-Welch算法是EM算法的一个特例，可以用于估算HMM模型中的概率参数。其缺失的数据是路径 $\pi$ 。所以，Baum-Welch算法求解的是：
$Q(\theta|\theta^t) = \sum_\pi P(\pi|x,\theta^t)\mathrm{log} \, P(x,\pi|\theta) \tag{8}$

我们知道，HMM模型中，对于一条给定的路径 $\pi$ ，模型中的参数（如转移概率 $a_{kl}$ 、发射概率 $e_k(b)$ 等）都会在计算 $\mathrm{log}P(x,\pi|\theta)$ 的式子中出现多次。假设 $a_{kl}$ 出现的次数为 $A_{kl}$ ，而 $e_k(b)$ 出现的次数为 $E_k(b)$ ，那么对于路径 $\pi$ ：
$\displaystyle P(x,\pi|\theta) = \prod_{k=0}^M \prod_{l=1}^M a_{kl}^{A_{kl}(\pi)} \prod_{k=1}^M \prod_{b}[e_k(b)]^{E_k(b, \pi)} \tag{9}$

所以：
$\displaystyle \mathrm{log} \, P(x,\pi|\theta) = \sum_{k=0}^M \sum_{l=1}^M A_{kl}(\pi) \mathrm{log} \, a_{kl} + \sum_{k=1}^M \sum_{b} E_k(b, \pi) \mathrm{log} \, e_k(b) \tag{10}$

由此，公式（8）变成：
$\begin{align} Q(\theta|\theta^t) & = \sum_\pi P(\pi|x,\theta^t) \times \\ & \bigg[\sum_{k=0}^M \sum_{l=1}^M A_{kl}(\pi) \mathrm{log} \, a_{kl} + \sum_{k=1}^M \sum_{b} E_k(b, \pi) \mathrm{log} \, e_k(b)\bigg] \tag{11} \end{align}$

改变公式（11）的求和顺序，我们可以得到：
$Q(\theta|\theta^t) = \sum_{k=0}^M \sum_{l=1}^M A_{kl} \mathrm{log} \, a_{kl} + \sum_{k=1}^M \sum_{b} E_k(b) \mathrm{log} \, e_k(b) \tag{12}$

其中：
$A_{kl} = \sum_\pi P(\pi|x,\theta^t) A_{kl}(\pi) \\ E_k(b) = \sum_\pi P(\pi|x,\theta^t) E_k(b, \pi) \tag{12.1}$

公式（12）可以这样推导得到：
公式（11）先对 $\pi$ 求和，可以得到：
$\begin{align} Q(\theta|\theta^t) & = \sum_{k=0}^M \sum_{l=1}^M \sum_\pi P(\pi|x,\theta^t) A_{kl}(\pi) \mathrm{log} \, a_{kl} \\ & \ \ \ \ \ + \sum_{k=1}^M \sum_{b} \sum_\pi P(\pi|x,\theta^t) E_k(b, \pi) \mathrm{log} \, e_k(b)\\ & = \sum_{k=0}^M \sum_{l=1}^M \mathrm{log} \, a_{kl} \sum_\pi P(\pi|x,\theta^t) A_{kl}(\pi) \\ & \ \ \ \ \ + \sum_{k=1}^M \sum_{b} \mathrm{log} \, e_k(b) \sum_\pi P(\pi|x,\theta^t) E_k(b, \pi) \\ & = \sum_{k=0}^M \sum_{l=1}^M \mathrm{log} \, a_{kl} A_{kl} + \sum_{k=1}^M \sum_{b} \mathrm{log} \, e_k(b) E_k(b) \tag{12.2} \end{align}$

现在关键就是怎么样取 $a_{kl}$ 以及 $e_k(b)$ 的值，可以最大化 $Q(\theta|\theta^t)$ 。我们令
$a_{kl}^0 = \frac{A_{kl}}{\sum_{l'} A_{kl'}} \qquad e_k^0(b) = \frac{E_k(b)}{\sum_{b'} E_{k}(b')} \tag{13}$

那么此时对应的 $Q(\theta|\theta^t)$ 为其最大值。

公式（13）的结论证明如下：
假设 $a_{kl}^0$ 对应的 $Q(\theta|\theta^t)$ 为 $Q^0(\theta|\theta^t)$ ，其它任意 $a_{kl}$ 对应的 $Q(\theta|\theta^t)$ 就记为 $Q(\theta|\theta^t)$ 。那么：
$\begin{align} & Q^0(\theta|\theta^t) - Q(\theta|\theta^t) \\ & = \bigg[\sum_{k=0}^M \sum_{l=1}^M A_{kl} \mathrm{log} \, a_{kl}^0 + \sum_{k=1}^M \sum_{b} E_k(b) \mathrm{log} \, e_k^0(b)\bigg] \\ & \ \ \ \ \ - \bigg[\sum_{k=0}^M \sum_{l=1}^M A_{kl} \mathrm{log} \, a_{kl} + \sum_{k=1}^M \sum_{b} E_k(b) \mathrm{log} \, e_k(b)\bigg] \\ & = \bigg[\sum_{k=0}^M \sum_{l=1}^M A_{kl} \mathrm{log} \, a_{kl}^0 - \sum_{k=0}^M \sum_{l=1}^M A_{kl} \mathrm{log} \, a_{kl}\bigg]\\ & \ \ \ \ \ + \bigg[\sum_{k=1}^M \sum_{b} E_k(b) \mathrm{log} \, e_k^0(b) - \sum_{k=1}^M \sum_{b} E_k(b) \mathrm{log} \, e_k(b)\bigg]\\ & = \sum_{k=0}^M \sum_{l=1}^M A_{kl} \mathrm{log} \, \frac{a_{kl}^0}{a_{kl}} + \sum_{k=1}^M \sum_{b} E_k(b) \mathrm{log} \, \frac{e_k^0(b)}{e_k(b)}\\ & = \sum_{k=0}^M \bigg[\sum_{l'}A_{kl'}\bigg] \sum_{l=1}^M \frac{A_{kl}}{\sum_{l'}A_{kl'}} \mathrm{log} \, \frac{a_{kl}^0}{a_{kl}} \\ & \ \ \ \ \ + \sum_{k=1}^M \bigg[\sum_{b'}E_{k}(b')\bigg] \sum_{b} \frac{E_k(b)}{\sum_{b'}E_{k}(b')} \mathrm{log} \, \frac{e_k^0(b)}{e_k(b)}\\ & = \sum_{k=0}^M \bigg[\sum_{l'}A_{kl'}\bigg] \sum_{l=1}^M a_{kl}^0 \mathrm{log} \, \frac{a_{kl}^0}{a_{kl}} + \sum_{k=1}^M \bigg[\sum_{b'}E_{k}(b')\bigg] \sum_{b} e_k^0(b) \mathrm{log} \, \frac{e_k^0(b)}{e_k(b)}\\ & \geq 0 \tag{13.1} \end{align}$

最后一个不等式中 $\sum_{l=1}^M a_{kl}^0 \mathrm{log} \, \frac{a_{kl}^0}{a_{kl}}$ 可以看作 $a_{kl}^0$ 对于 $a_{kl}$ 的相对熵； $\sum_{b} e_k^0(b) \mathrm{log} \, \frac{e_k^0(b)}{e_k(b)}$ 可以看作 $e_k^0(b)$ 相对于 $e_k(b)$ 的相对熵。所以二者都不小于0。

小结

本文的推导过程基于《生物序列分析》第11章中的内容，笔者所做的工作就是将书中简要的推导过程补充详细。

（公众号：生信了）

序列比对(十五)——EM算法以及Baum-Welch算法的推导
原创：hxj7 本文介绍了 EM算法和 Baum-Welch算法的推导过程。一般地，估算概率模型参数可以用最大...
EM算法推导
推导EM算法，并证明收敛性。 Jensen’s inequality 定理：若是凸函数，是随机变量，我们有：若是...
序列比对(十六)——Baum-Welch算法估算HMM参数
原创：hxj7 本文介绍了如何用Baum-Welch算法来估算HMM模型中的概率参数。 Baum-Welch算法应...
3️⃣ 多序列比对(2):方法
序列比对和序列特征分析总目录多序列比对的算法很多，有 1 动态规划求解 2 渐进策略算法 3 基于一致性的方法 ...
EM算法
1. 前提推导 2. EM算法 EM算法就是含有隐变量的概率模型参数的极大似然估计法。延森不等式（Jensen'...
EM 算法
参考: 从最大似然到EM算法浅解（EM算法）The EM Algorithm EM算法及其推广学习笔记 EM算法...
全局比对与局部比对
Needleman-Wunsch算法和Smith-Waterman算法分别对应全局比对和局部比对，这两个生物序列比...
序列比对(三)——局部联配Smith-Waterman算法
原创：hxj7 关于全局联配的介绍可参见前文：序列比对（一）全局比对Needleman-Wunsch算法序列比对（...
算法文献阅读14：Waterman's new algorith
Efficient Sequence Alignment Algorithms 高效的序列比对算法 Journal...
生物信息学
一、生物信息学研究方向： 1.序列比对： BLAST算法、FASTA算法。 2.蛋白质比对。 3.基因...

序列比对(十五)——EM算法以及Baum-Welch算法的推导

EM算法

Baum-Welch算法

小结

相关文章

序列比对(十五)——EM算法以及Baum-Welch算法的推导

EM算法推导

序列比对(十六)——Baum-Welch算法估算HMM参数

3️⃣ 多序列比对(2):方法

EM算法

EM 算法

全局比对与局部比对

序列比对(三)——局部联配Smith-Waterman算法

算法文献阅读14：Waterman's new algorith

生物信息学

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读