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

原创：hxj7

本文介绍了如何用Baum-Welch算法来估算HMM模型中的概率参数。

Baum-Welch算法应用于HMM的效果

前文《序列比对（15）EM算法以及Baum-Welch算法的推导》介绍了EM算法和Baum-Welch算法的推导过程。Baum-Welch算法是EM算法的一个特例，用来估算HMM模型中的概率参数。其具体步骤如下：

图片引自《生物序列分析》

本文给出了Baum-Welch算法的C代码，还是以投骰子为例，估算出了转移概率以及发射概率。

具体效果如图：
（下面几张图中的 Real 表示真实的转移概率以及发射概率，而Baum-Welch表示用Baum-Welch算法估算的转移概率以及发射概率。）
首先是当若干条序列总长度为300时：


然后是当若干条序列总长度为30000时：


可以看出总长度为30000时已经很接近真实值了。但是，Baum-Welch算法的结果时一个局部最优值，很依赖初始值的设定。所以，当初始值不同时，也有可能会出现这种结果：


小结一下：

• Baum-Welch算法通过多次迭代来估算HMM模型中的概率参数。
• 本文代码设定了迭代的终止条件：当“归一化后的平均对数似然”的变化小于预先设定的阈值时或者迭代次数超出最大迭代次数时，迭代终止。
• Baum-Welch算法的最终结果非常依赖初始值的设定。
• 本文代码中的初始值是随机值。
• 在计算期望次数时，使用了伪计数。

代码中所用公式及其推导

其中的$A_{kl}$指的是$a_{kl}$在所有训练序列中出现的期望次数，而$E_k(b)$指的是$e_k(b)$在所有训练序列中出现的期望次数。用符号表示就是（其中$x^j$表示第j条符号序列）：
$$\begin{array}{rl}{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 )\end{array}\text{\hspace{1em}(1.1)}$$
(1.2) E k ( b ) = ∑ j ∑ π P ( π j ∣ x j , θ ) E k ( b , π j ) = ∑ j ∑ i P ( π i j = k , x i j = b ∣ x j , θ ) = ∑ j ∑ { i ∣ x i j = b } P ( π i j = k ∣ x j , θ ) \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{1.2} \end{aligned} Ek​(b)​=j∑​π∑​P(πj∣xj,θ)Ek​(b,πj)=j∑​i∑​P(πij​=k,xij​=b∣xj,θ)=j∑​{ i∣xij​=b}∑​P(πij​=k∣xj,θ)​(1.2)

我们可以推导出，对某一条序列$x^j$有如下结论：
$$P({\pi }_i=k,{\pi }_{i+1}=l|x,\theta )={\tilde{f}}_k(i)a_{kl}{e}_l(x_{i+1}){\tilde{b}}_l(i+1)\text{\hspace{1em}(2.1)}$$
$$P({\pi }_i=k|x,\theta )={\tilde{f}}_k(i){\tilde{b}}_k(i)s_i\text{\hspace{1em}(2.2)}$$

其中${\tilde{f}}_k(i)$、${\tilde{b}}_k(i)$ 以及 $s_i$ 的定义在前文《序列比对（12）：计算后验概率》中已经给出（下文给出了计算公式）。

公式（2.1）的推导如下：
P ( π i = k , π i + 1 = l ∣ x , θ ) = P ( π i = k , π i + 1 = l , x ∣ θ ) P ( x ∣ θ ) = P ( π i = k , π i + 1 = l , x 1 , . . . , x i , x i + 1 , . . . , x L ∣ θ ) P ( x ∣ θ ) = P ( x 1 , . . . , x i , π i = k ∣ θ ) P ( x i + 1 , . . . , x L , π i + 1 = l ∣ x 1 , . . . , x i , π i = k , θ ) P ( x ∣ θ ) = f k ( i ) P ( x i + 1 , . . . , x L , π i + 1 = l ∣ x 1 , . . . , x i , π i = k , θ ) P ( x ∣ θ ) = f k ( i ) P ( x i + 1 , . . . , x L , π i + 1 = l ∣ π i = k , θ ) P ( x