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

原创：hxj7

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

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

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

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

EM算法

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

其具体推导过程如下：
最终目的是找到最大对数似然对应的参数。
$$\hat{\theta }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_{\theta }\mathrm{l}\mathrm{o}\mathrm{g}P(x|\theta )\text{\hspace{1em}(1)}$$

假设$y$为缺失数据，我们首先可以得到：
$$\mathrm{l}\mathrm{o}\mathrm{g}P(x|\theta )=\mathrm{l}\mathrm{o}\mathrm{g}P(x,y|\theta )-\mathrm{l}\mathrm{o}\mathrm{g}P(y|x,\theta )\text{\hspace{1em}(2)}$$

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

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

为了得到${\theta }^{t+1}$，我们首先变换得到：
$$\begin{array}{rl}\mathrm{l}\mathrm{o}\mathrm{g}P(x|\theta )=& \sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(x,y|\theta )\\ & -\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(y|x,\theta )\end{array}\text{\hspace{1em}(4)}$$

公式（4）可以这样推导得到：
将公式（2）等式两边乘上$P(y|x,{\theta }^t)$，再对所有的$y$求和，得到：
$$\begin{array}{rl}\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(x|\theta )=& \sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(x,y|\theta )\\ & -\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(y|x,\theta )\end{array}\text{\hspace{1em}(4.1)}$$
很容易看出等式左边：
$$\begin{array}{rl}\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(x|\theta )& =\mathrm{l}\mathrm{o}\mathrm{g}P(x|\theta )\sum _{y}P(y|x,{\theta }^t)\\ & =\mathrm{l}\mathrm{o}\mathrm{g}P(x|\theta )\end{array}\text{\hspace{1em}(4.2)}$$
由公式（4.1）和（4.2）可以得到公式（4）。

如果令
$$Q(\theta |{\theta }^t)=\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(x,y|\theta )\text{\hspace{1em}(5)}$$

那么我们可以得到：
$$\begin{array}{rl}& \mathrm{l}\mathrm{o}\mathrm{g}P(x|\theta )-\mathrm{l}\mathrm{o}\mathrm{g}P(x|{\theta }^t)\\ & =Q(\theta |{\theta }^t)-Q({\theta }^t|{\theta }^t)+\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}\frac{P(y|x,{\theta }^t)}{P(y|x,\theta )}\end{array}\text{\hspace{1em}(6)}$$

公式（6）的推导也很容易：
$$\begin{array}{rl}& \mathrm{l}\mathrm{o}\mathrm{g}P(x|\theta )-\mathrm{l}\mathrm{o}\mathrm{g}P(x|{\theta }^t)\\ & =[\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(x,y|\theta )-\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(y|x,\theta )]\\ & -[\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(x,y|{\theta }^t)-\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(y|x,{\theta }^t)]\\ & =[Q(\theta |{\theta }^t)-\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(y|x,\theta )]\\ & -[Q({\theta }^t|{\theta }^t)-\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(y|x,{\theta }^t)]\\ & =[Q(\theta |{\theta }^t)-Q({\theta }^t|{\theta }^t)]\\ & +[\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(y|x,{\theta }^t)-\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}P(y|x,\theta )]\\ & =Q(\theta |{\theta }^t)-Q({\theta }^t|{\theta }^t)+\sum _{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}\frac{P(y|x,{\theta }^t)}{P(y|x,\theta )}\end{array}\text{\hspace{1em}(6.1)}$$

其中${\sum }_{y}P(y|x,{\theta }^t)\mathrm{l}\mathrm{o}\mathrm{g}\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{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_{\theta }Q(\theta |{\theta }^t)\text{\hspace{1em}(7)}$$

我们定义概率分布$P(x)$相对于概率分布$Q(x)$的相对熵为：
$$H(P||Q)=\sum _{i}P(x_i)\mathrm{l}\mathrm{o}\mathrm{g}\frac{P(x_i)}{Q(x_i)}\text{\hspace{1em}(7.1)}$$
那么$H(P||Q)\ge 0$，证明如下：
我们知道：
$$\mathrm{l}\mathrm{o}\mathrm{g}x\le x-1,whenx>0\text{\hspace{1em}(7.2)}$$
那么：
$$-\mathrm{l}\mathrm{o}\mathrm{g}x\ge 1-x,whenx>0\text{\hspace{1em}(7.3)}$$
所以：
$$\begin{array}{rl}\mathrm{l}\mathrm{o}\mathrm{g}\frac{P(x_i)}{Q(x_i)}& =-\mathrm{l}\mathrm{o}\mathrm{g}\frac{Q(x_i)}{P(x_i)}\\ & \ge 1-\frac{Q(x_i)}{P(x_i)}\end{array}\text{\hspace{1em}(7.4)}$$
因而：
$$\begin{array}{rl}H(P||Q)& =\sum _{i}P(x_i)\mathrm{l}\mathrm{o}\mathrm{g}\frac{P(x_i)}{Q(x_i)}\\ & \ge \sum _{i}P(x_i)(1-\frac{Q(x_i)}{P(x_i)})\\ & =\sum _{i}P(x_i)-\sum _{i}Q(x_i)\\ & =1-1\\ & =0\end{array}\text{\hspace{1em}(7.5)}$$
等号成立当且仅当对所有的$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{l}\mathrm{o}\mathrm{g}P(x,\pi |\theta )\text{\hspace{1em}(8)}$$

我们知道，HMM模型中，对于一条给定的路径$\pi$，模型中的参数（如转移概率$a_{kl}$、发射概率$e_k(b)$等）都会在计算$\mathrm{l}\mathrm{o}\mathrm{g}P(x,\pi |\theta )$的式子中出现多次。假设$a_{kl}$出现的次数为$A_{kl}$，而$e_k(b)$出现的次数为$E_k(b)$，那么对于路径$\pi$：
$$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 )}\text{\hspace{1em}(9)}$$

所以：
$$\mathrm{l}\mathrm{o}\mathrm{g}P(x,\pi |\theta )=\sum _{k=0}^M\sum _{l=1}^M{A}_{kl}(\pi )\mathrm{l}\mathrm{o}\mathrm{g}{a}_{kl}+\sum _{k=1}^M\sum _b{E}_k(b,\pi )\mathrm{l}\mathrm{o}\mathrm{g}{e}_k(b)\text{\hspace{1em}(10)}$$

由此，公式（8）变成：
$$\begin{array}{rl}Q(\theta |{\theta }^t)& =\sum _{\pi }P(\pi |x,{\theta }^t)\times \\ & [\sum _{k=0}^M\sum _{l=1}^M{A}_{kl}(\pi )\mathrm{l}\mathrm{o}\mathrm{g}{a}_{kl}+\sum _{k=1}^M\sum _b{E}_k(b,\pi )\mathrm{l}\mathrm{o}\mathrm{g}{e}_k(b)]\end{array}\text{\hspace{1em}(11)}$$

改变公式（11）的求和顺序，我们可以得到：
$$Q(\theta |{\theta }^t)=\sum _{k=0}^M\sum _{l=1}^M{A}_{kl}\mathrm{l}\mathrm{o}\mathrm{g}{a}_{kl}+\sum _{k=1}^M\sum _b{E}_k(b)\mathrm{l}\mathrm{o}\mathrm{g}{e}_k(b)\text{\hspace{1em}(12)}$$

其中：
$$\begin{gathered} 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 )\text{\hspace{1em}(12.1)} \end{gathered}$$

公式（12）可以这样推导得到：
公式（11）先对$\pi$求和，可以得到：
$$\begin{array}{rl}Q(\theta |{\theta }^t)& =\sum _{k=0}^M\sum _{l=1}^M\sum _{\pi }P(\pi |x,{\theta }^t)A_{kl}(\pi )\mathrm{l}\mathrm{o}\mathrm{g}{a}_{kl}\\ & +\sum _{k=1}^M\sum _b\sum _{\pi }P(\pi |x,{\theta }^t)E_k(b,\pi )\mathrm{l}\mathrm{o}\mathrm{g}{e}_k(b)\\ & =\sum _{k=0}^M\sum _{l=1}^M\mathrm{l}\mathrm{o}\mathrm{g}{a}_{kl}\sum _{\pi }P(\pi |x,{\theta }^t)A_{kl}(\pi )\\ & +\sum _{k=1}^M\sum _b\mathrm{l}\mathrm{o}\mathrm{g}{e}_k(b)\sum _{\pi }P(\pi |x,{\theta }^t)E_k(b,\pi )\\ & =\sum _{k=0}^M\sum _{l=1}^M\mathrm{l}\mathrm{o}\mathrm{g}{a}_{kl}{A}_{kl}+\sum _{k=1}^M\sum _b\mathrm{l}\mathrm{o}\mathrm{g}{e}_k(b)E_k(b)\end{array}\text{\hspace{1em}(12.2)}$$

现在关键就是怎么样取$a_{kl}$以及$e_k(b)$的值，可以最大化$Q(\theta |{\theta }^t)$。我们令
$$a_{kl}^0=\frac{A_{kl}}{{\sum }_{l^{\prime }}{A}_{k{l}^{\prime }}}{e}_k^0(b)=\frac{E_k(b)}{{\sum }_{b^{\prime }}{E}_k(b^{\prime })}\text{\hspace{1em}(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{array}{rl}& Q^0(\theta |{\theta }^t)-Q(\theta |{\theta }^t)\\ & =[\sum _{k=0}^M\sum _{l=1}^M{A}_{kl}\mathrm{l}\mathrm{o}\mathrm{g}{a}_{kl}^0+\sum _{k=1}^M\sum _b{E}_k(b)\mathrm{l}\mathrm{o}\mathrm{g}{e}_k^0(b)]\\ & -[\sum _{k=0}^M\sum _{l=1}^M{A}_{kl}\mathrm{l}\mathrm{o}\mathrm{g}{a}_{kl}+\sum _{k=1}^M\sum _b{E}_k(b)\mathrm{l}\mathrm{o}\mathrm{g}{e}_k(b)]\\ & =[\sum _{k=0}^M\sum _{l=1}^M{A}_{kl}\mathrm{l}\mathrm{o}\mathrm{g}{a}_{kl}^0-\sum _{k=0}^M\sum _{l=1}^M{A}_{kl}\mathrm{l}\mathrm{o}\mathrm{g}{a}_{kl}]\\ & +[\sum _{k=1}^M\sum _b{E}_k(b)\mathrm{l}\mathrm{o}\mathrm{g}{e}_k^0(b)-\sum _{k=1}^M\sum _b{E}_k(b)\mathrm{l}\mathrm{o}\mathrm{g}{e}_k(b)]\\ & =\sum _{k=0}^M\sum _{l=1}^M{A}_{kl}\mathrm{l}\mathrm{o}\mathrm{g}\frac{a_{kl}^0}{a_{kl}}+\sum _{k=1}^M\sum _b{E}_k(b)\mathrm{l}\mathrm{o}\mathrm{g}\frac{e_k^0(b)}{e_k(b)}\\ & =\sum _{k=0}^M\left[\sum _{l^{\prime }}{A}_{k{l}^{\prime }}\right]\sum _{l=1}^M\frac{A_{kl}}{{\sum }_{l^{\prime }}{A}_{k{l}^{\prime }}}\mathrm{l}\mathrm{o}\mathrm{g}\frac{a_{kl}^0}{a_{kl}}\\ & +\sum _{k=1}^M[\sum _{b^{\prime }}{E}_k(b^{\prime })]\sum _b\frac{E_k(b)}{{\sum }_{b^{\prime }}{E}_k(b^{\prime })}\mathrm{l}\mathrm{o}\mathrm{g}\frac{e_k^0(b)}{e_k(b)}\\ & =\sum _{k=0}^M\left[\sum _{l^{\prime }}{A}_{k{l}^{\prime }}\right]\sum _{l=1}^M{a}_{kl}^0\mathrm{l}\mathrm{o}\mathrm{g}\frac{a_{kl}^0}{a_{kl}}+\sum _{k=1}^M[\sum _{b^{\prime }}{E}_k(b^{\prime })]\sum _b{e}_k^0(b)\mathrm{l}\mathrm{o}\mathrm{g}\frac{e_k^0(b)}{e_k(b)}\\ & \ge 0\end{array}\text{\hspace{1em}(13.1)}$$

最后一个不等式中${\sum }_{l=1}^M{a}_{kl}^0\mathrm{l}\mathrm{o}\mathrm{g}\frac{a_{kl}^0}{a_{kl}}$可以看作$a_{kl}^0$对于$a_{kl}$的相对熵；${\sum }_b{e}_k^0(b)\mathrm{l}\mathrm{o}\mathrm{g}\frac{e_k^0(b)}{e_k(b)}$可以看作$e_k^0(b)$相对于$e_k(b)$的相对熵。所以二者都不小于0。

小结

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

（公众号：生信了）