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

2019-11-14 本文已影响0人生信了

原创：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算法的推导

EM算法

Baum-Welch算法

小结

猜你喜欢

热点阅读