第二季 · 第 53 章

看不见的生命状态

《算法的故事》第二季 · 第五篇：破解生命的密码——从鲍姆的隐马尔可夫模型，到基因预测、Forward–Backward与Baum–Welch算法

1994年。一台计算机收到了一段来自大肠杆菌的DNA序列。屏幕上密密麻麻地排列着四种字母。

A C G T A A G C C T G A C G T ...

对于计算机而言，它们只是一串字符。但对于生物学家而言，这些字符背后存在更复杂的结构。某些片段可能属于能够编码蛋白质的基因。某些片段位于基因之间。某些位置可能对应基因开始和结束的信号。还有一些区域具有特定的重复模式。问题在于：测序仪能够告诉我们，一个位置上读取到的碱基是什么。却不会在每个碱基旁边自动写上：“这里是基因。”“这里不是基因。”这些标签需要科学家通过实验、已有知识和计算分析进行判断。如果一段DNA长达几十万个字符，逐个位置人工检查显然十分困难。更麻烦的是：一个碱基A，并不只会出现在非编码区域。一个碱基G，也不只会出现在编码区域。相同的碱基可以出现在完全不同的生物学上下文里。因此，计算机不能简单地根据单个字符决定它究竟属于什么区域。它必须结合前后的信息，推断一段连续的隐藏结构。这一年，三位研究者发表了一篇论文：A Hidden Markov Model That Finds Genes in E. coli DNA《一种能够在大肠杆菌DNA中寻找基因的隐马尔可夫模型》。作者是：安德斯·克罗（Anders Krogh）。I. Saira Mian。戴维·豪斯勒（David Haussler）。论文于1994年6月21日提交，11月11日正式发表在《Nucleic Acids Research》。他们构造了一套包含密码子、基因间序列、调控信号以及可能测序错误的概率状态模型。随后，将它用于预测大肠杆菌基因的位置。

训练数据大约包含100万个经过注释的DNA碱基。独立测试数据约有32.5万个碱基。根据论文报告，这套模型准确定位了大约80%的已知基因，并为另外大约10%的基因给出了近似位置。它还指出了一些可能存在新基因或测序错误的位置。 ^{[PubMed][Oxford Academic]}这些数字不是今天所有基因预测工具都能达到的准确率。它们来自1994年那篇论文的具体模型和测试集。但在当时，它们展示了一件非常重要的事情。计算机不仅可以比较两条已经获得的DNA序列。还可以尝试根据一条序列的统计特征，推断那些并没有直接出现在数据中的内部状态。问题由此发生了变化。上一章BLAST问的是：数据库里有什么序列与它相似？这一章HMM问的是：我们看到的这些字符，可能是由怎样的隐藏过程产生的？而使这项研究成为可能的数学方法，诞生于更早的年代。不是1994年。也不是1990年。我们必须把时间重新拨回1966年。

第一幕：1966年，鲍姆与彼得里提出了一个看不见内部状态的数学问题

1966年12月，《The Annals of Mathematical Statistics》发表了一篇研究论文。题目是：Statistical Inference for Probabilistic Functions of Finite State Markov Chains《有限状态马尔可夫链概率函数的统计推断》。作者：Leonard E. Baum（伦纳德·鲍姆）。Ted Petrie（特德·彼得里）。论文刊登于第37卷第6期，第1554—1563页。^[JSTOR]这篇论文的重要背景，是概率论中的一个问题。假设一个系统内部拥有有限数量的状态。这些状态随着时间变化。但观察者不能直接看到它们。只能看到与内部状态有关的随机输出。那么：我们能否根据外部观测，推断系统的概率结构？能否估计状态转移规律？这就是后来隐马尔可夫模型发展的重要数学起点之一。

为什么要叫作马尔可夫？

因为我们需要先理解一种特殊的随机过程。马尔可夫链（Markov Chain）。它的核心条件可以写成：

\boxed{ P(Z_t\mid Z_1,\ldots,Z_{t-1}) = P(Z_t\mid Z_{t-1}) }

其中：

Z_t

代表第

t

个时刻的状态。这项条件说明：如果已经知道当前状态，那么预测下一个状态时，不需要再额外知道更早的全部状态历史。这通常被称为：一阶马尔可夫性质。但不要将它误解为：“现实世界没有历史。”它的真正含义是：在我们定义的模型中，当前状态已经包含预测下一步所需的相关历史信息。如果当前状态没有保留足够的信息，那么这项模型假设就可能不合适。

一个看得见状态的例子

假设一台机器只有两种工作状态：正常。故障。每天，它都可能保持原状态，也可能切换。例如：

今天	明天正常	明天故障
正常	0.9	0.1
故障	0.3	0.7

这就是一个简单的马尔可夫状态转移矩阵。如果我们可以直接查看机器内部诊断结果，那么每天的真实状态是可见的。问题相对简单。

如果状态无法直接观察呢？

现在，假设我们无法查看机器内部。只能听到它发出的声音。正常时：可能发出平稳声音。也可能偶尔发出异常声音。故障时：更容易发出异常声音。但也不是每一次都必然如此。于是，观察到：

平稳、平稳、异常、平稳、异常

并不能直接告诉我们每天机器真正处于什么状态。我们需要根据：状态转移的可能性。以及不同状态产生不同声音的可能性。一起推断隐藏的状态序列。这便是：隐马尔可夫模型（Hidden Markov Model，HMM）。

1970年：四位研究者继续解决参数估计问题

1966年的数学研究并没有结束这个问题。1970年2月，Baum、Petrie、George Soules与Norman Weiss发表：A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains《马尔可夫链概率函数统计分析中的一种最大化方法》。刊登于《The Annals of Mathematical Statistics》第41卷第1期，第164—171页。^{[ScienceOpen]}这项研究是后来通常称作：Baum–Welch算法的重要历史基础。需要纠正一个常见的过度简化说法：不能把1970年这篇四人署名论文，写成“只有Baum和Welch两个人发表的论文”。实际上，Baum–Welch这一名称来自一系列相关研究与后来的算法命名，1970年这篇重要论文的署名作者是：Baum、Petrie、Soules、Weiss。此外，HMM的统计估计理论，也不是到了1970年才突然从无到有。1966年的研究已经具有重要地位。

当时他们不是为了基因预测而研究HMM

1966年与1970年的这些概率论文，重点是随机过程和统计推断。它们不是专门为了：寻找大肠杆菌基因。预测人类外显子。或者分析蛋白质家族。而写的。后来，研究者发现这些数学方法能够解决生物序列中的重要问题。这正是基础数学工具的一项深刻价值。一个模型可以在某个领域被发展出来，又在几十年后成为另一个学科的重要研究方法。

第二幕：从通信到生命——为什么第49章的维特比算法还不够？

在第49章，我们已经学习了：Viterbi Algorithm（维特比算法）。它可以在一个具有有限状态的模型中，寻找最佳完整状态路径。对于HMM而言，它常用于寻找：

\boxed{ \arg\max_{z_{1:T}}P(z_{1:T}\mid x_{1:T}) }

也就是：给定观测序列后，最可能的一条完整隐藏状态路径。但如果我们想真正理解隐马尔可夫模型，还需要解决另外两个问题。

三个问题，其实完全不同

设观测序列：

X=x_1x_2\cdots x_T

模型参数：

\theta

那么：

问题	数学目标	典型算法
这段观测在模型下有多可能？	$P(X\mid\theta)$	Forward
哪一条隐藏状态路径最可能？	$\arg\max_ZP(Z\mid X,\theta)$	Viterbi
不知道模型参数时，怎样从数据估计？	$\max_\theta P(X\mid\theta)$	Baum–Welch

图53-1：同一段DNA，对应三个不同的计算问题

观察到的DNA

AATATGCGGCTATA

Forward把所有可能隐藏路径的概率贡献加起来，计算整段观测的总概率。

Viterbi只寻找概率最高的一条完整隐藏状态路径。

Baum–Welch反复估计隐藏状态的概率分布，再利用期望计数更新模型参数。

这三项任务有一条非常清晰的关系。Forward负责求和。Viterbi负责取最大值。Baum–Welch则利用Forward与Backward计算出的后验期望，进行参数更新。1989年，AT&T贝尔实验室研究者Lawrence R. Rabiner发表了一篇著名教程：A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition。论文将HMM研究整理成评估、解码和参数学习三大类问题，并系统介绍相关算法。^{[数字对象标识符]}这篇教程对后来HMM的教学和工程应用产生了深远影响。但必须记住：1989年Rabiner并没有发明1966年Baum与Petrie已经研究的概率模型，也没有发明1967年的维特比算法。他的重大贡献之一，是将相关理论与实现方法整理成非常系统的教程，并介绍它们在语音识别中的应用。接下来，我们要把这三项任务放在同一个DNA实验中，真正理解它们的差别。

第三幕：一个只有两种隐藏状态的DNA世界

我们先建立一套最小的DNA隐马尔可夫模型。它只有两种隐藏状态。第一种：

N=\text{AT偏好区域}

第二种：

C=\text{GC偏好区域}

我们选择这两个名称，是为了清楚地描述数学模型，而不是预先断言它们对应真实的非编码区与编码区。在某些基因组里，不同功能区域确实可能具有不同的碱基组成。但真实编码区与非编码区不能仅凭GC含量作出可靠区分。本章的教学模型，只模拟一种具有不同碱基组成的隐藏区域结构。

第一项参数：初始状态概率

假设DNA序列的第一个位置，处于两种隐藏状态的概率相同。

P(Z_1=N)=0.5

P(Z_1=C)=0.5

写成向量：

\boxed{ \pi= \begin{bmatrix} 0.5&0.5 \end{bmatrix} }

第二项参数：状态转移概率

假设：AT偏好区域有90%的概率，在下一个位置仍属于AT偏好区域。GC偏好区域也有90%的概率继续保持原状态。只有10%的概率切换。因此：

\boxed{ A= \begin{bmatrix} 0.9&0.1\\ 0.1&0.9 \end{bmatrix} }

第一行表示从N出发。第二行表示从C出发。这个矩阵有两个特点。每个元素都在0与1之间。每一行之和都是1。因为从某个状态出发，下一步一定要进入某个合法状态。

第三项参数：发射概率

现在，真正体现DNA统计差异的参数出现了。假设处于N状态时，更容易观察到A或T。处于C状态时，更容易观察到G或C。我们设计：

隐藏状态	A	C	G	T
N：AT偏好	0.40	0.10	0.10	0.40
C：GC偏好	0.10	0.40	0.40	0.10

这便是：发射概率（Emission Probability）。记作：

b_j(x)=P(X_t=x\mid Z_t=j)

例如：

P(G\mid C)=0.4

P(G\mid N)=0.1

也就是说，看到G时，GC偏好状态相对于AT偏好状态更容易解释这项观测。但千万不要将这句话理解成：看到G，就一定处于C状态。因为：

P(G\mid N)=0.1

并不等于零。

两个状态与四个DNA字符之间的关系

图53-2：两状态DNA-HMM

查看 SVG 配图原件

53-005-44e281d8af26.svg · 点击查看大图

这个模型怎样生成一段DNA？

假设初始状态选择N。那么，按照N状态的发射概率，随机生成一个字符。例如：A。接着，按照状态转移矩阵决定下一步。可能继续N。也可能进入C。如果进入C，那么下一个字符更有机会是G或者C。重复这一过程，便可以产生一段DNA样式的字符序列。但如果我们只拿到最终生成的字符，却看不到内部状态选择过程，问题就反过来了。我们需要推断：这些观测最有可能来自怎样的一串隐藏状态？这便是HMM解码问题。

第四幕：一段只有三个字符的DNA，背后竟然有八条可能历史

现在，我们使用刚才的模型，观察三个DNA字符：

\boxed{AGG}

它们非常短。但如果每一个位置都可能属于N或C，那么隐藏状态路径的数量就是：

2^3=8

全部候选路径为：

NNN
NNC
NCN
NCC
CNN
CNC
CCN
CCC

这八条路径都可能产生：

AGG

但概率不同。

计算第一条路径：NNN

首先，初始状态是N：

P(Z_1=N)=0.5

第一个字符是A：

P(A\mid N)=0.4

第二步仍然为N：

P(N\mid N)=0.9

第二个字符是G：

P(G\mid N)=0.1

第三步仍然为N：

P(N\mid N)=0.9

第三个字符是G：

P(G\mid N)=0.1

因此，这条隐藏路径与观测共同发生的概率是：

P(AGG,NNN) =0.5\times0.4\times0.9\times0.1\times0.9\times0.1

得到：

\boxed{P(AGG,NNN)=0.00162}

再计算：CCC

第一步处于C：

0.5

但C状态发射A的概率只有：

0.1

因此，它一开始看起来并不理想。不过，后面连续两个字符都是G。而C状态发射G的概率为：

0.4

如果状态保持为C，每次状态转移还可以获得：

0.9

的概率。所以：

P(AGG,CCC) =0.5\times0.1\times0.9\times0.4\times0.9\times0.4

得到：

\boxed{P(AGG,CCC)=0.00648}

这比NNN的联合概率更大。虽然第一个字符A并不符合C状态的主要偏好，但后面的连续G以及较高的状态保持概率，使整条CCC路径获得了更高的联合概率。

八条路径的完整比较

图53-3：AGG对应的八条隐藏路径

隐藏路径	联合概率
CCC	0.006480
NCC	0.002880
NNN	0.001620
NNC	0.000720
CCN	0.000180
NCN	0.000080
CNN	0.000045
CNC	0.000020

模型参数为π=(0.5,0.5)，状态保持概率0.9，N状态A/T各0.4，C状态C/G各0.4。全部为精确模型计算结果，并非真实基因组的实验频率。

最可能的一条完整隐藏路径是：

\boxed{CCC}

其联合概率：

0.00648

但请注意：这还没有告诉我们：观察到AGG这件事，在整个模型下发生的总概率是多少。因为除了CCC，其他七条路径也都有可能产生相同的AGG。要得到观测总概率，我们必须将八条路径的联合概率相加。

P(AGG) = \sum_{Z\in\{N,C\}^3}P(AGG,Z)

得到：

\boxed{P(AGG)=0.012025}

这便是Forward算法要计算的目标。而刚才找到CCC，则是Viterbi算法的目标。一个求和，一个取最大值。这项区别，将决定后面所有的公式。

第五幕：Forward算法——为什么不必真的列出八条、八百条或八亿条路径？

对于三个字符：

AGG

我们可以枚举八条隐藏路径。但如果观察序列长达1000个字符，而且有两个隐藏状态，那么可能路径数量为：

2^{1000}

根本不可能全部枚举。于是，我们需要再次使用动态规划。

定义Forward变量

令：

\alpha_t(j) = P(x_1,\ldots,x_t,Z_t=j)

它表示：前

t

个观测已经产生，并且第

t

个隐藏状态为

j

的联合概率。注意：它包含了所有能在第

t

个时刻抵达状态

j

的历史路径。不是只包含其中概率最大的一条。

第一步：初始化

对于第一个观测A：

\alpha_1(N) = P(Z_1=N)P(A\mid N)

所以：

\alpha_1(N)=0.5\times0.4=0.20

同理：

\alpha_1(C)=0.5\times0.1=0.05

第二步：计算第二个观测G

现在，我们希望得到：

\alpha_2(N)

第二个位置处于N，有两种可能的历史来源。上一时刻也是N。或者上一时刻是C。因此：

\alpha_2(N) = [\alpha_1(N)P(N\mid N) +\alpha_1(C)P(N\mid C)] P(G\mid N)

代入：

=(0.20\times0.9+0.05\times0.1)\times0.1

得到：

\boxed{\alpha_2(N)=0.0185}

再计算：

\alpha_2(C) = (0.20\times0.1+0.05\times0.9)\times0.4

得到：

\boxed{\alpha_2(C)=0.026}

第三步：计算第三个观测G

同样的递推：

\alpha_3(N) = (0.0185\times0.9+0.026\times0.1)\times0.1

得到：

\boxed{\alpha_3(N)=0.001925}

而：

\alpha_3(C) = (0.0185\times0.1+0.026\times0.9)\times0.4

得到：

\boxed{\alpha_3(C)=0.0101}

最后，将两个终点状态相加

P(AGG) = \alpha_3(N)+\alpha_3(C)

所以：

0.001925+0.0101=0.012025

这与刚才枚举全部八条路径的结果完全相同。

Forward的核心递推公式

对于一般HMM：

\boxed{ \alpha_t(j) = b_j(x_t) \sum_i \alpha_{t-1}(i)a_{ij} }

其中：

a_{ij}=P(Z_t=j\mid Z_{t-1}=i)

最终：

\boxed{ P(X)=\sum_j\alpha_T(j) }

它只需要反复计算每个时刻、每个状态的联合概率。如果有：

K

个隐藏状态。观测长度：

T

那么，朴素实现的时间复杂度为：

\boxed{O(TK^2)}

而不是：

O(K^T)

这里的突破，与维特比有什么异同？

维特比每次计算：

\max

只保留抵达当前状态的最佳历史。Forward每次计算：

\sum

把全部抵达当前状态的历史概率合并。它们使用同一张状态转移图。但计算目标完全不同。

算法	状态递推中的关键运算	最终得到
Forward	求和	全部隐藏路径贡献的观测概率
Viterbi	最大值	最可能的单条完整隐藏路径

这项差别非常重要。因为：

\boxed{ \max_Z P(X,Z) \neq \sum_ZP(X,Z) }

两者不应该被混淆。

第六幕：Backward算法——为什么后面发生的事，可以改变我们对前面状态的判断？

在刚才的Forward算法中，计算方向是：从左到右。从第一个DNA字符开始，不断向后处理。但如果我们已经观察到整条序列，是否可以反过来利用后面的信息，帮助判断较早位置的隐藏状态？答案是：可以。这便是：Backward Algorithm（后向算法）。

定义Backward变量

令：

\beta_t(i) = P(x_{t+1},\ldots,x_T\mid Z_t=i)

它表示：假设第

t

个位置处于状态

i

，那么后面剩余的观测序列出现的概率是多少？注意，Backward并不重复计算当前位置的发射概率。它从：

t+1

开始，处理后续观测。

从末尾开始

假设我们已经到达最后一个观测位置：

T

后面没有任何字符。因此，在没有额外终止状态约束的基本HMM中：

\boxed{\beta_T(i)=1}

随后，向前递推。

\boxed{ \beta_t(i) = \sum_j a_{ij}b_j(x_{t+1})\beta_{t+1}(j) }

这与Forward的结构很相似。但方向相反。

仍然使用AGG

对于：

X=AGG

我们已经知道：

P(X)=0.012025

Backward从最后一位开始。

\beta_3(N)=\beta_3(C)=1

然后计算第二个位置。如果第二个位置属于N，那么第三个位置可能是N，也可能是C。所以：

\beta_2(N)= 0.9\times0.1+0.1\times0.4

即：

\boxed{\beta_2(N)=0.13}

如果第二个位置属于C：

\beta_2(C)= 0.1\times0.1+0.9\times0.4

即：

\boxed{\beta_2(C)=0.37}

继续向前：

\boxed{\beta_1(N)=0.0265}

\boxed{\beta_1(C)=0.1345}

图53-4：AGG的Forward与Backward结果

位置	观测	α(N)	α(C)	β(N)	β(C)
1	A	.200000	.050000	.0265	.1345
2	G	.018500	.026000	.1300	.3700
3	G	.001925	.010100	1.0000	1.0000

这里使用未归一化的Forward和Backward概率，以便与手算对应。长序列实际编程时，应使用缩放或对数域计算，避免数值下溢。

如何将两个方向的信息合并？

我们真正希望知道的，不是单独的：

\alpha_t(i)

或者：

\beta_t(i)

而是：已经看到整段DNA后，第 $t$ 个位置属于状态 $i$ 的概率。这就是：

\gamma_t(i)=P(Z_t=i\mid X)

利用贝叶斯公式和前后向概率：

\boxed{ \gamma_t(i) = \frac{\alpha_t(i)\beta_t(i)} {P(X)} }

例如，第一个位置属于N的后验概率：

\gamma_1(N)= \frac{0.2\times0.0265}{0.012025}

得到：

\boxed{\gamma_1(N)\approx0.4407}

而：

\boxed{\gamma_1(C)\approx0.5593}

这是不是有些反直觉？第一个字符是：A。但在看完后面连续两个G后，模型认为第一个位置属于C状态的后验概率，反而略高于属于N状态。为什么？因为C状态拥有较高的保持概率。如果整段一直处于C，虽然第一个A不太典型，但后面两个G都比较符合C状态。而频繁切换隐藏状态，也需要付出较低的转移概率代价。因此，后面的观测能够改变我们对前面隐藏状态的判断。

整段AGG的逐位置后验概率

位置	观测	P(N \| AGG)	P(C \| AGG)
1	A	0.4407	0.5593
2	G	0.2000	0.8000
3	G	0.1601	0.8399

所以，三个位置里：C状态的后验概率都超过0.5。而上一幕的维特比最优路径：

CCC

也恰好一致。但这只是这个具体例子的一项结果。一般情况下：逐位置后验概率最大的状态序列，可能与维特比找到的最优完整路径不同。它们优化的是不同目标。

第七幕：概率最大的一条路径，与每个位置最可能的状态，为什么可能不一样？

这项区别尤其容易出现在较长的序列里。假设有两种候选完整路径：

NNNNCCCC

以及：

NNNCCCCC

它们只在一个位置不同。但在概率模型里，可能还有许多其他路径。例如：

NNCCCCCC
NCCCCCCC
NNNNNCCC

每条路径都拥有不同的概率。维特比算法寻找：

\boxed{ Z^*=\arg\max_ZP(Z\mid X) }

它只关心：哪一条完整路径的联合后验概率最大？而Forward–Backward算法能够计算：

\boxed{ P(Z_t=j\mid X) }

它把所有在第

t

个位置处于状态

j

的路径的概率贡献汇总起来。因此，某个位置的后验概率：

P(Z_t=C\mid X)

并不只来自维特比选中的那一条路径。它来自所有在该位置处于C的可能路径。

什么时候应该使用Viterbi？

如果我们希望获得：一条内部逻辑连贯、整体概率最高的状态路径。例如，在一个已经定义了合法基因结构转移规则的模型里，输出一组一致的隐藏状态注释。那么，Viterbi非常自然。

什么时候应该使用Forward–Backward？

如果我们更关心：某个位置属于某状态的置信程度。某个状态边界是否不确定。或者希望进一步进行概率参数估计。那么，Forward–Backward提供了更加丰富的信息。

为什么生物学家尤其需要这种区分？

因为基因边界并不总是非常明确。一段DNA可能存在多个可能的结构解释。如果程序只给出：

CCCCCCCCNNNNNN

研究者可能误以为：第8个位置与第9个位置之间的边界已经被确定。但如果后验概率显示：第7位属于C的概率是0.58。第8位是0.53。第9位是0.47。那说明：这个边界附近仍然存在很大不确定性。而如果概率是：0.999。0.998。0.002。那么，模型对边界附近的判断就更加集中。不过，这些数值仍然是在指定模型及参数条件下的后验概率。如果模型忽略了重要的真实生物学机制，后验概率很高也不意味着实际预测一定正确。

第八幕：Baum–Welch——如果我们连状态转移概率都不知道，怎样让机器自己估计？

到现在，我们已经解决两个问题。Forward能够计算观测序列的总概率。Viterbi能够寻找最优完整状态路径。Forward–Backward还能够计算逐位置的状态后验概率。但还有一个问题。我们之前直接规定：

P(N\rightarrow N)=0.9

P(C\rightarrow C)=0.9

也直接规定：

P(A\mid N)=0.4

P(G\mid C)=0.4

可是，真正分析一份新的生物数据时，研究者未必已经知道这些参数。

如果隐藏状态已经全部标注，问题相对容易

假设我们拥有一份训练数据。每个DNA位置，都已经由可靠的外部信息标注了隐藏状态。例如：

DNA:    A A T G G C T
状态:   N N N C C C N

那么，我们可以直接统计：N状态出现多少次A？多少次T？N到C发生了多少次？C到C发生了多少次？然后用频数估计相应概率。例如：

\hat a_{NC} = \frac{\text{N转移到C的次数}} {\text{所有从N出发的转移次数}}

这叫作监督式参数估计的一种基本方式。

但是，如果真实状态根本没有标注呢？

现在，我们只知道：

A A T G G C T

却不知道：

N N N C C C N

是否正确。我们无法直接数出：究竟有多少次N转移到C？究竟有多少个A来自N状态？因为这些都是隐藏信息。这就出现了一个看起来具有循环依赖的问题。要推断隐藏状态，需要知道模型参数。要估计模型参数，又需要知道隐藏状态。那么，该怎样开始？

一项反复逼近的方法

Baum–Welch使用的核心思想是：先给模型参数一个初始估计。然后利用当前参数，计算隐藏状态的后验分布。再利用这些概率，估计各类状态转移与字符发射的期望次数。随后，用期望次数重新估计参数。不断重复。

图53-5：Baum–Welch的参数学习循环

初始化模型先给初始概率、转移概率和发射概率赋予合法数值

E步：估计隐藏结构利用Forward–Backward计算状态与状态转移的后验期望

M步：重新估计参数利用期望计数更新转移概率与发射概率

重新计算数据似然检查当前参数下观测序列的对数似然

继续迭代或停止直到满足指定的收敛或迭代次数条件

E步究竟估计什么？

首先，计算每个位置属于某状态的后验概率：

\gamma_t(i)=P(Z_t=i\mid X)

同时，还需要计算：第

t

个位置处于状态

i

，第

t+1

个位置处于状态

j

的联合后验概率。记作：

\xi_t(i,j) = P(Z_t=i,Z_{t+1}=j\mid X)

其计算公式为：

\boxed{ \xi_t(i,j)= \frac{ \alpha_t(i)a_{ij}b_j(x_{t+1})\beta_{t+1}(j) }{ P(X) } }

于是：

\sum_{t=1}^{T-1}\xi_t(i,j)

可以被解释成：从状态

i

转移到状态

j

的期望次数。注意，不是某条单独维特比路径上的实际计数。而是对所有可能隐藏路径，按照其后验概率加权求出的期望计数。

M步怎样更新转移矩阵？

使用期望次数：

\boxed{ a_{ij}^{\mathrm{new}} = \frac{ \sum_{t=1}^{T-1}\xi_t(i,j) }{ \sum_{t=1}^{T-1}\gamma_t(i) } }

分子是：

i\rightarrow j

的期望次数。分母是：从状态

i

发出某次转移的期望总次数。因此，更新后的每一行仍然可以归一化为1。

发射概率怎样更新？

假设我们想重新估计：N状态生成A的概率。那么，可以将所有观察到A的位置，属于N状态的后验概率加起来。再除以N状态在整段观测中出现的期望次数。一般形式：

\boxed{ b_j^{\mathrm{new}}(v) = \frac{ \sum_{t:x_t=v}\gamma_t(j) }{ \sum_{t=1}^{T}\gamma_t(j) } }

其中：

v\in\{A,C,G,T\}

这相当于用概率权重，对DNA字母进行软计数。

为什么这叫作“软计数”？

因为我们没有强行宣布：第一个位置一定属于N。第二个一定属于C。而是允许：某个位置有70%的概率属于N。30%的概率属于C。那么，在估计发射频率时，这个位置就分别提供0.7与0.3的期望贡献。这种方法比只根据一条维特比路径进行硬计数，保留了更多隐藏状态的不确定性信息。

Baum–Welch与EM算法是什么关系？

现代统计学通常将Baum–Welch理解为：EM（Expectation–Maximization）算法在HMM参数估计中的一种实现。EM包含两步：E-step：根据当前参数，计算隐藏变量的后验期望。M-step：选择能够最大化相应期望完全数据对数似然的参数。对于固定结构的标准离散HMM，前面的转移与发射概率更新，便是相应M步的常见形式。^{[Wisc Pages]}

它是否保证找到全球最佳参数？

不能。在恰当实现且数值计算稳定、精确完成EM步骤的条件下，Baum–Welch的观测数据对数似然通常具有非下降性质。但这不代表它一定找到：

\arg\max_\theta P(X\mid\theta)

所对应的全局最佳参数。因为不同初始化，可能收敛到不同的局部最优或者其他驻点。而且，HMM还可能存在状态标签交换的问题。例如，把模型里的N和C全部互换，只要相关参数也同步交换，观测序列的概率就不会改变。因此，仅凭无监督训练结果，我们甚至不一定能确定哪个内部状态应该被赋予某个具体的生物学名称。Baum–Welch能够从数据中估计概率结构，但不会自动赋予隐藏状态正确的生物学含义。这项区别非常关键。

第九幕：互动实验室——让36个DNA字符显露隐藏状态

现在，我们构造一段人工DNA。它由三个区段组成。前12个字符，以A、T为主。中间12个字符，以G、C为主。最后12个字符，再次以A、T为主。完整序列是：

AATTAATATATA
GCGCGGCCGCGC
TTATAATAATAT

这段序列共有：

36

个字符。它专门为教学构造，不是从真实物种的基因组截取的。我们希望观察：一个只有两种隐藏状态的HMM，能否从碱基组成的变化中，推断出几个不同区域？

实验：Forward、Viterbi与逐位置后验推断

图53-6：36位DNA隐藏状态推断实验台模型只知道碱基组成概率与状态转移规则，不知道我们构造的三个区段边界。观察改变模型参数后，推断结果如何变化。

显示哪一种状态判断？

每个格子的上半部分是可观测DNA，下半部分是推断的隐藏状态。蓝色N代表AT偏好，绿色C代表GC偏好。

第1—12位第13—24位第25—36位

隐藏状态保持概率

92%

这个值越高，模型越倾向于保持一段连续的状态，不轻易发生N与C的切换。

状态对应字母组的发射偏好

80%

例如设为80%，意味着N发射A或T的合计概率为80%，C发射G或C的合计概率为80%。

检查第几个DNA位置18 / 36

当前观测碱基

G

Viterbi路径判断

C · GC偏好

N状态后验概率0.4%

C状态后验概率99.6%

当前路径的状态切换数

2

观测序列的对数似然

-40.333

本实验同时执行带缩放的Forward–Backward与对数域Viterbi。所有结果均在当前两状态教学模型下计算。蓝色与绿色不对应真实基因的确定标签，模型也没有编码起始密码子、终止密码子或剪接结构。

实验一：为什么默认模型能够识别三个区段？

恢复默认参数。观察：第1—12位。第13—24位。第25—36位。由于这段人工DNA的组成特征非常明显，模型通常能够将中间的GC富集区域与两侧AT富集区域区分开来。但它不是因为程序提前知道：第12位之后一定发生切换。而是因为：相应状态的发射概率能够更好地解释那些碱基。同时，状态保持概率倾向于让同类字符组成连续区段。

实验二：为什么保持概率越高，边界越不容易变化？

将状态保持概率从：

0.92

增加到：

0.98

这意味着，模型更加不愿意切换状态。如果某个位置只出现一两个与当前状态不太匹配的字符，模型可能选择继续保持原状态。因为切换需要支付较低的转移概率代价。

实验三：如果碱基偏好不明显会发生什么？

点击：“高度保持、弱偏好”。此时两个状态的发射分布开始变得接近。例如：N与C对于某个字符的概率差距可能很小。于是，观察到：A。或者：G。能够提供的状态区分信息也随之减少。你会看到，许多位置的状态后验概率更加接近：

50\%

而Viterbi路径可能因为较高的保持概率，选择更少的状态切换。这非常重要。当观测数据不足以区分隐藏状态时，算法应当表现出不确定性，而不是凭空制造确定结论。

第十幕：互动实验室——让模型在没有状态标签的情况下，逐轮学习参数

现在，我们进一步提高难度。上一项实验里，所有参数都是我们事先指定的。这一次，我们不把正确状态标签交给训练程序。只提供人工DNA：

AATTAATATATA
GCGCGGCCGCGC
TTATAATAATAT

然后给出一套并不十分准确的初始参数。模型起初只是稍微偏向：N产生A/T。C产生G/C。我们希望通过Baum–Welch训练，观察这些偏好能否从数据中逐渐增强。

初始模型

我们从：

P(N\rightarrow N)=0.78

P(C\rightarrow C)=0.78

开始。N产生A或T的合计概率为：

0.62

C产生G或C的合计概率也为：

0.62

然后执行E步与M步。反复更新。

实验：观察参数怎样逐轮变化

图53-7：Baum–Welch参数学习实验台同一段36字符的人工DNA；初始参数固定。滑动迭代轮次，观察观测对数似然、状态保持概率和发射分布如何改变。

当前完成的EM迭代次数

0 / 8

训练过程中的观测对数似然对数似然是负数；曲线向上意味着观测数据在当前模型下具有更高概率。它不是预测准确率曲线，也不能证明模型具有正确的生物学解释。

当前观测对数似然

-48.303

相对初始值改善

0.000

参数	初始值	当前值
N→N保持概率	78.0%	78.0%
C→C保持概率	78.0%	78.0%
N产生A/T的合计概率	62.0%	62.0%
C产生G/C的合计概率	62.0%	62.0%

当前模型的维特比状态路径

实验使用标准离散HMM的EM期望计数更新，无额外先验或伪计数。对一条短而刻意分段的序列重复拟合，发射概率可能趋近0或1

第三幕：一段DNA，为什么需要两种看不见的状态？

现在，我们暂时离开复杂的真实基因组。构造一个能够亲手计算的最小模型。假设一条人工DNA由两种隐藏状态控制：状态N：AT偏好区域。状态G：GC偏好区域。我们假设，在状态N下，A和T更常见；在状态G下，G和C更常见。但请注意：这只是一个根据碱基组成设计的教学模型。真实DNA的编码区域、非编码区域并不能简单地由AT含量与GC含量二分。高GC区域不必然是基因，低GC区域也不必然是基因。

模型由三类参数组成

第一类是初始状态概率：

\boxed{\pi=(0.5,\ 0.5)}

表示一开始进入N或G的概率各为50%。第二类是状态转移矩阵：

\boxed{ A= \begin{bmatrix} 0.9&0.1\\ 0.1&0.9 \end{bmatrix} }

行表示当前状态，列表示下一个状态。例如：

P(N\rightarrow N)=0.9

P(N\rightarrow G)=0.1

这意味着模型倾向于在同一状态中连续停留，而不是每读取一个碱基就频繁切换。第三类是发射概率（Emission Probability）。

当前隐藏状态	A	C	G	T
N：AT偏好	0.40	0.10	0.10	0.40
G：GC偏好	0.10	0.40	0.40	0.10

每一行概率之和都是1。

图53-2：最简单的双状态DNA隐马尔可夫模型

查看 SVG 配图原件

53-017-241e9345ab76.svg · 点击查看大图

一段DNA究竟怎样从模型中产生？

想象模型准备输出一段长度为12的DNA。第一步，它按照初始概率选择状态N或G。第二步，根据当前状态的发射概率，随机选择一个碱基。第三步，按照转移矩阵选择下一个状态。第四步，再根据新状态的发射概率选择下一个碱基。如此重复。在某次人工模拟中，隐藏状态可能是：

N N N N G G G G G N N N

观测到的DNA可能是：

A T A A G C G C C T A T

这只是模型可能生成的一种结果，不是唯一结果。即使在状态N中，也有20%的概率输出G或C。所以，看到一个G，不意味着状态一定是G。这便是隐藏状态推断的困难所在。

为什么状态序列要倾向连续？

因为很多真实生物学区域不是每个碱基位置都随机更换功能类别。例如，一段连续的基因间区域可能包含许多碱基；一个编码区也会持续一段长度。所以，如果模型希望寻找连续区域，就需要适当地表达状态持续性。在我们的教学模型里：

P(N\rightarrow N)=0.9

意味着在N状态下，每一步离开N的概率为0.1。在模型假设下，一个连续N状态片段的期望长度是：

E[L_N]=\frac{1}{0.1}=10

个碱基位置。同样：

E[L_G]=10

这来自几何分布。但真实基因、外显子或内含子的长度分布往往更加复杂，因此真正的基因预测模型需要更精细的状态结构与长度建模。现在，我们已经构造了：初始状态概率。转移概率。发射概率。以及一段观察到的DNA。接下来，就可以提出第一项数学问题：如果我们只知道DNA是AGG，这段观测在当前模型下出现的总概率是多少？

第四幕：Forward算法——不是选一条路，而是把所有可能的路加起来

观察序列：

\boxed{X=AGG}

只有三个碱基。每个位置可以属于两种隐藏状态。所以，全部可能的隐藏状态路径有：

2^3=8

条。例如：

N N N
N N G
N G N
N G G
G N N
G N G
G G N
G G G

我们现在不是寻找哪条路径最好。而是要问：这八条路径加起来，总共给观测AGG带来了多少概率？

先直接计算一条路径

假设隐藏状态为：

N N N

观察为：

A G G

联合概率为：

P(NNN,AGG) = P(N)\,P(A\mid N)\,P(N\mid N)\,P(G\mid N)\,P(N\mid N)\,P(G\mid N)

代入：

=0.5\times0.4\times0.9\times0.1\times0.9\times0.1

得到：

\boxed{0.00162}

这只是其中一条路径的概率。为了得到全部观测的概率，我们需要对所有可能状态路径求和：

\boxed{ P(X)=\sum_Z P(X,Z) }

如果序列长度为：

T=100

两状态路径数量就是：

2^{100}

直接枚举将变得极其昂贵。所以，需要动态规划。

Forward的关键状态

定义：

\boxed{ \alpha_t(j)=P(x_1,\ldots,x_t,Z_t=j) }

它表示：已经观测到前

t

个字符，并且第

t

个隐藏状态为

j

的联合概率。注意，它把所有可能的前驱历史都加在了一起。初始化：

\boxed{ \alpha_1(j)=\pi_jb_j(x_1) }

递推：

\boxed{ \alpha_t(j)= b_j(x_t)\sum_i\alpha_{t-1}(i)a_{ij} }

最终：

\boxed{ P(X)=\sum_j\alpha_T(j) }

第一轮：观察A

初始进入N的概率：

0.5

N发出A的概率：

0.4

因此：

\alpha_1(N)=0.5\times0.4=0.2

同理：

\alpha_1(G)=0.5\times0.1=0.05

第二轮：观察G

考虑第二个状态为N。可能来自前一个N，也可能来自前一个G。所以：

\alpha_2(N) = [0.2\times0.9+0.05\times0.1]\times0.1

得到：

\boxed{\alpha_2(N)=0.0185}

如果第二个状态为G：

\alpha_2(G) = [0.2\times0.1+0.05\times0.9]\times0.4

得到：

\boxed{\alpha_2(G)=0.026}

第三轮：再观察一个G

重复递推：

\alpha_3(N)= [0.0185\times0.9+0.026\times0.1]\times0.1

得到：

\boxed{0.001925}

对于状态G：

\alpha_3(G)= [0.0185\times0.1+0.026\times0.9]\times0.4

得到：

\boxed{0.0101}

现在，把两者相加：

P(AGG)=0.001925+0.0101

最终：

\boxed{P(AGG)=0.012025}

图53-3：Forward逐轮计算结果

时刻	观测	$\alpha_t(N)$	$\alpha_t(G)$
1	A	0.200000	0.050000
2	G	0.018500	0.026000
3	G	0.001925	0.010100

观测序列的总概率

0.012025

这是一项怎样的计算突破？

如果有：

K

个隐藏状态，以及：

T

个时间步。枚举全部隐藏路径的数量是：

K^T

但Forward算法每个时刻只需对各目标状态汇总前驱状态。对于稠密转移矩阵，其时间复杂度为：

\boxed{O(TK^2)}

如果：

K=2

它随序列长度T线性增长。这与上一章序列比对、以及第49章维特比算法的动态规划思想完全一致。大量不同历史共享一个当前状态，就可以在当前状态上合并计算。但Forward合并的是所有路径的概率之和。维特比保留的，则是其中概率最大的一条路径。这项区别非常重要。

第五幕：Backward算法——从未来反过来计算，为什么能够帮助判断过去？

现在，我们已经从左到右计算出了：

P(AGG)=0.012025

但是，假设我们希望知道：第一个A背后究竟更可能是N状态还是G状态？仅仅知道整段观测的总概率，还不够。因为：

\alpha_1(N)=0.2

\alpha_1(G)=0.05

这两个数字只考虑了第一个字符A。没有利用后面观察到的两个G。而我们真正需要的，是同时考虑完整观测序列：

AGG

之后，第一个位置属于N或者G的概率。于是，我们需要第二套递推。Backward Algorithm（后向算法）。

Backward的数学含义

定义：

\boxed{ \beta_t(i) = P(x_{t+1},\ldots,x_T\mid Z_t=i) }

它表示：如果我们已经知道第

t

个时刻处于状态

i

，那么此后剩余的观测序列出现的概率是多少？注意，它不包含当前位置之前的观测。也不包括当前位置

x_t

的发射概率。因为那部分已经由Forward变量负责。

为什么可以从末尾开始？

当我们来到最后一个观测位置：

t=T

后面已经没有任何字符。因此：

\boxed{\beta_T(i)=1}

这里的1并不是“最终状态一定正确”。它是空的后续观测序列，在已知当前状态的条件下，所对应的乘法单位概率。然后，从右向左递推：

\boxed{ \beta_t(i)= \sum_j a_{ij}b_j(x_{t+1})\beta_{t+1}(j) }

它与Forward很相似。只是计算方向相反。Forward把过去的路径合并起来。Backward把未来的路径合并起来。

对AGG计算Backward

最后一个位置：

\beta_3(N)=\beta_3(G)=1

第二个位置之后，还有一个G。从N出发：

\beta_2(N) = 0.9\times0.1+0.1\times0.4

所以：

\boxed{\beta_2(N)=0.13}

从G出发：

\beta_2(G) = 0.1\times0.1+0.9\times0.4

得到：

\boxed{\beta_2(G)=0.37}

继续向前：

\beta_1(N) = 0.9\times0.1\times0.13 + 0.1\times0.4\times0.37

得到：

\boxed{\beta_1(N)=0.0265}

同样：

\boxed{\beta_1(G)=0.1345}

完整结果如下。

位置	观测	$\beta_t(N)$	$\beta_t(G)$
1	A	0.0265	0.1345
2	G	0.1300	0.3700
3	G	1.0000	1.0000

Forward与Backward合起来，能得到什么？

它们可以计算：

\boxed{ \gamma_t(i)=P(Z_t=i\mid X) }

也就是：在已知完整观测序列X的情况下，第t个位置属于状态i的后验概率。根据条件概率，可以得到：

\boxed{ \gamma_t(i)= \frac{\alpha_t(i)\beta_t(i)} {P(X)} }

对于第一个位置：

\gamma_1(N) = \frac{0.2\times0.0265}{0.012025}

得到约：

\boxed{0.4407}

而：

\gamma_1(G) = \frac{0.05\times0.1345}{0.012025}

得到约：

\boxed{0.5593}

一个有趣的现象出现了。如果只观察第一个字符A，我们会更倾向于状态N。因为：

P(A\mid N)=0.4

而：

P(A\mid G)=0.1

可是，当后面连续出现两个G以后，我们重新考虑整个序列。第一个位置属于G状态的后验概率，竟然略微超过了N。这是因为模型具有很强的状态持续倾向。如果整段处于G状态，虽然第一个A有些不寻常，但后面两个G很容易得到解释。相比之下，如果第一个位置属于N，而后面又变成G，就需要支付一次状态切换的概率代价。

这是不是意味着未来真的改变了过去？

不是。第一个碱基在序列生成时对应的真实状态，并没有因为我们后来看到两个G而发生变化。改变的是：我们对过去状态的概率判断。当我们拥有更多观测证据后，可以重新计算过去某个隐藏状态的后验概率。这便是Forward–Backward方法的价值。在语音识别、时间序列分析和生物序列注释中，这种利用整段观测重新评估每个位置隐藏状态的能力，非常重要。Rabiner在1989年的经典教程里，明确区分了这种逐时刻后验状态估计与维特比的完整路径最优估计。^{[Proceedings of the IEEE]}

第六幕：Viterbi与Forward–Backward，为什么会对同一个A给出不同的判断？

现在，我们用完全相同的模型与观测序列：

AGG

执行维特比算法。回忆：

\delta_t(j)

表示：所有抵达状态j并产生前t个观测的路径中，概率最高的那一条路径的联合概率。初始化：

\delta_1(j)=\pi_jb_j(x_1)

递推：

\boxed{ \delta_t(j)= b_j(x_t)\max_i \left[\delta_{t-1}(i)a_{ij}\right] }

看起来和Forward极其相似。区别就在：Forward：

\boxed{\sum}

Viterbi：

\boxed{\max}

AGG的维特比结果

我们按照递推计算，得到：

位置	观测	最佳N路径概率	最佳G路径概率
1	A	0.20000	0.05000
2	G	0.01800	0.01800
3	G	0.00162	0.00648

最终最大值是：

\boxed{0.00648}

从它保存的前驱指针向前回溯，得到一条最佳完整路径：

观测：A G G
状态：G G G

这看起来有些反直觉。第一个字符是A。而G状态产生A的概率只有0.1。为什么整条最优路径仍然选择GGG？因为如果整条路径持续保持G，就不需要支付额外的状态切换概率代价。直接计算：

P(GGG,AGG) = 0.5\times0.1\times0.9\times0.4\times0.9\times0.4

得到：

\boxed{0.00648}

它是所有八条完整隐藏状态路径中联合概率最高的一条。

这里至少有三个不同的概率

图53-4：同一个AGG，三种完全不同的统计量

观测总概率Forward对全部8条路径的概率求和

0.012025

最佳路径联合概率Viterbi只取其中最大的一条

0.006480

第一个位置属于G的后验概率Forward–Backward利用完整AGG判断第一个状态

55.93%

第一项：

P(X)

是观测序列出现的总概率。第二项：

\max_ZP(X,Z)

是最佳完整路径的联合概率。第三项：

P(Z_1=G\mid X)

是某个具体位置的隐藏状态后验概率。它们不是同一个量。

最佳路径占全部观测概率的多少？

我们甚至还可以计算：

P(Z=GGG\mid X=AGG)

根据贝叶斯公式：

=\frac{0.00648}{0.012025}

约为：

\boxed{53.89\%}

这说明：GGG虽然是最佳的单条完整路径，但其他七条路径合起来，仍然贡献了相当一部分概率。因此，不能把维特比找到GGG，解释成：“模型对整条路径具有100%的确定性。”它只是：在指定的完整路径候选中，这一条概率最大。

如果只选择每个位置后验概率最高的状态呢？

那是另一个不同的目标。逐位置后验最大化会选择：

\hat z_t=\arg\max_i\gamma_t(i)

而维特比则求：

\boxed{ \hat Z= \arg\max_ZP(Z\mid X) }

在某些模型和观测中，它们可能给出不同的完整状态序列。原因是，逐位置后验最大化分别优化每个位置的判断。维特比则同时要求整条路径具有最高的联合概率，并严格遵守模型中的合法状态转移。如果某些转移概率为零，独立选出的逐位置最佳状态甚至可能构成一条非法路径。这便是为什么，现代生物序列分析不能只输出某个位置“最可能属于哪个状态”。有时候还需要保证：整个预测片段具有一致的结构。例如，一个基因内部的密码子阅读框，不能在没有相应生物学事件的情况下随意跳转。这需要更复杂的状态设计。

第七幕：Baum–Welch——如果连状态转移概率都不知道，算法又该怎样学习？

前面三项计算，都有一个共同前提。我们已经知道HMM的全部参数。例如：

P(N\rightarrow N)=0.9

P(G\rightarrow G)=0.9

以及：

P(A\mid N)=0.4

但真正的科研工作中，这些数字并不会自动出现在DNA旁边。研究者可能只有一批序列。他们并不知道：不同隐藏状态出现的概率。状态之间的转移频率。以及不同状态偏好哪些碱基。那么，能不能让计算机根据观测数据，反过来估计这些参数？这正是：Baum–Welch算法。

一个看起来会陷入循环的难题

如果我们知道每个位置的真实状态，就容易计算参数。例如，在全部N状态对应的位置里，数一数：A出现多少次。C出现多少次。G出现多少次。T出现多少次。然后估计N状态的发射概率。但现在，真实状态是隐藏的。我们并不知道每个A到底来自N还是G。所以，我们需要先推断隐藏状态。然而，要推断隐藏状态，又需要提前知道模型参数。于是问题似乎陷入一个循环。没有参数，就无法推断状态。不知道状态，就无法估计参数。Baum–Welch提供的办法，是让这两项工作交替进行。

第一步：先猜一套参数

假设我们暂时设定：

P(N\rightarrow N)=0.75

P(G\rightarrow G)=0.75

N状态偏好AT。G状态偏好GC。但这种偏好并不特别强。例如：

P(A\mid N)=P(T\mid N)=0.325

P(C\mid N)=P(G\mid N)=0.175

我们不知道这些参数是否准确。但它们提供了一个可以开始计算的模型。

第二步：根据当前参数，估计隐藏状态的后验概率

我们使用Forward–Backward，计算：

\gamma_t(i)=P(Z_t=i\mid X,\theta)

它告诉我们：在当前模型下，每个位置属于状态i的概率是多少。例如，某个位置观测到A。程序可能计算：

P(Z_t=N\mid X)=0.82

P(Z_t=G\mid X)=0.18

于是，我们不必硬性断言这个A一定来自N。而是把它视为：0.82个N状态的期望贡献。0.18个G状态的期望贡献。这便是：软分配（Soft Assignment）。

第三步：还需要计算状态转移的期望次数

仅有：

\gamma_t(i)

还不够。因为我们还希望更新状态转移矩阵。所以，需要定义：

\boxed{ \xi_t(i,j) = P(Z_t=i,Z_{t+1}=j\mid X,\theta) }

它表示：在完整观测序列已知的条件下，第t个状态为i、下一个状态为j的后验联合概率。根据Forward与Backward，可以计算：

\boxed{ \xi_t(i,j)= \frac{ \alpha_t(i)a_{ij}b_j(x_{t+1})\beta_{t+1}(j) }{ P(X) } }

在实际实现中，为避免数值下溢，我们通常使用经过缩放的等价公式。

第四步：用期望计数更新转移概率

设当前模型共有T个观测。我们希望估计：状态i转移到状态j的概率。可以使用：

\boxed{ a'_{ij} = \frac{ \sum_{t=1}^{T-1}\xi_t(i,j) }{ \sum_{t=1}^{T-1}\gamma_t(i) } }

分子：在全部时间中，从i转移到j的期望次数。分母：在能够发生下一步转移的各时刻，处于状态i的期望次数。因此，分子除以分母，就形成新的转移概率估计。

第五步：用期望计数更新发射概率

对于某个碱基：

k\in\{A,C,G,T\}

我们可以写成：

\boxed{ b'_j(k) = \frac{ \sum_{t=1}^{T}\gamma_t(j)\mathbf 1[x_t=k] }{ \sum_{t=1}^{T}\gamma_t(j) } }

其中：

\mathbf1[x_t=k]

是指示函数。如果当前位置观测字符是k，就取1。否则取0。这项公式的含义很直接。我们不再直接统计硬性判定为G状态的位置。而是把每个位置根据：

\gamma_t(G)

的大小，部分计入G状态的统计量。

第六步：重新计算，再更新

现在，我们已经获得一套新的参数：

\theta'

于是：用新参数重新计算Forward–Backward。重新得到隐藏状态的期望分布。再估计新的转移与发射概率。重复进行。直到：观测数据的对数似然提升很小。或者达到预设的训练轮数。这就是Baum–Welch的基本工作方式。

图53-5：Baum–Welch的两阶段循环

E步骤：Expectation

用当前参数执行Forward–Backward，计算每个状态和状态转移的后验期望。

M步骤：Maximization

利用期望统计量重新估计初始概率、转移矩阵与发射概率。

重新执行E步骤，直到满足停止条件

Baum–Welch与EM算法的关系

在现代机器学习语言中，Baum–Welch是：Expectation–Maximization（EM，期望最大化）在离散HMM上的一个经典应用。它并不是随机尝试各种参数，遇到更高分就保留。它利用隐藏变量的后验分布，计算与模型参数有关的期望统计量，再更新参数。在标准的精确EM条件下，完整执行E步骤和M步骤不会降低训练数据的观测对数似然。但需要明确：对数似然不下降，并不代表一定找到全局最优参数。因为HMM的似然函数可能具有多个局部最优点。初始参数不同，最终可能得到不同结果。而如果两个隐藏状态的标签完全互换，模型还可能产生统计上等价的参数表示。这叫作：标签置换问题（Label Switching）。在无监督学习中，模型自己不会天然知道：哪一个状态应该叫作“编码区”。哪一个状态应该叫作“非编码区”。它学到的是具有某种统计性质的状态。是否可以赋予特定的生物学含义，需要独立证据。

1970年代的一项数学思想，怎样走入1990年代的生命科学？

现在，我们终于能够理解历史之间的联系。1960年代与1970年代的Baum等人，研究隐藏状态随机过程的统计估计。1989年的Rabiner，将HMM的评估、解码和训练方法系统整理。1994年的Krogh、Mian与Haussler，则将这些思想与DNA的密码子、基因间序列和具体生物学信号结合起来。因此，HMM进入基因预测的历史，不是：“生物学家发现一个DNA字符串，于是随手套用机器学习。”而是：概率模型、动态规划、统计参数学习与领域知识经过长期发展，最终在一个具体的生命科学问题上结合起来。

第八幕：互动实验室——让一段DNA的隐藏状态浮出水面

现在，我们把刚才全部的数学推导做成一套可以直接操作的实验。我们将使用一个双状态HMM。N状态倾向产生A和T。G状态倾向产生G和C。你可以修改：DNA序列。状态保持概率。以及AT/GC的偏好强度。程序会同时运行：Forward。Backward。Viterbi。然后画出每个位置属于G状态的后验概率，以及维特比选中的整条最佳状态路径。

互动实验：同一段DNA，两种状态推断

图53-6：双状态DNA隐马尔可夫模型实验台修改下方参数，观察整条序列的隐藏状态推断怎样变化。黄色代表更偏向GC状态，蓝色代表更偏向AT状态。

输入人工DNA序列（最多40位）输入仅接受A、C、G、T；本工具不支持真实基因注释。

状态保持概率90%

N状态发出A或T的总概率80%

DNA长度22 bp

观测GC含量36.4%

观测对数似然 ln P(X)-26.3051

每个位置属于G状态的后验概率

条形越高，Forward–Backward认为该位置属于G状态的后验概率越高。灰色表示尚未在下方进度中显示的位置。维特比解码出的完整状态路径

N状态

G状态

逐位置检查22 / 22

第22位观测

A

P(G状态 | 完整DNA)

3.6%

Viterbi完整路径在此位置N（AT偏好）

计算采用缩放Forward–Backward和对数域Viterbi。Forward–Backward的条形图显示逐位置后验概率；Viterbi色块显示全局最佳路径，二者并不保证每一位都完全一致。参数仅用于人工DNA模拟，不代表真实基因的生物学分布。

实验一：为什么单个G不能决定当前状态？

选择：“三段式DNA”。你会看到中间连续的：

GCGCGCGC

通常更容易被模型判定为G状态。而两端的：

AATATAT

和：

ATATATA

更容易倾向N状态。但模型不会保证每个G都被标注为G，也不会保证每个A都被标注为N。因为当前状态不仅取决于当前字符。还取决于整条序列以及转移概率。

实验二：把状态保持概率调到99%

现在，将：

P(N\rightarrow N)

和：

P(G\rightarrow G)

调高到：

0.99

你会发现，模型更不愿意频繁改变状态。即使中间出现一些与当前状态不太一致的字符，也可能继续维持原来的状态。为什么？因为一次状态切换的概率只有：

0.01

与保持状态的：

0.99

相比，代价非常大。这说明：对隐藏状态持续性的先验假设，会影响最终的分段结果。

实验三：把状态保持概率降低到55%

当：

P(\text{保持})=0.55

时，模型没有那么强的状态持续倾向。它更容易根据局部碱基组成变化而调整状态。结果可能出现更多状态切换。但这并不意味着：状态切换越频繁，结果就越符合真实基因结构。我们必须结合实际生物学区域的长度分布，判断模型是否合理。

实验四：手算AGG

选择：“手算例子AGG”。恢复：状态保持概率90%。AT偏好总概率80%。此时，交互模型应与我们刚才的数学推导一致。整段观测概率：

0.012025

而Viterbi路径：

G G G

第一个位置属于G的后验概率约：

55.93\%

这使互动计算和手工数学推导形成了直接对应。

第九幕：第二个互动实验——看着模型把未知参数一点点学出来

现在，我们不再提前告诉模型准确的参数。而是让它从一段人工DNA中，逐步调整状态转移与发射概率。训练序列是：

AATATATATGCGCGCGCTATATATA

它一共有25个碱基。肉眼可以看出，两端偏向A/T，中间偏向G/C。但模型不会提前获得这些区域的真实标签。它只能根据字符和当前概率参数执行训练。初始状态保持概率设为75%。N状态发出A或T的总概率为65%。G状态发出G或C的总概率也为65%。然后执行六轮Baum–Welch。

实验：观察EM训练过程

图53-7：25碱基人工DNA的Baum–Welch训练过程每轮都完整执行Forward–Backward期望计算与参数重新估计。下方数值是使用本章Python算法实际计算并四舍五入得到的记录。

训练数据 · 未提供隐藏状态真值AATATATATGCGCGCGCTATATATA

EM迭代轮次

0 / 6

EM迭代轮次

当前观测对数似然

-33.4507

N状态保持概率75.00%

G状态保持概率75.00%

N状态发出A/T的概率65.00%

G状态发出G/C的概率65.00%

轮次	ln P(X)	N→N	G→G
0	-33.451	0.750	0.750
1	-29.491	0.831	0.740
2	-26.021	0.890	0.829
3	-24.507	0.929	0.874
4	-24.189	0.936	0.877
5	-24.110	0.937	0.876
6	-24.081	0.937	0.876

训练使用一个两状态、四字母、非零初始参数的离散HMM，并对初始概率、转移概率和发射概率进行无平滑最大似然重新估计。这里显示的是固定实验的真实计算轨迹，而不是保证所有数据都产生相同变化。

第一次训练，发生了什么？

最初的对数似然：

\ln P(X)=-33.45075

第一轮重新估计后：

\ln P(X)=-29.49127

由于：

-29.49127>-33.45075

说明训练数据在更新后模型下的概率提高了。继续迭代：

-26.021

-24.50709

-24.18893

-24.11048

最后：

\boxed{-24.08112}

在这六次更新中，观测对数似然持续上升。

模型学到了什么？

N状态对A/T的偏好不断增强。G状态对G/C的偏好也逐渐增强。同时，状态保持概率被重新调整。这些变化说明，模型正在尝试用两种统计上不同的状态，解释训练序列中存在的组成差异。但这并不意味着它“发现了真正的编码区与非编码区”。因为我们的训练数据只是一串人工构造的字符。它没有真正的基因结构。没有密码子阅读框。也没有启动、终止或剪接信号。

还有一个值得警惕的现象

随着训练推进，N状态发出A/T的概率，已经接近：

100\%

这可能反映当前模型很好地适应了这份短训练数据。但也提醒我们：短样本上的最大似然训练，可能产生非常极端的参数。如果拿它去解释新的、具有不同字符分布的DNA，结果未必可靠。因此，在真实研究中，参数学习往往需要考虑：更多训练数据。参数平滑。先验知识。独立验证集。以及合适的模型复杂度。这也是为什么，真正成熟的生物学HMM不能只依靠几轮EM迭代就宣称已经学会了生命的规律。

第十幕：Python实战——完整实现Forward、Backward、Viterbi与Baum–Welch

现在，我们把本章的四种计算整理成一份完整程序。不使用机器学习框架。也不调用专业HMM库。只使用Python标准库。程序能够执行：Forward概率计算。Backward概率计算。逐位置隐藏状态后验推断。Viterbi最优完整路径搜索。以及Baum–Welch参数学习。另外，我们使用缩放技巧，避免较长序列的概率乘积很快下溢。

第一部分：定义HMM与缩放Forward


import math

DNA = "ACGT"
STATES = ("N", "G")


def make_hmm(stay=0.9, bias=0.8):
    """
    两状态DNA教学模型。

    N状态偏好A/T；
    G状态偏好G/C。

    stay:
        状态保持概率

    bias:
        当前状态偏好碱基组的总概率
    """

    pi = [0.5, 0.5]

    trans = [
        [stay, 1 - stay],
        [1 - stay, stay],
    ]

    emit = []

    for j in range(2):
        probs = {}

        for base in DNA:
            favored = (
                (base in "AT") == (j == 0)
            )

            probs[base] = (
                bias / 2
                if favored
                else (1 - bias) / 2
            )

        emit.append(probs)

    return pi, trans, emit


def forward_scaled(seq, pi, trans, emit):
    """
    缩放版Forward。

    返回：
        alpha: 逐时刻归一化Forward变量
        scales: 每一时刻的缩放因子
        loglik: 观测对数似然
    """

    if not seq:
        raise ValueError("序列不能为空")

    k = len(pi)

    alpha = []
    scales = []

    for t, base in enumerate(seq):
        raw = []

        for j in range(k):
            if t == 0:
                prior = pi[j]
            else:
                prior = sum(
                    alpha[-1][i] * trans[i][j]
                    for i in range(k)
                )

            raw.append(
                prior * emit[j][base]
            )

        c = sum(raw)

        if c <= 0:
            raise ValueError(
                "当前模型无法产生该观测序列"
            )

        alpha.append([
            value / c for value in raw
        ])

        scales.append(c)

    loglik = sum(
        math.log(c)
        for c in scales
    )

    return alpha, scales, loglik

为什么需要缩放？因为如果每个位置产生对应碱基的概率大约是：

0.25

那么一段长度为1000的序列，其观测概率可能非常小。直接在普通浮点数里连续相乘，容易发生数值下溢。所以，我们将每一步的Forward值归一化。并保存每一步的缩放因子：

c_t

最终：

\boxed{ \log P(X)=\sum_{t=1}^{T}\log c_t }

这里的：

c_t

可以解释为在已经观察到前面字符的条件下，当前观测字符的预测概率。在该缩放约定下，这些因子的乘积等于完整观测概率。

第二部分：Backward与后验概率


def backward_scaled(seq, trans, emit, scales):
    """
    与forward_scaled配套的缩放Backward。
    """

    k = len(trans)
    n = len(seq)

    beta = [
        [1.0] * k
        for _ in range(n)
    ]

    # 最后位置后面已经没有观测
    for t in range(n - 2, -1, -1):
        for i in range(k):
            beta[t][i] = sum(
                trans[i][j]
                * emit[j][seq[t + 1]]
                * beta[t + 1][j]
                for j in range(k)
            ) / scales[t + 1]

    return beta


def posteriors(seq, pi, trans, emit):
    """
    返回：
        gamma[t][i]:
            位置t属于状态i的后验概率

        xi[t][i][j]:
            从位置t的状态i转移到t+1状态j
            的后验联合概率

        loglik:
            观测对数似然
    """

    alpha, scales, loglik = forward_scaled(
        seq, pi, trans, emit
    )

    beta = backward_scaled(
        seq, trans, emit, scales
    )

    k = len(pi)

    gamma = []

    for t in range(len(seq)):
        row = [
            alpha[t][i] * beta[t][i]
            for i in range(k)
        ]

        z = sum(row)

        gamma.append([
            value / z
            for value in row
        ])

    xi = []

    for t in range(len(seq) - 1):
        matrix = []

        for i in range(k):
            row = []

            for j in range(k):
                value = (
                    alpha[t][i]
                    * trans[i][j]
                    * emit[j][seq[t + 1]]
                    * beta[t + 1][j]
                    / scales[t + 1]
                )

                row.append(value)

            matrix.append(row)

        # 避免浮点误差导致总和略偏离1
        total = sum(
            sum(row) for row in matrix
        )

        xi.append([
            [v / total for v in row]
            for row in matrix
        ])

    return gamma, xi, loglik

这段程序最值得理解的是：

gamma

和：

xi

的区别。gamma[t][i]表示某一个位置处于状态i的概率。xi[t][i][j]表示某个位置从状态i转移到下一个状态j的联合概率。它们刚好对应Baum–Welch需要的两类统计量。

第三部分：对数域Viterbi


def viterbi(seq, pi, trans, emit):
    """
    对数域Viterbi。

    返回：
        最佳状态路径
        最佳路径的对数联合概率
    """

    if not seq:
        raise ValueError("序列不能为空")

    def ln(p):
        return (
            math.log(p)
            if p > 0
            else float("-inf")
        )

    k = len(pi)

    scores = [
        ln(pi[i]) + ln(emit[i][seq[0]])
        for i in range(k)
    ]

    parents = []

    for base in seq[1:]:
        next_scores = []
        previous_states = []

        for j in range(k):
            candidates = [
                scores[i] + ln(trans[i][j])
                for i in range(k)
            ]

            best_i = max(
                range(k),
                key=lambda i: candidates[i]
            )

            next_scores.append(
                candidates[best_i]
                + ln(emit[j][base])
            )

            previous_states.append(best_i)

        scores = next_scores
        parents.append(previous_states)

    last = max(
        range(k),
        key=lambda i: scores[i]
    )

    best_logprob = scores[last]

    path = [last]

    for previous_states in reversed(parents):
        last = previous_states[last]
        path.append(last)

    path.reverse()

    return path, best_logprob

这里再次出现动态规划与回溯。但是，程序不需要计算：

P(X)

的总和来找到最佳路径。因为：

\arg\max_ZP(Z\mid X) = \arg\max_ZP(X,Z)

对于固定的观测X，分母：

P(X)

对所有候选状态路径都相同。因此，在路径比较阶段可以省略。

第四部分：Baum–Welch完整训练


def baum_welch(
    seq,
    pi,
    trans,
    emit,
    iterations=6,
):
    """
    单条DNA序列上的Baum-Welch学习。

    完全重新估计：
        pi
        trans
        emit

    教学实现不使用平滑与先验。
    对很短的数据可能出现极端参数。
    """

    pi = pi[:]

    trans = [
        row[:] for row in trans
    ]

    emit = [
        row.copy() for row in emit
    ]

    history = []

    for r in range(iterations + 1):

        gamma, xi, loglik = posteriors(
            seq, pi, trans, emit
        )

        history.append({
            "round": r,
            "loglik": loglik,
            "pi": pi[:],
            "trans": [
                row[:] for row in trans
            ],
            "emit": [
                row.copy() for row in emit
            ],
        })

        if r == iterations:
            break

        k = len(pi)

        # 更新初始状态概率
        pi = gamma[0][:]

        # 更新状态转移矩阵
        new_trans = []

        for i in range(k):
            denom = sum(
                gamma[t][i]
                for t in range(len(seq) - 1)
            )

            row = []

            for j in range(k):
                numer = sum(
                    xi[t][i][j]
                    for t in range(len(seq) - 1)
                )

                row.append(numer / denom)

            new_trans.append(row)

        # 更新发射概率
        new_emit = []

        for i in range(k):
            denom = sum(
                gamma[t][i]
                for t in range(len(seq))
            )

            probs = {}

            for base in DNA:
                numer = sum(
                    gamma[t][i]
                    for t, x in enumerate(seq)
                    if x == base
                )

                probs[base] = numer / denom

            new_emit.append(probs)

        trans = new_trans
        emit = new_emit

    return history

第五部分：用AGG验证全部公式


pi, trans, emit = make_hmm(
    stay=0.9,
    bias=0.8,
)

seq = "AGG"

alpha, scales, loglik = forward_scaled(
    seq, pi, trans, emit
)

gamma, xi, loglik2 = posteriors(
    seq, pi, trans, emit
)

path, best_logprob = viterbi(
    seq, pi, trans, emit
)

print("观测:", seq)
print("Forward观测概率:", math.exp(loglik))
print("Viterbi状态路径:", "".join(
    STATES[i] for i in path
))
print("最佳路径联合概率:", math.exp(best_logprob))

for t, base in enumerate(seq):
    print(
        f"第{t+1}位 {base}: "
        f"P(G|完整观测)={gamma[t][1]:.4f}"
    )

assert abs(
    math.exp(loglik) - 0.012025
) < 1e-10

assert path == [1, 1, 1]

assert abs(
    math.exp(best_logprob) - 0.00648
) < 1e-10

print("Forward与Viterbi测试通过")

主要输出：

观测: AGG
Forward观测概率: 0.012025
Viterbi状态路径: GGG
最佳路径联合概率: 0.00648

第1位 A: P(G|完整观测)=0.5593
第2位 G: P(G|完整观测)=0.8000
第3位 G: P(G|完整观测)=0.8399

第六部分：验证Baum–Welch的对数似然变化


training_seq = "AATATATATGCGCGCGCTATATATA"

pi, trans, emit = make_hmm(
    stay=0.75,
    bias=0.65,
)

history = baum_welch(
    training_seq,
    pi,
    trans,
    emit,
    iterations=6,
)

print("轮次  loglik    N->N    G->G")

for record in history:
    a = record["trans"]

    print(
        f"{record['round']:>2} "
        f"{record['loglik']:>9.4f} "
        f"{a[0][0]:>7.4f} "
        f"{a[1][1]:>7.4f}"
    )

# 此实验使用精确EM更新。
# 检查观测对数似然没有下降。
for old, new in zip(
    history[:-1],
    history[1:]
):
    assert (
        new["loglik"]
        >= old["loglik"] - 1e-10
    )

print("Baum-Welch训练测试通过")

它与刚才的训练交互实验使用完全相同的人工DNA与初始化参数。所以，结果应当一致。

HMM算法的复杂度

假设：

T

是观测长度。

K

是隐藏状态数量。对于稠密状态转移矩阵：Forward：

\boxed{O(TK^2)}

Backward：

\boxed{O(TK^2)}

Viterbi：

\boxed{O(TK^2)}

Baum–Welch每次迭代，也需要类似量级的Forward–Backward与统计更新。如果需要：

R

轮训练，那么主要计算量约为：

\boxed{O(RTK^2)}

这里还要考虑参数发射统计与模型结构。如果状态转移是稀疏的，可以只遍历合法边，使计算量与实际转移边数有关。这对大规模生物学HMM尤其重要。但最核心的复杂度突破已经清楚：不需要枚举全部 $K^T$ 条隐藏状态路径，就能计算观测似然、状态后验和最优路径，并进行参数重新估计。

第十一幕：1994年，HMM真正开始尝试阅读大肠杆菌的基因结构

现在，让我们回到本章开头。1994年。Krogh、Mian与Haussler发表了那篇关于大肠杆菌基因预测的论文。但我们已经有能力理解：为什么他们需要使用HMM？以及他们真正的模型，与我们的两状态教学模型有什么不同？

真实基因结构不能只用AT与GC含量区分

我们的教学模型只有：N。G。两个状态。但如果要识别真正的蛋白质编码基因，仅靠两种碱基组成偏好显然不够。因为细菌的蛋白质编码区域，具有许多更加明确的生物学特征。其中非常重要的一项是：密码子（Codon）。

三个碱基，决定一个密码子

在标准遗传密码中：三个连续核苷酸构成一个密码子。例如：

ATG

是DNA编码链中的一种常见起始密码子。但基因并不只依赖一个ATG，就能被确认。研究者还需要考虑：正确的阅读框。终止密码子。周围的调控信号。编码区域内部的密码子使用频率。以及更广泛的生物学背景。不同物种对不同密码子的使用偏好，也可能具有统计差异。因此，一个合理的基因预测模型，不能只记住：“当前区域GC含量比较高。”它还可能需要记住：当前位于一个密码子的第几个位置。

为什么要将编码区域拆成多个状态？

假设DNA是一段：

ATGAAACGT...

如果从第一个位置开始，按照三个字符分组：

ATG AAA CGT ...

就形成一种阅读框。如果从第二个位置开始分组：

TGA AAC GT...

阅读框就变了。同一段DNA，不同阅读框可能对应完全不同的密码子结构。因此，HMM可以使用不同的状态来表达编码区域中相位或密码子位置的变化。例如概念上：

编码第1位 → 编码第2位 → 编码第3位 → 编码第1位

这些状态的转移不是任意的。它们必须遵守当前模型定义的合法阅读框顺序。这便是有限状态模型能够表达生物学约束的一个重要例子。

1994年的论文还考虑了什么？

Krogh等人使用了EcoSeq6数据库，由Kenn Rudd维护的大肠杆菌基因组序列资料。他们的模型不仅考虑编码区域的密码子及其频率。还考虑了基因间区域的序列特征。例如：重复性基因间回文序列。Shine–Dalgarno相关的核糖体结合信号。以及可能由测序错误造成的单碱基插入或删除、阅读框偏移等情况。^{[PMC][Oxford Academic]}这说明，真正的生物学HMM并不只是将一个普通马尔可夫链套在DNA字符上。它需要将大量领域知识转化成：状态。转移限制。发射概率。以及特殊的序列结构。

原始论文的训练数据与测试数据

图53-8：1994年大肠杆菌HMM研究的实验规模

参数估计所用注释序列

约100万

核苷酸

独立测试片段

约32.5万

核苷酸

已知基因位置精确识别比例

约80%

另外近似定位的已知基因

约10%

根据1994年Krogh、Mian、Haussler论文摘要整理。统计对象为论文中已知基因的定位表现，不代表每个碱基的准确率，也不应与现代其他数据集的性能直接比较。

这不是只有一小段人工DNA的演示

我们在前面使用：

AATATATGCGCGCGCATATATA

只有22个碱基。而1994年的研究涉及超过百万核苷酸的注释资料，以及独立测试数据。此外，论文中的参数是依据已有注释DNA等资料进行估计的。不能将这项工作描述成：“研究者只把一百万个没有任何标签的字符交给Baum–Welch，然后模型自动学会了全部生物学知识。”既有基因注释与生物学结构假设，都是模型成功的重要条件。

1994年的结果为什么具有意义？

因为它展示了：在已有生物学知识与统计模型之间建立联系，可以帮助科学家从长DNA序列中预测具有内部一致性的基因结构。而不是每个位置孤立地判断：这个字符是不是G？这个字符是不是C？研究者可以用有限状态结构，表达一整段DNA的可能生成过程。然后让动态规划算法寻找符合模型的高概率路径。但原论文发现的候选新基因与可能测序错误，仍然需要其他证据验证。预测，并不等于实验确认。这也是生物信息学模型必须保持的基本科学边界。

第十二幕：1997年，真核生物的基因结构让HMM面临更复杂的挑战

大肠杆菌属于细菌。但人类等真核生物的基因结构，通常更加复杂。一个常见的蛋白质编码基因，可能包含多个外显子（Exon）。这些外显子之间，还可能存在内含子（Intron）。在RNA加工过程中，内含子可以被剪接去除，外显子按照相应规则连接起来。因此，基因预测不再只是判断：这里像不像一段连续的编码序列？还需要推断：外显子从哪里开始？在哪里结束？内含子位于哪里？有哪些剪接信号？多个外显子能否构成一致的阅读框？这对HMM提出了更高要求。

1997年的GENSCAN

1997年4月25日：克里斯·伯奇（Chris Burge）。与：塞缪尔·卡林（Samuel Karlin）。发表：Prediction of Complete Gene Structures in Human Genomic DNA《人类基因组DNA中完整基因结构的预测》。发表于《Journal of Molecular Biology》第268卷，第78—94页。他们介绍了：GENSCAN。这是一套面向真核基因结构预测的概率模型与程序。^[PubMed]

它考虑了哪些生物学信号？

GENSCAN模型综合描述：转录与翻译信号。剪接供体与受体信号。外显子与内含子的长度分布。基因间区域的组成特征。以及不同GC组成区域之间的差异。程序还能够考虑同一DNA序列中的多个基因，处理部分基因结构，以及DNA两条链上的候选结构。^{[Computer Science | Rice University]}这比我们的两状态模型复杂得多。

为什么普通HMM的状态长度分布可能不够用？

在一个简单HMM里，假设某个状态每一步保持的概率为：

p

那么，该状态持续长度服从几何分布：

\boxed{ P(L=\ell)=(1-p)p^{\ell-1} }

这项分布有一个性质：每继续停留一步，离开该状态的条件概率保持不变。但真实外显子、内含子的长度分布，未必适合用单一几何分布表达。因此，研究者可以使用更加一般化的状态持续长度模型。相关方法常被称为：广义隐马尔可夫模型（Generalized HMM）。或者在某些概率建模框架中，与：半马尔可夫模型（Semi-Markov Model）密切相关。它们能够直接表达更加灵活的状态持续长度概率。

模型越复杂，就一定越好吗？

不一定。如果我们增加：更多状态。更多转移。更多发射参数。以及更复杂的持续时间分布。模型可能更加贴近真实生物学结构。但同时也需要更多数据估计参数。如果样本不足，就可能过拟合。而且，状态空间越大，动态规划计算也可能更加昂贵。因此，模型设计始终需要在：表达能力。计算成本。以及统计可学习性之间权衡。1997年的GENSCAN论文，在当时标准测试集上报告，约75%—80%的外显子能够被准确识别。但这同样属于论文中的特定评估结果，不能直接视为今天所有人类基因预测任务的统一准确率。^[PubMed]

第十三幕：1998年，隐马尔可夫模型又重新回到蛋白质序列比对

到现在，我们主要将HMM用于：给DNA每个位置寻找隐藏区域状态。但这不是HMM在生物信息学中的唯一方向。还有另一条非常重要的道路。它与上一章BLAST、以及第51章的序列比对直接相连。

假设我们已经知道一个蛋白质家族

经过实验与比较，科学家发现：某些蛋白质虽然来自不同物种，但可能拥有共同进化来源。它们的序列长度不一定相同。某些位置高度保守。某些位置则允许大量替换。还有一些成员可能多出或缺少几个残基。我们可以先对这些蛋白质进行多序列比对。从中找出：哪些位置普遍存在某种氨基酸？哪些位置变化较大？哪些位置容易发生插入？哪些位置容易发生删除？然后，将这些信息转换成一个概率模型。这便是：Profile Hidden Markov Model（Profile HMM，谱型隐马尔可夫模型）。

为什么需要Match、Insert与Delete状态？

在Profile HMM中，一种常见结构会包含：Match状态。Insert状态。Delete状态。Match表示某个共识位置上的氨基酸分布。Insert表示相对于共识模型出现的额外氨基酸。Delete表示跳过某个共识位置。需要注意：Delete状态通常是静默状态（Silent State）。它不发射一个实际氨基酸字符。因此，Profile HMM与我们前面那种“每个状态每一步一定发射一个碱基”的简单模型，在状态与观测位置的对应方式上存在差异。

图53-9：一个简化Profile HMM的共识位置

查看 SVG 配图原件

53-027-e8f2240dd42c.svg · 点击查看大图

M为匹配状态，I为插入状态，D为删除状态。此图仅解释三类状态，不是完整可运行的Profile HMM拓扑。

1994年，Krogh团队已经研究蛋白质Profile HMM

在本章开头提到的E. coli基因预测论文之外，Krogh、Michael Brown、I. Saira Mian、Kimmen Sjölander和David Haussler，还在1994年2月发表了：Hidden Markov Models in Computational Biology: Applications to Protein Modeling《计算生物学中的隐马尔可夫模型：蛋白质建模应用》。这项研究是Profile HMM发展的重要早期文献。^[PubMed]因此，HMM进入生物信息学并不是沿着：先基因预测，后蛋白质模型。这样单一的时间顺序发展。两条研究方向在1990年代初期已经相互推进。

1998年：Sean Eddy总结Profile HMM的发展

1998年，研究者：肖恩·埃迪（Sean R. Eddy）。发表：Profile Hidden Markov Models《谱型隐马尔可夫模型》。发表于《Bioinformatics》第14卷第9期，第755—763页。这篇综述系统介绍了Profile HMM如何把多序列比对转换成位置特异性评分模型，用于从数据库中寻找可能具有远缘共同进化来源的蛋白质。^[PubMed]

这与上一章BLAST有什么关系？

BLAST常用一条查询序列或相应评分模型，在数据库中寻找高分局部片段。Profile HMM则能够将一个蛋白质家族中多个位置的统计信息，组织成更复杂的状态模型。于是，模型能够体现：某些位置高度保守。某些位置容易发生替换。某些位置允许插入。某些位置可以删除。这为远缘蛋白质家族识别提供了另一种有力方法。HMMER等工具则将Profile HMM研究发展成实用的软件体系。不过，Profile HMM同样不是在所有数据库搜索任务中都无条件优于BLAST。它们的用途、输入资料、评分模型与计算成本可能不同。

第五篇的算法又一次发生了连接

现在我们可以看出：第51章的序列比对，告诉我们怎样寻找序列对应。第52章的BLAST，告诉我们怎样在大型数据库里快速寻找高分候选。第53章的HMM，进一步告诉我们怎样把序列家族的统计结构，表示成包含隐藏状态的概率模型。这三种技术不是三座彼此孤立的小岛。它们共同推动了1990年代生物信息学的发展。

第十四幕：模型能够找到最可能的隐藏状态，但它真的读懂了生命吗？

到现在为止，我们已经看到一项令人惊讶的能力。只要提供：隐藏状态。转移概率。发射概率。以及观测序列。计算机就可以推断：最可能发生过的状态路径。甚至在某些条件下，利用数据重新估计未知的模型参数。但这还不意味着HMM能够自动理解所有生命结构。

第一项局限：模型看到的统计规律，不一定就是基因

我们在教学模型里设置：N状态偏好A/T。G状态偏好G/C。即使模型成功把一段GC含量较高的区域标记成G，也只能说明：这段区域在当前统计模型下，更符合G状态的发射分布。不能因此认定：这里一定存在蛋白质编码基因。因为真正的基因预测还需要考虑：阅读框。起始与终止信号。剪接结构。密码子使用倾向。以及大量生物学证据。状态名称是研究者赋予的，状态含义必须经过验证。

第二项局限：简单HMM假设当前状态决定发射概率

在我们的模型中：

P(x_t\mid Z_t)

只取决于当前隐藏状态。但真实DNA里的相邻碱基，不一定在给定一个简单状态之后就完全独立。例如：某些二核苷酸组合可能具有特殊的统计倾向。蛋白质编码序列受到密码子三联体结构影响。某些调控信号还具有更长距离的上下文依赖。如果希望表达这些关系，就需要：更细致的隐藏状态。更高阶的条件依赖。或者其他更复杂的概率结构。

第三项局限：模型可能过拟合

刚才25个碱基的训练实验中，N状态对A/T的偏好，最终接近100%。这是一个很好的教学现象。但如果训练数据非常短，这种极端概率也可能是过拟合的表现。一个参数很多的HMM，可以非常好地解释某一份训练数据。却未必能够正确解释新的DNA序列。因此，科学家通常需要：独立测试数据。参数正则化或先验。合理的模型结构。以及对真实生物学标签的验证。

第四项局限：不同参数可能解释同一份数据

假设两个状态分别叫作：N。G。如果我们将两个状态的名称全部交换。并同时交换初始概率、转移矩阵的行列及发射参数。模型仍然可能给出完全相同的观测概率。这再次说明：HMM的隐藏状态标签并不是自然界直接刻在DNA上的名称。研究者必须借助已知标注或外部证据，确定每个统计状态与生物学概念的对应关系。

第五项局限：最可能的路径，并不等于完全确定的真相

假设某个位置：

P(Z_t=N\mid X)=0.51

P(Z_t=G\mid X)=0.49

那么，该位置的状态判断显然存在很大不确定性。如果我们只输出：

就会掩盖这种不确定性。因此，在一些实际分析中，逐位置后验概率与候选路径的置信程度，可能比单纯输出一条最佳路径更加重要。一项负责任的科学推断，不仅应当知道自己的最佳猜测，还应当尽量保留这项猜测背后的不确定性。

第十五幕：从第49章到第53章，四种看似不同的算法其实共享一套计算结构

现在，我们已经拥有足够的知识，观察它们之间更深的数学联系。第49章：维特比译码。第51章：Needleman–Wunsch与Smith–Waterman。第52章：BLAST的种子搜索与局部延伸。第53章：Forward、Backward、Viterbi与Baum–Welch。这些算法在形式上并不完全相同。但其中许多核心计算，都能够归入：有限状态上的递推计算。

为什么Forward使用加法？

因为它希望统计所有可能路径的概率贡献。

\alpha_t(j)= b_j(x_t)\sum_i\alpha_{t-1}(i)a_{ij}

这里的：

\sum

表示不同隐藏历史的概率相加。

为什么Viterbi使用最大值？

因为它只需要其中最佳的一条历史。

\delta_t(j)= b_j(x_t)\max_i[\delta_{t-1}(i)a_{ij}]

于是：

\sum

变成：

\max

并且需要保存最优前驱关系。

为什么序列比对也会出现max？

因为它同样希望找到得分最高的一条合法比对路径。例如：

F(i,j)=\max(\text{对角},\text{上},\text{左})

这与维特比的状态最优化，在动态规划结构上具有密切联系。但序列比对的状态表示和路径评分机制不同，不能认为它们只是同一个程序改了名字。

为什么Baum–Welch又重新使用求和？

因为它不仅关心一条最佳路径。它需要计算隐藏状态和状态转移的后验期望。因此，需要考虑所有相关路径的概率贡献。之后，才能利用期望计数更新参数。

图53-10：三类重要的动态规划计算目标

目标	核心操作	代表算法
累积全部路径概率	$\sum$	Forward–Backward
选择最高分路径	$\max$	Viterbi、序列最优比对
利用隐变量期望更新参数	后验期望 + 最大化	Baum–Welch

这张表实际上包含了整个第五篇的一条重要主线。当我们面对大量无法直接枚举的候选历史时，可以利用状态结构，将计算组织成局部递推。有时，我们要将全部候选概率相加。有时，要选择其中最优的一条。有时，还要利用这些概率反过来更新模型本身。这些方法共同构成了现代概率序列分析的重要基础。

尾声：生命没有把答案写在每一个字母旁边

1966年。Baum与Petrie在研究有限状态马尔可夫链的统计推断。他们关注的是随机过程背后那些无法直接观察的内部状态。1970年。Baum、Petrie、Soules和Weiss继续研究，怎样从观测数据中改进概率模型的参数。随后，隐马尔可夫模型的理论和算法不断发展。1989年。Rabiner发表系统性的教程，让更多研究者能够掌握HMM的概率计算、路径推断与参数学习方法。1994年。Krogh、Mian和Haussler将这些数学思想与大肠杆菌的真实DNA序列结构结合起来。模型不再只是在数学纸面上移动概率。它开始帮助科学家寻找真实DNA中的基因。此后，HMM又进一步进入真核基因预测、蛋白质家族识别和其他生物序列分析任务。一项来自概率论的研究，逐渐成为理解生命分子序列的重要计算工具。

如果我们只能记住本章的四个公式

第一个是Forward：

\boxed{ \alpha_t(j)= b_j(x_t)\sum_i\alpha_{t-1}(i)a_{ij} }

它计算：所有隐藏路径共同产生观测序列的概率贡献。第二个是Backward：

\boxed{ \beta_t(i)= \sum_ja_{ij}b_j(x_{t+1})\beta_{t+1}(j) }

它从后向前考虑剩余观测。第三个是Viterbi：

\boxed{ \delta_t(j)= b_j(x_t)\max_i[\delta_{t-1}(i)a_{ij}] }

它寻找：一条完整的最高概率状态路径。第四个是Baum–Welch的典型转移参数更新：

\boxed{ a'_{ij}= \frac{ \sum_{t=1}^{T-1}\xi_t(i,j) }{ \sum_{t=1}^{T-1}\gamma_t(i) } }

它利用隐藏状态转移的期望计数，重新估计模型参数。这四项计算，共同回答了我们在开篇提出的问题。一段DNA序列，能够告诉我们它包含哪些字符。但如果我们想知道：这些字符背后可能存在怎样的结构？哪些隐藏状态最有可能出现？当前模型能否解释这段观测？模型参数又应该怎样通过数据改进？那么，我们就需要超越单纯的字符串比较。进入：概率推断。

第五篇的前三章

《破解生命的密码》· 第51—53章

第51章

序列比对怎样寻找两条DNA或蛋白质序列之间最有意义的对应关系？

第52章

BLAST怎样在巨大的生物序列数据库中，快速发现值得研究的高分局部片段？

第53章

隐马尔可夫模型怎样利用观测序列推断隐藏结构，并从数据中学习概率模型？

我们已经从序列比较，走向了大规模搜索。又从大规模搜索，进一步走向隐藏状态推断。但生命科学还有一个更困难的问题。如果科学家拿到的，不是一条完整DNA，而是成千上万条破碎的小片段，应该怎么办？这正是下一章的故事。

历史档案：隐马尔可夫模型与生物信息学

年份	研究成果	意义
1966年	Baum与Petrie研究有限状态马尔可夫链的概率函数	HMM统计推断的重要理论起点
1967年	Viterbi发表卷积码译码算法	发展有限状态最佳路径动态规划
1970年	Baum、Petrie、Soules、Weiss发表最大化方法	后来Baum–Welch训练思想的重要基础
1989年	Rabiner发表HMM经典教程	系统整理Forward、Backward、Viterbi和参数学习
1994年2月	Krogh等发表蛋白质建模HMM研究	Profile HMM在计算生物学中的重要早期工作
1994年11月	Krogh、Mian、Haussler发表大肠杆菌基因预测研究	HMM用于实际细菌基因预测
1997年	Burge与Karlin发表GENSCAN	更复杂的真核基因结构概率建模
1998年	Sean Eddy发表Profile HMM综述	总结蛋白质家族概率序列模型与数据库搜索方法

主要参考文献与历史依据

1966年 · HMM的早期统计理论Baum & Petrie — Statistical Inference for Probabilistic Functions of Finite State Markov Chains《The Annals of Mathematical Statistics》第37卷第6期，第1554—1563页。用于核查HMM相关统计推断理论的早期历史。

1970年 · Baum等人的最大化研究A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov ChainsBaum、Petrie、Soules、Weiss四人署名的原始论文，第41卷第1期，第164—171页。

1989年 · Rabiner经典教程A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition（PDF）IEEE正式书目建议重点阅读论文中关于前向后向变量、Viterbi回溯和Baum–Welch重新估计的章节。本章的主要递推公式与该教程一致。

1994年 · 大肠杆菌基因预测HMMKrogh、Mian、Haussler — A Hidden Markov Model That Finds Genes in E. coli DNA完整开放原文。包含真实的编码状态模型、密码子频率、基因间序列建模、独立测试与预测结果。

1997年 · GENSCANBurge & Karlin — Prediction of Complete Gene Structures in Human Genomic DNA研究真核基因结构中的外显子、内含子、剪接信号、区域长度与GC组成，展示概率模型如何与基因结构知识结合。

1998年 · Profile HMMSean R. Eddy — Profile Hidden Markov Models解释多序列比对怎样转化成位置特异性概率模型，以及Profile HMM在蛋白质家族分析中的作用。

**史实与教学重构说明：**本章关于1966年Baum与Petrie、1970年四位作者、1989年Rabiner教程、1994年大肠杆菌基因预测、1997年GENSCAN及1998年Profile HMM的历史叙述，依据上述论文与书目。开头的DNA分析情景是结合当时科研问题的叙事重构，不代表一项具有完整现场记录的具体操作。本章的N/G两状态HMM、AGG手算、25位DNA训练数据、互动图与Python程序，全部是为了教学构造的现代实验。它们不是真实基因预测器，也不代表1994年论文所使用的状态数量、基因结构或具体统计参数。Baum–Welch的观测似然改善属于模型内训练性质，不等同于独立测试准确率提高，更不能自动证明学到的状态具有指定的生物学含义。

下一章预告 · 第54章

《把破碎的生命之书重新拼起来：基因组组装与德布鲁因图》

2001年，一项古老的图论思想，再一次走进生命科学。科学家想知道一个生物的完整DNA序列。但测序设备往往只能一次读取其中的一小段。于是，原本连续的基因组，被转换成大量短小的读段。有些读段之间存在重叠。有些读段可能包含测序错误。某些DNA区域还会重复出现。现在，算法必须根据这些碎片，尝试恢复原来的完整序列。如果我们将每条读段看作一个节点，并根据序列重叠情况建立连接，就会得到一张复杂的图。而如果把短词之间的前后缀关系表示成另一种图结构，问题又会发生变化。2001年8月14日，三位研究者在《Proceedings of the National Academy of Sciences》发表：An Eulerian Path Approach to DNA Fragment Assembly《一种基于欧拉路径的DNA片段组装方法》。作者是：帕维尔·佩夫兹纳（Pavel A. Pevzner）。唐海旭（Haixu Tang）。Michael S. Waterman。他们研究了利用欧拉路径思想，重新组织DNA片段组装问题的方法，并尝试改善重复序列导致的组装困难。 ^[PMC]下一章，我们将回到真实的基因组测序历史，从早期重叠读段组装，一直讲到De Bruijn图与欧拉路径。我们将亲手把：

ATGCGAT

分成多个短小的k-mer。再把这些k-mer转换成图中的边。最后让算法沿着一条欧拉路径，将原本破碎的DNA重新拼接起来。这一章，我们根据完整DNA推断看不见的结构；下一章，我们将反过来，试着从不完整的DNA碎片中恢复整条序列。

——《算法的故事》第二季 · 第五篇《破解生命的密码》 · 第53章完

本章配图

第53章配图