跳转至

第 2 章:概率模型(Probabilistic Models)

章节作者:Mykel J. Kochenderfer 章节定位:理论篇第二章,系统介绍不确定性在概率框架下的表示、推理、参数学习与结构学习四大主题;引入贝叶斯网络(Bayesian network)作为变量间条件独立关系的紧凑表示,并讨论在该表示下进行精确与近似推理的方法,最后讨论从数据中学习网络参数与网络结构的两类问题。

2.1 Representation

作者开宗明义指出,理性的决策需要我们对"不确定性"和"目标"进行推理;本章聚焦于"把不确定性表示为概率分布"这一议题。现实问题往往要求在多个变量之上同时进行分布推断——例如同时考虑电力系统状态、推进器状态、轨迹偏差、通信中断等多个随机变量,因此本章将系统讨论如何构造这些多变量概率模型以及如何用它们做推理。本节先回答"为什么要把不确定性形式化"以及"用什么样的语言来形式化",从 2.1.1 节起逐层展开。整章可以视为一条从"信念度公理化"到"贝叶斯网络表示"再到"推理与学习算法"的递进链——这种递进结构与全书"从基本语言到应用案例"的大尺度结构是同构的,也是我们从"决策理论"过渡到"机器学习"的桥梁。

不确定性可以源自对世界状态的不完全信息。设想我们正在监测一颗距地面数千公里的卫星;该卫星已稳定地通过下行链路向地面传送任务数据与遥测数据数月,但突然之间通信中断。多种事件都能导致通信中断:星载电源系统的故障、星载通信系统的故障,或者地面监测系统本身的失效。在掌握手头信息的情况下,要做出唯一确定的诊断是不可能的——任何一种可能的原因都不能被完全排除,也不能被完全确认。这是"观察不完备"所导致的不确定性,也是工程上"传感器冗余"与"多源信息融合"两大思想的根本动因。

不确定性可以源自对世界状态的不完全信息。设想我们正在监测一颗距地面数千公里的卫星;该卫星已稳定地通过下行链路向地面传送任务数据与遥测数据数月,但突然之间通信中断。多种事件都能导致通信中断:星载电源系统的故障、星载通信系统的故障,或者地面监测系统本身的失效。在掌握手头信息的情况下,要做出唯一确定的诊断是不可能的——任何一种可能的原因都不能被完全排除,也不能被完全确认。这是"观察不完备"所导致的不确定性,也是工程上"传感器冗余"与"多源信息融合"两大思想的根本动因。

不确定性也可以源自我们对未来事件进行预测时面临的实践与理论层面的限制。例如,要准确预测一名人类操作员如何响应决策支持系统给出的建议,需要一个关于人脑工作机制的详细模型——而这样的详细模型在实践上几乎不可能得到;即便是卫星轨迹的预测也并非易事:尽管牛顿力学在原理上允许对卫星轨迹做高精度的预测,但姿控推进器的突发故障仍可能让卫星显著偏离标称轨道,使得长期轨迹预测带上有意义的误差。这是"预测能力受限"所导致的不确定性,与观察不完备所导致的不确定性在形式上同构——它驱动了"鲁棒控制"与"模型预测控制 + 反馈校正"等领域的兴起。

一个稳健的决策系统必须把世界当前状态与未来结果中的不确定性来源都纳入考量。在计算层面要刻画这些不确定性,需要一种形式化的表示——这就是本章余下篇幅要展开的概率模型。本节先讨论单个变量上的概率分布,下一节则讨论如何比较两个命题的"可信度"以及如何把这种比较形式化为概率论的基本公理;2.1.2–2.1.7 节讨论各种具体的分布形式与图模型表示(贝叶斯网络、混合贝叶斯网络、时序模型)。从结构上看,本章剩下的小节构成一个"逐层抽象"的递进:单个变量 → 多变量静态 → 多变量时序 → 推理 → 参数学习 → 结构学习。这一递进结构在工程上对应"先建模、再推算、再学习"的实践工作流。

2.1.1 Degrees of Belief and Probability

在涉及不确定性的问题中,能够比较不同陈述的可信度是基本需求。举例而言,我们希望能够把"卫星上存在电气异常的信念"表达为强于"卫星上存在推进器异常的信念"。令命题 E 表示"卫星上存在电气异常"、命题 T 表示"卫星上存在推进器系统异常",则可记为 E ≻ T;当我们以同等强度持有 E 与 T 时,则记 E ∼ T。在给定某些信息时也比较信念的有用性同样存在:例如我们希望表达"在通信中断这一条件下存在电气异常的可能性高于存在推进器异常的可能性"。若以 C 表示通信中断,则记为 (E | C) ≻ (T | C)。这种"对命题可信度的比较"是后续一切概率运算的语义基础,也是决策理论把"信念"与"效用"作为两个独立维度来构建"期望效用"框架的起点。

作者希望对由 ≻ 与 ∼ 算子所诱导的关系做出一组假设。普遍可比性(universal comparability)要求:在给定 C 的条件下,恰好有 (A | C) ≻ (B | C)、(A | C) ∼ (B | C) 或 (A | C) ≺ (B | C) 之一成立——也就是说,对任意两个条件化命题,我们总能在"更可信"、"同等可信"、"更不可信"三者之间给出唯一判决。这条假设与数学中"任意两个实数之间可以比较大小"的性质同构。传递性(transitivity)则要求:若 (A | D) ≻ (B | D) 且 (B | D) ≻ (C | D),则 (A | D) ≻ (C | D)——这与数值排序中的传递性完全一致。作者援引文献 [1] 指出:普遍可比性与传递性两条假设合起来就足以保证信念度可以用一个实值函数来表示。也就是说,我们可以使用一个函数 P,它严格满足下面两条性质:

\(P(A|C) > P(B|C) \quad \text{当且仅当} \quad (A|C) \succ (B|C)\)

\(P(A|C) = P(B|C) \quad \text{当且仅当} \quad (A|C) \sim (B|C)\)

如果我们再对 P 的形式做一组额外的假设(这些假设在多数情境下都自然成立),则可证明 P 必须满足概率论的基本公理;从而 \(0 \le P(A|B) \le 1\)。当我们在条件 B 下对 A 完全确信时,\(P(A|B) = 1\);当认为 (A | B) 不可能时,\(P(A|B) = 0\)。(A | B) 为真的不确定性则由介于两个极端值之间的数值来表达——0 与 1 之间的实数正是我们用来编码"对 (A | B) 的不同强度信念"的语言。值得指出的是:从"对命题可信度的排序"出发推导概率公理,是 Jaynes 等贝叶斯学派的核心工作;这与 Kolmogorov 那种"从测度论出发直接给出概率公理"形成对照——两条路径在数学上等价,但前者更直接地反映了"概率表达的是信念"这一语义,并直接为贝叶斯公式 (2.3) 提供了认识论上的合法性。这一认识论视角也使得"主观概率"(subjective probability)成为可能——只要两个人的"信念度排序"满足普遍可比性与传递性,他们各自的概率赋值都会满足概率公理;不同主体之间的概率可以不同,但同一主体内部的概率必须自洽。这一观点在贝叶斯决策理论中尤其重要:决策者的"主观先验"与其所掌握的证据被一视同仁地纳入推理,从而把"决策"还原为"信念 + 偏好"两个维度的联合优化。

本书并不对概率论做综合性回顾,但重述两条重要的性质,因为它们在后续章节中会反复出现。第一条是条件概率的定义,它把"在 B 条件下 A 的概率"用"联合概率"和"B 的边缘概率"表达出来:

\[ P(A|B) = \frac{P(A,B)}{P(B)} \tag{2.1} \]

其中 \(P(A,B)\) 表示 A 与 B 同时为真的概率。第二条重要性质是全概率公式(law of total probability):当 B 是一组互斥且穷举的命题时(即这些命题两两互斥、且其中必有一个为真),

\[ P(A|C) = \sum_{B \in \mathcal{B}} P(A|B,C) P(B|C) \tag{2.2} \]

由条件概率的定义可以容易地推出下式成立:

\[ P(A|B) = \frac{P(B|A) P(A)}{P(B)} \tag{2.3} \]

该式即贝叶斯公式(Bayes' rule),在本书后续将扮演重要角色——它把"已知 A 时 B 的概率"和"已知 B 时 A 的概率"联系起来,正是"根据证据更新信念"这一核心操作的数学表达。

2.1.2 Probability Distributions

假设我们有一个二元随机变量 A,其取值可以是 0 或 1。A 上的概率分布描述了 A 取各种值时的概率——具体而言即 \(P(A=0)\)\(P(A=1)\)。在下文中,作者将用小写字母配合上标作为讨论随机变量取值的简记:例如 \(P(a^0)\)\(P(A=0)\) 的简写。分布 \(P(A)\)\(P(a^0)\)\(P(a^1)\) 的取值共同决定,但由于 \(P(a^1) = 1 - P(a^0)\),所以该分布只需一个独立参数 \(P(a^0)\) 即可表示。如果 A 是一个取 n 个可能值的离散随机变量,则 \(P(A)\)\(n-1\) 个独立参数 \(P(a^1), \dots, P(a^{n-1})\) 给出,因为 \(P(a^n) = 1 - (P(a^1) + \dots + P(a^{n-1}))\)

如果 A 是连续随机变量,则对 A 上的分布进行表示要复杂一些:任何具体取值出现的概率都是无穷小的。考虑均匀分布 \(U(0, 10)\),它对 \((0, 10)\) 区间内的所有值赋予相等的概率——从这个分布中抽出的样本恰好等于常数 π 的概率实质上为零。但我们可以为"样本落入某区间"这一事件定义非零概率,例如落入 \((3, 4)\) 区间。如果 \(P(A) = U(0, 10)\),则样本 a 落在 3 与 4 之间的概率为 \(1/10\)

连续变量上的分布可以借由累积分布函数概率密度函数来表达。累积分布函数给出了"取值低于某一阈值"的概率质量。若 \(p(a)\) 是 A 上的概率密度函数,则 \(p(a) da\) 是 A 落入区间 \((a, a+da)\) 内的概率(\(da \to 0\))。累积分布函数 P 可以用概率密度函数表达为:

\[ P(a) = \int_{-\infty}^{a} p(a) \, da \tag{2.4} \]

假设我们想用密度函数 \(p(a)\) 来表示纽约 JFK 机场终端区飞机高度上的分布。我们首先需要选择一种分布形式,然后指定其参数。连续变量的一种常用分布是高斯分布(Gaussian distribution,亦称正态分布)。高斯分布由均值 \(\mu\) 与方差 \(\sigma^2\) 参数化:

\(p(w) = N(w | \mu, \sigma^2) \tag{2.5}\)

我们用 \(N(\mu, \sigma^2)\) 表示具有参数 \(\mu\)\(\sigma^2\) 的高斯分布;用 \(N(w | \mu, \sigma^2)\) 表示按如下定义在 w 处的密度:

\[ N(w | \mu, \sigma^2) = \frac{1}{\sigma} \varphi\!\left(\frac{w - \mu}{\sigma}\right) \tag{2.6} \]

式中的 \(\varphi\) 是标准正态密度函数:

\[ \varphi(x) = \frac{1}{\sqrt{2\pi}} \exp(-x^2/2) \tag{2.7} \]

尽管高斯分布常常由于其"只由两个参数定义、计算方便"而成为便利选择,但它存在一些限制,特别是在表示飞机高度分布时。飞机高度必须为正,但高斯分布却对负高度赋予了非零概率;同时它也对不切实际的高高度赋予了非零概率。补救办法是限制其支撑(support,即赋予非零概率的值域),从而得到截断高斯分布,其密度函数为

\[ N(w | \mu, \sigma^2, a, b) = \frac{\frac{1}{\sigma} \varphi\!\left(\frac{w-\mu}{\sigma}\right)}{\Phi\!\left(\frac{b-\mu}{\sigma}\right) - \Phi\!\left(\frac{a-\mu}{\sigma}\right)} \tag{2.8} \]

其中 w 落在区间 \([a, b]\) 内;\(\Phi\) 是标准正态累积分布函数:

\[ \Phi(x) = \int_{-\infty}^{x} \varphi(x) \, dx \tag{2.9} \]

在飞机高度这一例子里,我们可以把分布截断到 0 ft 以下与 65000 ft 以上。

高斯分布的另一项限制是单峰性(unimodality)——也就是说,分布中存在一个点,密度在该点一侧单调递增、另一侧单调递减。而 JFK 终端区的高度分布并不是单峰的:图 2.1 左上的图显示了根据 2011 年 8 月 1800 万条雷达报告估计出的 2000–10000 ft 区间内高度的分布,可以看到每 1000 ft 出现一个峰,这是由空域结构造成的。显然,对这样的分布使用高斯是不恰当的。

表示多峰连续分布有多种办法。一种办法是把一组单峰分布混合起来。高斯混合模型(GMM)就是不同高斯分布的加权和;其参数包括各高斯成分的参数 \(\mu_1, \sigma_1^2, \dots, \mu_n, \sigma_n^2\) 以及它们的权重 \(\rho_1, \dots, \rho_n\)。其密度为

\[ p(x | \mu_1, \sigma_1^2, \dots, \mu_n, \sigma_n^2, \rho_1, \dots, \rho_n) = \sum_{i=1}^{n} \rho_i N(x | \mu_i, \sigma_i^2) \tag{2.10} \]

且权重之和必须为 1。若用 GMM 来表示高度分布,我们可以让各个高斯成分分别以各峰为中心,并赋予合适的权重。

另一种表示多峰连续分布的方法是离散化(discretization)。例如我们可以每 100 ft 划分一个高度区间,把分布表示为分段均匀密度(如图 2.1 所示)。密度由区间端点决定,每一区间上分配相应的概率质量;图 2.1 同时显示了不同离散化方案对高度分布表示效果的影响:200 ft 区间尚可接受,而 1000 ft 区间则会丢失分布的若干重要特征。

2.1.3 Joint Distributions

在现实问题中表示不确定性的一大挑战是要处理多个变量上的联合分布。暂时假设我们要建模二元变量 A、B、C 上的联合分布。表 2.1 展示了一个例子——该表给出了 A、B、C 的所有 \(2^3 = 8\) 种可能赋值所对应的概率。这八个组合是 (0,0,0)、(0,0,1)、(0,1,0)、(0,1,1)、(1,0,0)、(1,0,1)、(1,1,0)、(1,1,1),每一组合对应一个非负实数 \(P(A, B, C)\)

表 2.1 共有 8 项。由于我们枚举了所有可能赋值(也就是对每一种 A、B、C 的 0/1 组合都给出了对应的概率),表中各项之和为 1——这是联合分布作为概率分布的归一化条件。尽管表中有 8 项,但其中只有 7 项是独立的:若以 \(\theta_i\) 表示表中第 i 行的概率,则只需 \(\theta_1, \dots, \theta_7\) 七个参数即可表示该分布,因为 \(\theta_8 = 1 - (\theta_1 + \dots + \theta_7)\)。换言之,n 元联合分布的自由度比"项数"少一——这是因为概率总和必须为 1 这一约束。每一项本身也必须落在 [0, 1] 区间内,从而对参数 \(\theta_i\) 的可行域施加了一组不等式约束——这组约束在工程实现上意味着"联合分布表"不能任意填写,而必须满足归一化与取值范围两条限制。在工程上,表 2.1 这样的"联合分布表"通常不是手工填写的——而是从历史数据中统计得出,或者从一个"专家先验"派生出——这两种来源各有其不确定性,必须与后续的推理/学习算法协调。值得注意的是,给定一个联合分布,我们可以导出所有"边缘分布"与"条件分布"——也就是说,联合分布是关于一组变量的"完整信息"载体,边缘化与条件化都是它的"投影"操作。例如,对表 2.1 的联合分布求 \(P(A, B)\),只需把 C = 0 与 C = 1 的行分别相加。对联合分布做边缘化是 2.2 节"推理"操作的核心步骤——也就是说,本小节介绍的"联合分布表"既是 2.1 节其他表示方法的概念锚点,也是 2.2 节精确推理算法的起点。从信息论的角度看,联合分布表 2.1 中的每一项都是一个非负实数;这意味着可以用一个 \(2^3\) 维的"概率向量"来表示这张表——这一表示在工程上对应"把联合分布展平为一个长向量"的操作,是后续近似推理算法中常用的数据布局。联合分布表也可视为对"系统全部可能状态"的一个完整列举——这种"枚举式"表示在原理上最简单,但代价是规模随变量数指数增长。

若有 n 个二元变量,则指定其联合分布所需的独立参数最多为 \(2^n - 1\) 个。这种参数数量随变量个数的指数式增长,使得不确定性的表示与概率模型的学习都变得困难。即使 n 仅为 20,参数规模就已达到百万量级;n = 30 时参数规模已超过十亿——这正是后面 2.1.4 节要引入"贝叶斯网络"这种紧凑表示的根本动机。贝叶斯网络通过把联合分布分解为各节点上"只依赖其父节点"的条件分布之积,把参数规模从指数级降到了与图结构相关的多项式级。换言之,n 个二元变量的联合分布可以用一个含"独立参数 = ∑ 父节点数"的多项式规模参数集来表达——这正是从"朴素的联合分布表"到"贝叶斯网络"这一表示跃迁的核心收益。从学习的角度,"参数多"也意味着"需要的数据多"——指数规模的参数需要指数规模的样本来可靠地估计。因此,"联合分布表"这种表示除了概念简单之外几乎没有工程价值,它存在的意义更多是"作为概念锚点"——所有更复杂的表示都可以被视为它的某种紧凑化或结构化版本。

2.1.4 Bayesian Network Representation

贝叶斯网络(Bayesian network)是联合分布的一种紧凑表示。网络的结构以图来表达,图的元素是节点(nodes)有向边(directed edges)。每个节点对应一个随机变量;有向边(有时也称为箭头)连接成对的节点,且图中不允许出现环(这是"有向无环图"这一术语的来源)。这些箭头指示直接的概率关系——若 A → B 是一条边,则 A 的取值会直接影响 B 上的分布。每个节点 \(X_i\) 都关联一个条件概率分布 \(P(X_i | Pa_{X_i})\),其中 \(Pa_{X_i}\)\(X_i\) 在图中的父节点(parents)——也就是那些有边直接指向 \(X_i\) 的节点集。

图 2.2 展示了一个用于卫星监测问题的贝叶斯网络例子,其中涉及五个二元变量。电池失效与太阳能板失效都很少发生,但太阳能板失效的可能性要高于电池失效(这是领域内对"哪种电源故障更常见"的工程经验)。两者中的任何一个失效都可能导致电气系统失效;而电气系统失效还可能由电池或太阳能板以外的原因造成,例如电源管理单元的问题——也就是说,电气系统失效这一节点有两个入度为 1 的父节点(B、S),但它们并不构成对该失效原因的完整枚举。电气系统失效可能造成两类可观测后果:轨迹偏差(用地面望远镜可以观测到)以及通信中断(这会中断向各个地面站下行的遥测与任务数据)。非电气系统原因也能造成轨迹偏差与通信中断——这正是轨迹偏差与通信中断两个节点都"指向电气系统"的原因:它们可以由电气系统失效引起,也可以由其他未在模型中显式表达的故障引起。

五个变量各与五个条件概率分布相关联(如图 2.3 所示)。因为 B 与 S 没有父节点,我们只需要指定 \(P(B)\)\(P(S)\)\(P(B)\) 可由一个独立参数 \(P(b^0)\) 给出,\(P(S)\) 同理由 \(P(s^0)\) 给出——即"电池失效"和"太阳能板失效"两个事件在没有任何外部信息时各自出现的先验概率。

与 E 关联的节点有两个父节点 B 与 S。表 2.2 表示 \(P(E | B, S)\),共有 \(2^3 = 8\) 行;其中仅一半行需要指定,原因是 \(P(e^1 | b, s) = 1 - P(e^0 | b, s)\),其中 b、s 表示 B、S 的任意赋值——也就是说,对每一组 (b, s) 我们只需指定 \(P(e^0 | b, s)\)\(P(e^1 | b, s)\) 由互补关系自动给出。剩下两张条件概率表 \(P(D | E)\)\(P(C | E)\) 各由两个独立参数表示。当变量是二元的,\(P(X | Pa_X)\) 可由 \(2^n\) 个独立参数表示,其中 n 是 X 的父节点数。表 2.2 中 E 的独立参数为 4(其他 4 行由互补关系给出);D 与 C 各自只有一个父节点 E,因此独立参数为 2——这与图 2.3 中标出的具体数值一致。

贝叶斯网络的链式法则(chain rule)说明了如何由各局部的条件概率分布构造联合分布。假设我们有变量 \(X_1, \dots, X_n\),并希望计算所有变量取某特定赋值的概率 \(P(x_1, \dots, x_n)\)。以 \(pa_{x_i}\) 表示 \(X_i\) 的父节点被赋予的特定取值。链式法则说:

\[ P(x_1, \dots, x_n) = \prod_{i=1}^{n} P(x_i | pa_{x_i}) \tag{2.11} \]

在卫星的例子中,假设我们想计算"一切正常"的概率 \(P(b^0, s^0, e^0, d^0, c^0)\)。由链式法则,

\[ P(b^0, s^0, e^0, d^0, c^0) = P(b^0) P(s^0) P(e^0 | b^0, s^0) P(d^0 | e^0) P(c^0 | e^0) \tag{2.12} \]

若要完整地指定五个变量 B、S、E、D、C 上的联合分布,则需要 \(2^5 - 1 = 31\) 个独立参数。而上面这个贝叶斯网络允许我们仅用 \(1 + 1 + 4 + 2 + 2 = 10\) 个独立参数即可指定其联合分布。10 与 31 之间的差异在参数数量上算不上特别显著,但在更大的贝叶斯网络中节省可以非常可观——10 节点时由 1023 降至数十,20 节点时由约百万降至百量级。贝叶斯网络的力量来自其降低指定联合分布所需参数数量的能力——这种降低是通过让网络结构编码"某些条件独立假设"来实现的(见下一小节)。换言之,贝叶斯网络并不是一种新的概率分布类型,而是一种紧凑的因式分解——它把指数规模的联合分布表用多项式规模的"局部条件分布 + 图形结构"来等价表达。

2.1.5 Conditional Independence

贝叶斯网络能以比通常所需更少的独立参数表示联合分布,原因在于其图形结构编码的条件独立(conditional independence)假设。若贝叶斯网络所做出的条件独立假设不成立,则有"未能恰当地建模联合分布"的风险——这一点将在 2.4 节进一步讨论;也就是说,网络结构必须与领域内真实的条件独立关系一致。

我们说变量 A 与 B 独立(independent) 当且仅当 \(P(A, B) = P(A) P(B)\)。A 与 B 独立的断言记作 \(A \perp B\)。由式 (2.1) 可知,\(A \perp B\) 当且仅当 \(P(A) = P(A | B)\)。换言之,关于 B 的信息不会给 A 带来任何额外信息,反之亦然。举例而言,假设卫星上的电池失效 (B) 与太阳能板失效 (S) 独立,那么知道发生了电池失效不会让我们提高或降低对"是否发生了太阳能板失效"的信念。我们可以用两个参数 \(P(b^0)\)\(P(s^0)\) 即可指定 \(P(B, S)\) 联合分布(见表 2.3)。事实上,若有 n 个独立的二元变量,则指定其联合分布只需 n 个独立参数——相比之下,若不能使用独立性假设,则需要 \(2^n - 1\) 个独立参数。

变量 A 与 B 在给定 C 的条件下条件独立(conditionally independent)当且仅当 \(P(A, B | C) = P(A | C) P(B | C)\)。A 与 B 在给定 C 的条件下条件独立的断言记作 \((A \perp B | C)\)。由该定义可以证明 \((A \perp B | C)\) 当且仅当 \(P(A | C) = P(A | B, C)\)。给定 C,关于 B 的信息不会给 A 带来任何额外信息,反之亦然。例如,假设卫星轨迹偏差 (D) 在给定"是否发生了电气系统失效" (E) 的条件下,与"是否发生了通信中断" (C) 条件独立。我们可记为 \((D \perp C | E)\)。若已知发生了电气系统失效,则观察到通信中断这一事实不会影响我们对"是否发生轨迹偏差"的信念。我们对发生轨迹偏差的预期可能有所提高,但那仅是因为我们已知电气系统失效已发生——而非通信中断本身提供了新信息。

我们可以用一组规则判断节点 A 与 B 在给定节点集 C 的条件下是否条件独立。若 \((A \perp B | C)\),则我们说 C d-分隔(d-separate)了 A 与 B(其中的"d"代表"directional",意指分隔的方向性源于有向图的结构)。我们也说 A 与 B 之间的一条路径被 C d-分隔,如果下列任一条件为真:

  1. 路径中包含一串链式节点(chain) \(X \to Y \to Z\),且 Y 在 C 中——也就是说,Y 是这条链上的"中转节点"且被观测到,从而阻断了 X 与 Z 之间的依赖;
  2. 路径中包含一个分叉(fork) \(X \leftarrow Y \to Z\),且 Y 在 C 中——Y 是 X 与 Z 的共同原因,观测到 Y 后 X 与 Z 之间的依赖被阻断;
  3. 路径中包含一个反向分叉(inverted fork,也称 v-结构) \(X \to Y \leftarrow Z\),且 Y 不在 C 中、Y 的后代也不在 C 中——这种情况下,Y 是 X 与 Z 的共同结果;只有当我们观测到 Y 或 Y 的某个后代时,X 与 Z 之间的依赖才会被"激活"。

若 A 与 B 之间的所有路径都被 C d-分隔,则 \((A \perp B | C)\)。"Markov 毯(Markov blanket)"这个术语有时被用来指代那个"能把某节点与所有其他节点 d-分隔开的最小节点集"——在吉布斯采样等算法中,Markov 毯的概念允许我们在更新某节点时只需考虑其局部邻居。

在图 2.2 的网络中,存在唯一的 v-结构 \(B \to E \leftarrow S\)。在缺乏关于 E、D 或 C 的信息时,B 与 S 是独立的。然而一旦给定 E、D 或 C,B 与 S 便不再独立:影响可以从 B 流到 S。例如,若已知发生了电气系统失效,则再知道发生了电池失效就会降低我们对"是否发生了太阳能板失效"的信念——因为电池失效"解释掉了"一部分电气系统失效的可能原因,使得太阳能板失效作为另一原因的可能性相应降低。这种通过 v-结构发生的影响有时被称作"解释消退(explaining away)"——电池失效的出现"解释掉了"电气系统失效的部分原因;在医学诊断、风险归因、工程故障归因等领域,"解释消退"是误诊与误归因的常见根源。

2.1.6 Hybrid Bayesian Networks

本章到目前为止的例子都只涉及二元变量,但贝叶斯网络可以同时包含离散变量与连续变量。同时具有离散与连续变量的贝叶斯网络常被称为混合贝叶斯网络(hybrid Bayesian networks)。图 2.4 展示了一个混合贝叶斯网络的例子,它刻画了飞机特性、其雷达截面积、以及雷达探测目标能力之间的关系。翼展较大的飞机往往具有较大的雷达截面积,其中截面积以相对于一平方米的分贝(dBsm)来度量——这是工程中常用的"对数化"雷达截面积单位。军用飞机有时被专门设计成具有较低的雷达截面积(低于 0 dBsm)以规避探测——这种设计被称为"低可观测性(low observability)"或"隐身(stealth)"设计。截面积较大的目标更可能被探测到,但其他因素也会影响探测(例如目标的相对姿态、当时的天气、对雷达的电子干扰等)。在图 2.4 中,M 与 D 自然为二元变量,W 与 C 自然为连续变量——这一混合性质正是"混合贝叶斯网络"这一名称的由来。

与任何贝叶斯网络一样,我们需要为每个节点指定条件分布。节点 W 没有父节点,因此我们只需为 W 指定一个分布;如 2.1.2 节所讨论的,我们将使用由参数 \(\mu\)\(\sigma^2\) 定义的高斯分布,尽管它可能对负的翼展和过大的翼展赋予小概率——这在工程上可以接受,因为真实飞机翼展的取值范围远离这两个尾部。变量 M 是二元的,所以我们可以仅用一个参数 \(P(m^0)\) 即可定义该分布。

截面积 C 同时依赖于连续变量 W 与二元变量 M。我们暂时忽略其对 M 的依赖,先定义密度 \(p(c | w)\)。定义"给定一个连续变量的另一个连续变量上的分布"的一种常见做法是使用线性高斯分布(linear Gaussian distribution)。例如

\(p(c | w) = N(c | \theta_1 w + \theta_2, \theta_3) \tag{2.13}\)

由上式可见,均值是 w 的线性函数,由参数 \(\theta_1\)\(\theta_2\) 定义;方差由 \(\theta_3\) 定义。因为我们希望"翼展越大、截面积越大",因此 \(\theta_1\) 必须取正。翼展极小的飞机的截面积也会极小,所以 \(\theta_2\) 大致应为 0。参数 \(\theta_3\) 控制 c 与 w 之间线性关系中的方差大小——\(\theta_3\) 越大则同等翼展下截面积的散布越宽。

在现实中,C 同时依赖于 W 与 M。我们可以直接让线性高斯分布中所用的参数依赖于 M:

\[ P(c | w, m) = \begin{cases} N(c | \theta_1 w + \theta_2, \theta_3) & \text{if } m^0 \\ N(c | \theta_4 w + \theta_5, \theta_6) & \text{if } m^1 \end{cases} \tag{2.14} \]

这种分布被称为条件线性高斯(conditional linear Gaussian)。在这个例子中我们需要六个参数来表示 \(p(c | w, m)\)。由于军用飞机被更多地设计为具有比非军用飞机更低的雷达截面积,我们会希望 \(\theta_4\) 小于 \(\theta_1\)——也就是说,在给定翼展的条件下,军用飞机的截面积随翼展增长的速度较慢。

最后我们需要定义条件分布 \(P(D | C)\)。我们想表达"雷达更可能探测到截面积较大的飞机"这一性质。当然,我们可以设置一个阈值 \(\theta\),并令 \(P(d^1 | c) = 0\)(若 \(c < \theta\))、\(P(d^1 | c) = 1\)(否则)。然而这种模型可能对那些"实际本可发生"的探测事件赋予零概率——例如截面积恰好略低于阈值的飞机本应有非零概率被探测到,但硬阈值模型把这种情形完全排除掉了。

与其使用硬阈值来定义 \(P(D | C)\),我们可以用一个软阈值(soft threshold)——在低于阈值时赋予低概率、高于阈值时赋予高概率。表示软阈值的一种方式是使用logit 模型,它产生一条"S"形的 sigmoid 曲线:

\[ P(d^1 | c) = \frac{1}{1 + \exp\!\left(-2 \frac{c - \theta_1}{\theta_2}\right)} \tag{2.15} \]

参数 \(\theta_1\) 决定阈值的位置,\(\theta_2\) 控制概率分布的"软度"或"展宽"——\(\theta_2\) 越小,曲线越陡峭,越接近硬阈值;\(\theta_2\) 越大,曲线越平缓,对阈值附近的概率过渡越不敏感。图 2.5 以实线绘出了 \(\theta_1 = 0\)\(\theta_2 = 1\) 时的 \(P(d^1 | c)\)

logit 模型的替代是probit 模型

\(P(d^1 | c) = \Phi((c - \theta_1)/\theta_2) \tag{2.16}\)

其中 \(\Phi\) 是 2.1.2 节介绍的标准正态累积分布函数。logit 模型与 probit 模型在形态上彼此非常接近(logit 的曲线与"标准差为 \(\pi/\sqrt{3}\) 的正态累积"几乎重合),如图 2.5 所示;二者之间的选择通常由计算便利性或具体应用领域决定。

2.1.7 Temporal Models

时序模型(temporal model)刻画一组变量随时间如何演化。一个简单的时序模型是马尔可夫链(Markov chain),其中时刻 t 的状态记为 \(S_t\)。例如,马尔可夫链可以刻画飞机随时间的位置与速度。图 2.6 展示了一个用于表示马尔可夫链的贝叶斯网络结构。图中只画出了前三个状态,但马尔可夫链可以无限延展——这是因为状态转移的规律对所有时刻都成立。初始分布由 \(P(S_0)\) 给出;条件分布 \(P(S_t | S_{t-1})\) 通常被称为状态转移模型(state transition model)。若状态转移分布不随 t 变化,则该模型被称为平稳(stationary)的——这一条件在多数实际应用中作为简化假设成立。

马尔可夫链中的状态不必是标量。例如,若我们要建模飞机随时间的随机行为,状态可以是向量 \(s = (h, \dot{h})\),其中 h 是飞机的高度、\(\dot{h}\) 是垂直速率。初始分布 \(P(S_0)\) 可由多元高斯分布表示,该多元高斯分布由均值向量 \(\mu\) 与协方差矩阵 \(\Sigma\) 参数化,其密度为

\(p(s) = N(s | \mu, \Sigma) \tag{2.17}\)

其中 \(N(s | \mu, \Sigma)\) 是式 (2.6) 中高斯分布的 k 维推广:

\[ N(s | \mu, \Sigma) = \frac{1}{(2\pi)^{k/2} |\Sigma|^{1/2}} \exp\!\left(-\frac{1}{2} (s - \mu)^\top \Sigma^{-1} (s - \mu)\right) \tag{2.18} \]

对于我们的飞机模型,k 取 2,\(\mu\) 是两元素的向量,\(\Sigma\) 是 2×2 矩阵。

飞机模型的状态转移分布可以由如下线性高斯表示:

\(p(s_t | s_{t-1}) = N(s_t | M s_{t-1} + b, \Sigma) \tag{2.19}\)

均值是上一状态的简单线性函数。若飞机在平均意义上保持直线飞行,那么均值的一个合理选择是

\[ M s_{t-1} + b = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix} s_{t-1} \tag{2.20} \]

即:高度按"上一时刻高度 + 上一时刻垂直速率"推进、垂直速率在一步之内不变。协方差矩阵 \(\Sigma\) 控制施加于飞机高度与垂直速率的随机扰动量——\(\Sigma\) 的对角元代表两个分量上"扰动的方差",非对角元则刻画两分量之间扰动的相关性。

马尔可夫链可以通过加入观测节点来扩展,如图 2.7 所示。时刻 t 的观测记为 \(O_t\)。观测节点以阴影表示,意指该节点上的值是已知的。若状态对应飞机的位置与速度,那么观测可以是对距离与方位角的有噪雷达测量。若状态变量是离散的,则模型称为隐马尔可夫模型(HMM);若状态变量是连续的且条件分布是线性高斯,则模型称为线性动态系统(linear dynamical system)。HMM 在语音识别、自然语言处理、生物信息学中有广泛应用;线性动态系统则与控制理论、信号处理密切相关。

线性动态系统的一个例子是上一飞机模型的扩展:在时刻 t 的真实状态由向量 \((h_t, \dot{h}_t)\) 表示。我们的观测是高度的有噪测量;垂直速率无法直接观测。时刻 t 的观测被建模为来自如下线性高斯分布:

\[ p(o_t | s_t) = N(o_t | \begin{bmatrix} 1 & 0 \end{bmatrix} s_t, \Sigma) \tag{2.21} \]

协方差矩阵 \(\Sigma\) 在此处是一个单元素矩阵,它控制测量噪声。\(\begin{bmatrix} 1 & 0 \end{bmatrix} s_t\) 表示"从状态向量中提取高度分量",意即我们只观测到高度、而垂直速率保持隐变量身份。

涉及多个状态变量的平稳时序模型可以用动态贝叶斯网络(dynamic Bayesian network) 紧凑地表示。动态贝叶斯网络由两个贝叶斯网络组成:一个表示初始分布,另一个表示转移分布。转移分布由一个具有两个时间片的贝叶斯网络表示:第一个时间片表示时刻 t 的变量,第二个时间片表示时刻 t+1 的变量。图 2.8 展示了一个具有四个状态变量的动态贝叶斯网络——通过这种方式,整个时间序列的联合分布被分解为"初始时间片的分布"加上"任意两相邻时间片之间的转移分布",从而把无穷多时间步上的联合分布表示为有限结构。动态贝叶斯网络是后续第 4 章马尔可夫决策过程、第 6 章部分可观测马尔可夫决策过程的图模型先导——把"决策节点"与"奖励节点"加入动态贝叶斯网络,即得到"决策动态贝叶斯网络";而把"决策"与"状态观测不完备"二者结合,即得到 POMDP。

2.2 Inference

上一节说明了如何表示概率分布;本节讨论如何使用这些概率表示来进行推理。推理(inference)指在给定一组观测变量的取值后,确定一个或多个未观测变量上的分布。例如,假设我们想用 2.1.4 节介绍的卫星贝叶斯网络来推断 \(P(B | d^1, c^1)\)。此时 B 是查询变量(query variable),D 与 C 是证据变量(evidence variables),S 与 E 是隐变量(hidden variables)。这一查询在工程上对应"我们观察到卫星轨迹偏差与通信中断同时发生,在这一证据下'卫星存在电池失效'的可能性有多大"——它把观察结果"翻译"为对系统状态的后验信念,是诊断与决策的基础。

在讨论几个"推理能派上用场"的例子之后,本节说明如何利用贝叶斯网络结构所固有的特性来高效地执行推理。具体而言,2.2.1 节讨论如何把推理用于分类任务、2.2.2 节讨论在时序模型中的推理、2.2.3–2.2.4 节讨论精确推理的算法与复杂度、2.2.5 节讨论近似推理方法。这五节构成本节从"问题动机"到"算法实现"的递进叙述。从更高的层次看,本节的核心问题是"在已知网络结构与参数的条件下,如何根据观测计算未观测变量上的后验分布"——这一问题的具体形式会随变量的"是否时序"、"是否离散"、"观测是否完备"等条件而变化,但背后的概率论核心(贝叶斯公式、全概率公式、链式法则)始终不变。推理问题在工程上意义重大——它把"已知的模型 + 当前观测"转化为"未观测变量的后验",是诊断、预测、决策三类应用共同的计算基础。本节的最后一个小节(2.2.5)将看到:当精确推理计算上不可行时,蒙特卡洛方法变分方法这两大类近似推理策略在工程上扮演着关键角色;这一观察也为后续章节(如第 6 章 POMDP 的近似求解)埋下伏笔。读者若熟悉机器学习,可以把本节的精确推理对应于"在封闭形式下可解的后验推断"(如线性回归、共轭先验),把近似推理对应于"在不可解析情形下的后验推断"(如神经网络中的随机梯度 MCMC、变分自编码器)——两者的概率论根基相同,但算法实现方式有显著差异。从软件工程的角度看,推理问题的实现往往涉及"稀疏矩阵运算 + 图遍历 + 缓存"三者的紧密结合;许多贝叶斯网络推理库(如 libDAI、pgmpy、BayesiaLab)的核心数据结构都是对这三者的精细化实现。从认知科学的角度看,人类对"在不确定下做推断"的直观能力远不如对"在确定下做计算"的能力——这正是为什么贝叶斯推理在工程上需要形式化算法支持的根本原因。精确推理与近似推理之间的关系是"理想"与"现实"的妥协——理想上我们想要精确后验,现实上常常不得不接受采样估计。理解这一权衡,是设计实用推理系统的关键。本节的另一个隐含主题是"算法与数据结构的协同"——精确推理的高效实现往往需要"问题结构 + 算法设计 + 数据结构"三者的紧密配合;离开这三者中的任何一环,再聪明的算法也无法在工程上发挥其威力。读者若熟悉数据库系统,可以把精确推理中的"中间表缓存"类比于"查询计划缓存"——两者都通过重用先前计算来减少重复工作量。从应用的角度看,推理问题几乎贯穿所有概率模型的应用场景——一个工程上可用的贝叶斯网络系统,至少应当包含一种精确推理算法与一种近似推理算法,以应对不同规模的问题。

2.2.1 Inference for Classification

推理可以用于分类(classification)任务——给定一组观测或特征,推断类别。例如,假设我们想根据雷达航迹的属性判断一个雷达目标是鸟还是飞机。此时类别为鸟或飞机,观测可包括对速度与航向在整个航迹持续期间内的波动量的测量。大多数飞机比大多数鸟飞得快,但也存在一定重叠,特别是较小、较低性能的飞机——这是"空域管理"中区分"目标是否构成威胁"时需要精细分类的重要原因。迁徙中的鸟通常会保持航向不变,这与执行机动动作的飞机形成对比——因此航向稳定性是一个具有判别力的特征。

朴素贝叶斯模型(naive Bayes model)是分类任务中常用的一种简单概率模型,其结构如图 2.10 所示;图 2.11 给出了一个等价的、更紧凑的表示——板(plate),以圆角矩形表示。板底的 \(i = 1:n\) 表示"变量名下标 i 从 1 重复到 n"——板是描述"重复结构"的一种图模型简记。

在朴素贝叶斯模型中,类别 C 是查询变量,观测特征 \(O_1, \dots, O_n\) 是证据变量。为简洁起见,本书在脚标中会偶尔使用冒号记法:例如 \(O_{1:n}\)\(O_1, \dots, O_n\) 的紧凑写法。朴素贝叶斯之所以称为"朴素",是因为它假设证据变量在给定类别的条件下条件独立。使用 2.1.5 节引入的记号,对所有 \(i \neq j\),可写为 \((O_i \perp O_j | C)\)。如果这些条件独立假设不成立,我们只需在观测特征之间加上必要的有向边即可——也就是说,朴素贝叶斯的"朴素性"在结构上体现为"证据变量之间无边"。

在朴素贝叶斯模型中,我们必须指定先验 \(P(C)\)类条件分布(class-conditional distribution) \(P(O_i | C)\)。在雷达目标分类问题中,先验表示我们在没有任何航迹信息时关于"目标是鸟还是飞机"的信念——例如鸟类迁徙季节与日常时段不同,这一先验应在工程实践中随时间调整。图 2.12 展示了由雷达数据估计出的空速类条件分布的示例。

我们可以通过对式 (2.11) 应用链式法则来推断朴素贝叶斯模型中的联合分布:

\[ P(c, o_{1:n}) = P(c) \prod_{i=1}^{n} P(o_i | c) \tag{2.22} \]

分类任务真正想要的是条件概率 \(P(c | o_1, \dots, o_n)\)。由式 (2.1) 中条件概率的定义可得

\[ P(c | o_{1:n}) = \frac{P(c, o_{1:n})}{P(o_{1:n})} \tag{2.23} \]

分母可借助联合分布与式 (2.2) 中的全概率公式轻易求出:

\[ P(o_{1:n}) = \sum_{c} P(c, o_{1:n}) \tag{2.24} \]

式 (2.23) 中的分母不是 C 的函数,因此可视作常数。于是可写为

\(P(c | o_{1:n}) = \kappa P(c, o_{1:n}) \tag{2.25}\)

其中 \(\kappa\) 是使 \(\sum_c P(c | o_{1:n}) = 1\) 的归一化常数。我们通常略去 \(\kappa\) 而简单写为

\(P(c | o_{1:n}) \propto P(c, o_{1:n}) \tag{2.26}\)

其中符号"\(\propto\)"表示左端"正比于"右端。这一记法上的简化之所以成立,是因为当我们最终需要的是 \(P(c | o_{1:n})\) 作为"对所有 c 取值的一个分布"时,分母 \(\kappa\) 会在归一化步骤中自动恢复。例如,假设由链式法则求得:

\(P(\text{bird, slow, little heading fluctuation}) = 0.03 \tag{2.27}\)

\(P(\text{aircraft, slow, little heading fluctuation}) = 0.01 \tag{2.28}\)

显然这些概率之和并不为 1。若我们想求"在证据条件下目标是鸟"的概率,则做如下计算:

\[ P(\text{bird} | \text{slow, little heading fluctuation}) = \frac{0.03}{0.03 + 0.01} = 0.75 \tag{2.29} \]

借助我们的模型与概率法则,我们已经求出目标是鸟的概率为 0.75、是飞机的概率为 0.25。但对很多应用来说,我们必须就一个具体类别做出选择。确定分类的一种常见方法是选择具有最大后验概率(即考虑了证据后概率最大的)的那个类别。然而,类别选择本身也是一个决策问题,往往应当把误分类的代价纳入考虑。例如,如果我们要用分类器为空中交通管制过滤掉非飞机目标,那么我们可以容忍偶尔漏放一些鸟或其他杂波航迹;但我们必须避免把任何真实的飞机过滤掉,否则可能导致碰撞。在这种情况下,我们大概只会把"后验概率接近 1"的航迹判为鸟。决策问题将在第 3 章讨论。

2.2.2 Inference in Temporal Models

许多重要应用(如语音识别、飞机跟踪、密码分析)都涉及在时序模型中做推理。四个常见的推理任务包括:

  • 滤波(filtering)\(P(S_t | O_{0:t})\)——根据直到当前时刻的所有观测,估计当前状态;
  • 预测(prediction)\(P(S_{t'} | O_{0:t})\),其中 \(t' > t\)——根据当前时刻为止的观测,预测未来某时刻的状态;
  • 平滑(smoothing)\(P(S_{t'} | O_{0:t})\),其中 \(t' < t\)——根据直到当前时刻的观测,估计过去某时刻的状态;
  • 最可能解释(most likely explanation)\(\arg\max_{S_{0:t}} P(S_{0:t} | O_{0:t})\)——找出最有可能产生整个观测序列的状态轨迹。

上面各项都使用了图 2.7 给出的隐马尔可夫结构与记法;变量 t 表示当前时间。滤波与最可能解释是实时应用中常用的两类任务,平滑则多用于离线的事后分析(例如雷达航迹的事后校准),预测则是控制与决策的前置。这四类任务在工程上对应不同的应用需求——例如语音识别关心最可能解释(解码整个发音序列)、雷达跟踪关心滤波(实时更新目标位置)、气候分析关心平滑(事后修正历史估计)、库存管理关心预测(提前备货)。这四类任务也对应不同的算法——例如卡尔曼滤波专用于"线性高斯 HMM"上的滤波、维特比算法专用于最可能解释、前向后向算法专用于平滑。值得指出的是,这四类任务之间的算法并不孤立——许多实际系统在实现时会共享子例程(如对前向变量 \(\alpha_t\) 与后向变量 \(\beta_t\) 的计算)。

为了说明时序模型中的推理,我们将聚焦于"在具有离散状态与离散观测的隐马尔可夫模型中做滤波"。由贝叶斯公式,

\(P(s_t | o_{0:t}) \propto P(o_t | s_t, o_{0:t-1}) P(s_t | o_{0:t-1}) \tag{2.30}\)

表示隐马尔可夫模型的贝叶斯网络结构允许我们做出条件独立假设 \((O_t \perp O_{0:t-1} | S_t)\),这意味着上式中的 \(P(o_t | s_t, o_{0:t-1})\) 等于 \(P(o_t | s_t)\)——也就是说,给定当前状态,过去的所有观测都不再为"当前观测"提供新信息。改写式 (2.30) 并对第二项施加全概率公式,可得

\[ P(s_t | o_{0:t}) \propto P(o_t | s_t) \sum_{s_{t-1}} P(s_t, s_{t-1} | o_{0:t-1}) \tag{2.31} \]

\(P(s_t, s_{t-1} | o_{0:t-1})\) 施加条件概率定义,可得

\[ P(s_t | o_{0:t}) \propto P(o_t | s_t) \sum_{s_{t-1}} P(s_t | s_{t-1}, o_{0:t-1}) P(s_{t-1} | o_{0:t-1}) \tag{2.32} \]

模型的结构保证 \((s_t \perp o_{0:t-1} | s_{t-1})\),因此上式可化简为

\[ P(s_t | o_{0:t}) \propto P(o_t | s_t) \sum_{s_{t-1}} P(s_t | s_{t-1}) P(s_{t-1} | o_{0:t-1}) \tag{2.33} \]

\(P(o_t | s_t)\)\(P(s_t | s_{t-1})\) 由模型直接给出;右端的 \(P(s_{t-1} | o_{0:t-1})\) 提示我们如何随时间推进、随新观测到来递归地更新状态分布——这是一个"用上一时刻的后验 + 当前时刻的观测 → 当前时刻的后验"的两步更新。算法 2.1 给出了这一递归贝叶斯估计(recursive Bayesian estimation)流程的具体过程。记时刻 t 的状态后验分布为 \(b_t\)。为减少脚标数量,算法 2.1 假设状态转移分布 \(P(S_t | S_{t-1})\) 与观测分布 \(P(O_t | S_t)\) 是平稳的——即不随时间变化。

算法 2.1 递归贝叶斯估计 1. 函数 RecursiveBayesianEstimation 2. 对所有 s,\(b_0(s) \leftarrow P(o_0 | s) P(s_0)\) 3. 归一化 \(b_0\) 4. 对 \(t \leftarrow 1\)\(\infty\) 5. 对所有 s,

\[ b_t(s) \leftarrow P(o_t | s) \sum_{s'} P(s | s') b_{t-1}(s') \]
  1. 归一化 \(b_t\)

若观测是连续而非离散,则 \(P(o | s)\) 将是概率密度而非概率质量。若状态是连续而非离散,则状态转移分布与 \(b\) 成为密度函数,第 5 行的求和变为积分——一般而言这是难以精确求解的。这一困难正是 2.2.5 节要讨论的近似推理方法的动机之一。

线性动态系统(linear dynamical system)——其状态转移与观测分布都是线性高斯——中做滤波,事实上可以做到精确。若 \(b_{t-1}\) 被表示为正态分布,则可以证明第 5 行的积分会得到高斯的后验 \(b_t\)卡尔曼滤波(Kalman filter)是用于线性动态系统的著名滤波算法,它能以更新 \(b_t\) 的均值与协方差的方式实现这一过程。在卡尔曼滤波下,每一步的均值与协方差更新都有一组闭式公式——这使得它在工程上既高效又精确,是目标跟踪、信号处理、控制系统中的核心算法。卡尔曼滤波的成功也提示了"在特定结构(线性高斯)下,推理可以保持精确"这一思想;脱离该结构则必须借助 2.2.5 节介绍的近似方法。

2.2.3 Exact Inference

让我们重访图 2.9 的网络并尝试精确地推断 \(P(b^1 | d^1, c^1)\)。由条件概率的定义,我们知道

\[ P(b^1 | d^1, c^1) = \frac{P(b^1, d^1, c^1)}{P(d^1, c^1)} \tag{2.34} \]

下面的讨论将聚焦于分子,因为计算分子的过程同样可以套用于分母——分母只是分子的一个特例(把 B 也边缘化掉)。由全概率公式,

\[ P(b^1, d^1, c^1) = \sum_{s} \sum_{e} P(b^1, s, e, d^1, c^1) \tag{2.35} \]

把隐变量求和消去的这一过程称为边缘化(marginalization)。由链式法则,

\[ P(b^1, d^1, c^1) = \sum_{s} \sum_{e} P(b^1) P(s) P(e | b^1, s) P(d^1 | e) P(c^1 | e) \tag{2.36} \]

精确推理的难点在于:必须把隐变量求和消去。在式 (2.36) 中我们只需消去两个变量,但在更大的网络中这一求和可能变得不可行。所需相加的项数会随隐变量数指数增长——尽管对许多贝叶斯网络而言,我们可以利用模型的结构来提升推理效率。隐变量数较多时,"暴力求和"在计算上不可行,但变量消去、连接树、置信传播等算法都利用了"网络结构"这一额外信息来削减计算量。

作为一个"特殊网络表示如何带来高效推理"的极端例子:假设我们有一个由二元变量 \(X_1, \dots, X_n\) 组成的贝叶斯网络,且网络中没有任何箭头。我们希望计算

\[ P(x_1^0) = \sum_{x_2} \cdots \sum_{x_n} P(x_1^0) P(x_2) \cdots P(x_n) \tag{2.37} \]

这里共有 \(2^{n-1}\) 项,每一项都是 n 个因子的乘积。当然,现实中我们没有理由先去构造联合分布再求和消去隐变量——\(P(x_1^0)\) 可以直接从指定 \(P(X_1)\) 的表里读出。

在更复杂的贝叶斯网络中,有多种方法可以执行高效推理。一种方法称为变量消除(variable elimination),它按序消去隐变量。下面我们用图 2.9 中贝叶斯网络的 \(P(B | d^1, c^1)\) 计算来说明变量消去算法。与网络中各节点关联的条件概率分布可以用下列各表表示:

\(T_1(B), T_2(S), T_3(E, B, S), T_4(D, E), T_5(C, E) \tag{2.38}\)

因为 D 与 C 是观测变量,最后两张表被替换为 \(T_6(E)\)\(T_7(E)\),仅保留 \(D = 1\)\(C = 1\) 的行。接着我们按序消去隐变量;不同策略可用于选择消去顺序,但本例中我们采用"先 E 后 S"的顺序。要消去 E,我们先收集所有涉及 E 的表:

\(T_3(E, B, S), T_6(E), T_7(E) \tag{2.39}\)

接着从这些表的乘积中把 E 求和消去,得到一张新表:

\[ T_8(B, S) = \sum_{e} T_3(e, B, S) T_6(e) T_7(e) \tag{2.40} \]

此时可以丢弃 \(T_3\)\(T_6\)\(T_7\),因为它们所含的全部信息都已经包含在 \(T_8\) 中。接下来消去 S:再次收集所有涉及 S 的表,从这些表的乘积中把 S 求和消去:

\[ T_9(B) = \sum_{s} T_2(s) T_8(B, s) \tag{2.41} \]

丢弃 \(T_2\)\(T_8\);此时剩下 \(T_1(B)\)\(T_9(B)\)。我们只需对这两张表的乘积做归一化,即可得到 \(P(B | d^1, c^1)\)

算法 2.2 给出了针对贝叶斯网络 \(\mathcal{B}\)、查询变量集 \(\mathcal{Q}\) 与观测值 \(\mathbf{o}\) 的变量消去算法梗概。对许多网络来说,变量消去使推理所需的时间与网络规模呈线性关系;但在最坏情况下它具有指数时间复杂度。影响计算量的关键是变量消去的顺序。选择最优的消去顺序——事实证明——是 NP-hard 的,也就是说在最坏情况下不能在多项式时间内完成(见 2.2.4 节)。即便我们找到了最优的消去顺序,变量消去仍可能需要指数次计算。变量消去的启发式一般试图最小化第 6 行生成的中间表中所涉及的变量数。

算法 2.2 贝叶斯网络中的变量消去 1. 函数 VariableElimination(\(\mathcal{B}, \mathcal{Q}, \mathbf{o}\)) 2. \(T \leftarrow\)\(\mathcal{B}\) 中各节点关联的条件概率表集合 3. 从 \(T\) 中所有表里移除与 \(\mathbf{o}\) 不一致的行 4. 对 \(i \leftarrow 1\)\(n\) 5. \(T' \leftarrow T\) 中所有涉及 \(X_i\) 的表 6. \(T \leftarrow\)\(X_i\)\(T'\) 中各表的乘积中消去 7. 从 \(T\) 中移除 \(T'\) 并把结果加入 \(T\) 8. \(T \leftarrow T\) 中剩余表的乘积 9. \(P(\mathcal{Q} | \mathbf{o}) \leftarrow\) 归一化 \(T\) 10. 返回 \(P(\mathcal{Q} | \mathbf{o})\)

一种被称为置信传播(belief propagation)的推理方法通过网络传递"消息"进行工作。置信传播需要线性时间,但仅在网络不含无向环的情况下给出精确答案。如果网络含有无向环,那么可以通过把多个变量合并为单个节点(使用所谓的连接树算法(junction tree algorithm))把它转化为树形结构。如果需要合并到任一节点中的变量数很小,那么推理就可以被高效地完成。

2.2.4 Complexity of Exact Inference

某些问题的求解难度可按复杂度类归类。本书后续将频繁出现的重要类包括:

  • P:能在多项式时间内求解的问题。
  • NP:其解能在多项式时间内验证的问题。
  • NP-hard:至少与 NP 中最难问题同等难度的问题。
  • NP-complete:既是 NP-hard 又属于 NP 的问题。

这些复杂度类的形式化定义相当复杂——它涉及"图灵机在多项式时间内能做什么"以及"问题之间的归约关系"等概念。一般认为 \(P \neq NP\),但这一结论尚未被证明,也是数学中最重要的未解问题之一。事实上,现代密码学依赖于"目前没有已知的多项式时间算法来求解 NP-hard 问题"这一事实——RSA、椭圆曲线等公钥密码系统的安全性都建立在这一假设之上。图 2.13 描绘了在 \(P \neq NP\) 假设下各复杂度类之间的关系:P 是 NP 的子集,NP-hard 与 NP 的交集即 NP-complete。

其中一个 NP-complete 问题——3SAT——来证明贝叶斯网络中的推理是 NP-hard。3SAT 是第一个已知的 NP-complete 问题:它要判断一个布尔公式是否可满足。该布尔公式由合取(\(\wedge\))、析取(\(\vee\))和取反(\(\neg\))组成,涉及 n 个布尔变量 \(x_1, \dots, x_n\)字面量(literal)是变量 \(x_i\) 本身或其取反 \(\neg x_i\)3SAT 子句是至多三个字面量的析取;例如 \(x_3 \vee \neg x_5 \vee x_6\)。一个 3SAT 公式是若干 3SAT 子句的合取,形如

\[ F(x_1, x_2, x_3, x_4) = (x_1 \vee x_2 \vee x_3) \wedge (\neg x_1 \vee \neg x_2 \vee x_3) \wedge (x_2 \vee \neg x_3 \vee x_4) \tag{2.42} \]

3SAT 的挑战在于:是否存在对变量真值的一种赋值使得公式为真。在式 (2.42) 中,

\(F(\text{true, false, false, true}) = \text{true} \tag{2.43}\)

因此该公式可满足。尽管对某些 3SAT 问题——有时仅靠快速观察——就能轻易找到可满足赋值,但 3SAT 一般而言是难以求解的。一种判断可满足赋值是否存在的方式是枚举所有 \(2^n\) 种变量真值组合。尽管"判断是否存在可满足真值赋值"是困难的,但"验证一个真值赋值是否满足公式"可以在线性时间内完成——这一"易验证 / 难求解"的非对称性正是 NP-complete 问题的本质特征,也是 Cook–Levin 定理 1971 年所证明的核心结论:所有 NP 问题都可以在多项式时间内归约到 3SAT。

从一个任意的 3SAT 问题构造贝叶斯网络是很容易的。图 2.14 是式 (2.42) 的一个贝叶斯网络表示。变量以 \(X_{1:4}\) 表示,子句以 \(C_{1:3}\) 表示。变量上的分布是均匀的。表示子句的节点以相关变量作为其父节点。由于这是 3SAT 问题,每个子句节点恰好有三个父节点。每个子句节点对"不满足该子句"的赋值赋予概率 0,对"满足该子句"的所有赋值赋予概率 1。剩下的节点在其所有父节点为真时对真赋予概率 1。原 3SAT 问题可满足当且仅当 \(P(y^1) > 0\)。因此贝叶斯网络中的推理至少与 3SAT 一样难——也就是说,任意一个 3SAT 实例都可以被"翻译"为一个贝叶斯网络上的推理问题;反之,一个高效的贝叶斯推理算法将立即给出高效的 3SAT 求解器。

从 2.2.4 节开始,作者把"贝叶斯网络中的精确推理是 NP-hard"这一事实给出显式证明——通过将 3SAT 归约为贝叶斯推理。这一"负面"结论在工程上不阻碍实际部署——许多网络类型仍可高效推理——但提醒我们"在最坏情况下没有通用高效算法"。因此 2.2.5 节转向近似推理:通过放弃一定的精度换取计算可行,这一策略在工业级贝叶斯网络系统中被广泛采用。

之所以要费力证明"贝叶斯网络中的推理是 NP-hard",是为了让我们清楚地知道:在所有贝叶斯网络上寻找一种高效、精确的推理算法是徒劳的。因此,过去几十年中,研究工作的重点转向了下一节要讨论的近似推理方法——这些方法牺牲一定的精度以换取计算可行。这一"理论上的负面结论"并不阻碍实际应用——对许多工程上关心的网络类型(如多叉树、低树宽网络、稀疏网络),精确推理仍是可行的;只是在最坏情况下我们需要退而求其次。

2.2.5 Approximate Inference

最简单的近似推理方法之一涉及从贝叶斯网络所表示的联合分布中采样(sampling)。第一步是为网络中的节点找到一个拓扑序(topological sort)。有向无环图中节点的拓扑序是一个有序列表,使得若存在边 \(A \to B\),则 A 在表中排在 B 之前。例如,对图 2.9 的网络,一种拓扑序是 B, S, E, D, C。拓扑序一定存在,但不一定唯一。另一种对图 2.9 的拓扑序是 S, B, E, C, D。算法 2.3 给出了对图 G 寻找拓扑序的算法。

算法 2.3 拓扑排序 1. 函数 TopologicalSort(\(G\)) 2. \(n \leftarrow\) \(G\) 中节点数 3. \(L \leftarrow\) 空列表 4. 对 \(i \leftarrow 1\)\(n\) 5. \(X \leftarrow\) 任何不在 \(L\) 中、且其所有父节点都已在 \(L\) 中的节点 6. 将 \(X\) 追加到 \(L\) 末尾 7. 返回 \(L\)

一旦我们有了拓扑序,就可以从条件概率分布中开始采样。假设拓扑序给出 \(X_{1:n}\)。算法 2.4 给出了从贝叶斯网络 \(\mathcal{B}\) 中采样的方法。在第 4 行,我们从 \(X_i\) 在给定"已赋值的父节点取值"的条件分布中抽取一个样本。由于 \(X_{1:n}\) 是拓扑序,我们知道 \(X_i\) 的所有父节点都已被实例化,因此可以完成这一采样。

算法 2.4 贝叶斯网络上的直接采样 1. 函数 DirectSample(\(\mathcal{B}\)) 2. \(X_{1:n} \leftarrow\) \(\mathcal{B}\) 中节点的拓扑序 3. 对 \(i \leftarrow 1\)\(n\) 4. \(x_i \leftarrow\)\(P(X_i | pa_{x_i})\) 中随机采样 5. 返回 \(x_{1:n}\)

表 2.4 展示了从图 2.9 的网络中采到的 10 个随机样本。我们感兴趣的是推断 \(P(b^1 | d^1, c^1)\)。10 个样本中只有 2 个(表中以箭头标出)与观测 \(d^1\)\(c^1\) 一致;其中一份样本的 B = 1,另一份的 B = 0。从这两份样本我们推断 \(P(b^1 | d^1, c^1) = 0.5\)。当然,实践中我们需要的样本量远不止 2 份才能准确估计 \(P(b^1 | d^1, c^1)\)

直接采样存在的问题是:我们可能浪费大量时间去生成与观测不一致的样本——尤其当观测本身不太可能发生时。另一种方法称为似然加权(likelihood weighting):它生成与观测一致的带权样本。我们先找到一个拓扑序,然后依次从条件分布中采样。似然加权与直接采样的唯一区别在于如何处理已观测变量:我们不通过采样决定其取值,而是将变量直接固定为观测值,并相应地调整样本的权重。样本的权重就是已观测节点上条件概率的乘积。算法 2.5 总结了对贝叶斯网络 \(\mathcal{B}\) 与观测 \(o_{1:n}\) 的过程。若 \(o_i\) 未被观测,则 \(o_i \leftarrow \text{nil}\)

算法 2.5 贝叶斯网络上的似然加权采样 1. 函数 LikelihoodWeightedSample(\(\mathcal{B}, o_{1:n}\)) 2. \(X_{1:n} \leftarrow\) \(\mathcal{B}\) 中节点的拓扑序 3. \(w \leftarrow 1\) 4. 对 \(i \leftarrow 1\)\(n\) 5. 若 \(o_i = \text{nil}\) 6. \(x_i \leftarrow\)\(P(X_i | pa_{x_i})\) 中随机采样 7. 否则 8. \(x_i \leftarrow o_i\) 9. \(w \leftarrow w \times P(x_i | pa_{x_i})\) 10. 返回 \((x_{1:n}, w)\)

表 2.5 展示了从图 2.9 网络采出的 5 份似然加权样本。我们从 \(P(B)\)\(P(S)\)\(P(E | B, S)\) 中采样——与直接采样相同。当遇到 D 与 C 时,我们令 D = 1、C = 1。若样本中 E = 1,则权重为 \(P(d^1 | e^1) P(c^1 | e^1)\);否则权重为 \(P(d^1 | e^0) P(c^1 | e^0)\)。若假设

\(P(d^1 | e^1) P(c^1 | e^1) = 0.95 \tag{2.44}\)

\(P(d^1 | e^0) P(c^1 | e^0) = 0.01 \tag{2.45}\)

则从表 2.5 中的样本可得

\[ P(b^1 | d^1, c^1) \approx \frac{0.95}{0.95 + 0.95 + 0.01 + 0.01 + 0.95} = 0.331 \tag{2.46,2.47} \]

虽然似然加权让所有样本都与观测一致,但它仍可能是浪费的。考虑图 2.15 所示的简单化学探测贝叶斯网络,假设我们探测到目标化学物质。我们想推断 \(P(c^1 | d^1)\)。因为网络很小,可以直接用贝叶斯公式精确计算该概率:

\[ P(c^1 | d^1) = \frac{P(d^1 | c^1) P(c^1)}{P(d^1 | c^1) P(c^1) + P(d^1 | c^0) P(c^0)} = \frac{0.999 \times 0.001}{0.999 \times 0.001 + 0.001 \times 0.999} = 0.5 \tag{2.48-2.50} \]

若使用似然加权,那么 99.9% 的样本将具有 C = 0、权重为 0.001;在我们抽到一个 C = 1 的样本(其权重为 0.999)之前,对 \(P(c^1 | d^1)\) 的估计将一直是 0。

另一种方法是吉布斯采样(Gibbs sampling),这是一种马尔可夫链蒙特卡洛(Markov chain Monte Carlo)技术。与前面讨论的采样方法不同,这种方法生成的样本并不独立:下一份样本在概率意义上依赖于当前样本,因此样本序列构成一个马尔可夫链。可以证明在极限下样本恰从给定观测下的未观测变量联合分布中抽取。

初始样本可以随机生成,已观测变量被设为观测值。算法 2.6 给出了从现有样本 \(x_{1:n}\) 出发、给定贝叶斯网络 \(\mathcal{B}\) 与观测 \(o_{1:n}\) 生成新样本 \(x'_{1:n}\) 的方法。与直接采样不同,我们可以使用任意节点顺序——不必是拓扑序。在这一顺序下,每次基于"其他变量的当前值"对一个变量进行更新。为生成 \(x'_i\) 的值,我们从 \(P(X_i | x'_{-i})\) 中采样,其中 \(x'_{-i}\) 表示除 \(X_i\) 之外所有其他变量的值。要在贝叶斯网络 \(\mathcal{B}\) 中计算 \(P(X_i | x'_{-i})\),可以使用算法 2.7。由于只需考虑 \(X_i\) 的 Markov 毯(见 2.1.5 节),该计算可以高效完成。

算法 2.6 贝叶斯网络上的吉布斯采样 1. 函数 GibbsSample(\(\mathcal{B}, o_{1:n}, x_{1:n}\)) 2. \(X_{1:n} \leftarrow\) \(\mathcal{B}\) 中节点的一个顺序 3. \(x'_{1:n} \leftarrow x_{1:n}\) 4. 对 \(i \leftarrow 1\)\(n\) 5. 若 \(o_i = \text{nil}\) 6. \(x'_i \leftarrow\)\(P(X_i | x'_{-i})\) 中随机采样 7. 否则 8. \(x'_i \leftarrow o_i\) 9. 返回 \(x'_{1:n}\)

算法 2.7 在其他节点均已观测的条件下某节点的分布 1. 函数 DistributionAtNode(\(\mathcal{B}, X_i, x_{-i}\)) 2. \(T \leftarrow\) \(\mathcal{B}\) 中所有涉及 \(X_i\) 的条件概率表 3. 从 \(T\) 的所有表中移除与 \(x_{-i}\) 不一致的行 4. \(T \leftarrow\) \(T\) 中剩余表的乘积 5. \(P(X_i | x_{-i}) \leftarrow\) 归一化 \(T\) 6. 返回 \(P(X_i | x_{-i})\)

图 2.16 对直接采样、似然加权与吉布斯采样下 \(P(c^1 | d^1)\) 估计的收敛过程做了对比。直接采样收敛最慢:当样本与观测不一致时,其估计长期不发生变化。似然加权采样在该例中收敛较快;当采到 C = 1 的样本时估计会出现尖峰,随后逐渐下降。吉布斯采样在该例中很快收敛到了真值 0.5。

如前所述,吉布斯采样与其他 MCMC 方法一样,是在极限下从目标分布中产生样本。实践中我们必须让吉布斯链先运行一段时间——称为老化期(burn-in period)——再让其收敛到稳态分布。老化期内产生的样本通常会被丢弃。此外,由于样本之间可能存在相关性,常用的做法是稀疏化(thinning)——即仅保留每隔 k 份的样本。

其他近似推理方法不涉及生成样本。例如,一种被称为环路置信传播(loopy belief propagation)的置信传播形式,可以用于具有无向环的网络以做近似推理。尽管不保证精确,环路置信传播在实践中效果通常很好,正在成为贝叶斯网络近似推理最流行的方法之一。

2.3 Parameter Learning

到目前为止,本章都假设概率模型的参数与结构是已知的。本节讨论如何从数据中学习模型的参数;2.4 节进一步讨论如何从数据中学习模型的结构。这两节构成本章的"学习"部分,与 2.1–2.2 节的"表示 + 推理"部分形成对照。具体而言,2.3.1 节讨论最大似然参数学习——一种把"模型对数据的拟合程度"作为目标函数的方法;2.3.2 节讨论贝叶斯参数学习——在最大似然框架之上引入参数先验;2.3.3 节讨论非参数学习——不预设固定参数个数、参数数量随数据增长的核密度估计。

学习问题的关键挑战是"在数据有限时如何避免过拟合"——最大似然估计在数据极少时容易给出极端值(0 或 1),而贝叶斯方法通过引入先验缓解了这一问题。另一个挑战是"模型复杂度"——过于复杂的模型会拟合噪声,过于简单的模型又拟合不足;这一权衡在 2.4 节的结构学习中通过"贝叶斯评分"得到系统处理。本节所讨论的"参数学习"在结构上假定模型形式是固定的——也就是说,我们已经知道变量之间的因果图,只是不确定每个节点的条件概率参数;而 2.4 节将进一步讨论"连图本身都未知"的情形。参数学习在工程上之所以重要,是因为现实问题中我们往往能获得大量数据样本,但参数的具体数值无法从物理原理直接推出——必须借助数据驱动的方法。例如在空管领域,飞行员的反应延迟分布需要从历史飞行数据中估计;这种估计的质量直接决定决策支持系统的可靠性。在机器学习领域,参数学习也是"模型训练"这一步骤的数学抽象——给定一个模型族,参数学习回答"如何在数据上挑出最佳参数"。本节余下三个小节将分别讨论最大似然、贝叶斯、非参数三种参数学习策略——前两种属于"参数方法",第三种属于"非参数方法";它们的根本差异在于"模型预设的灵活性"。从哲学的角度看,这三种方法对应三种不同的认识论立场:MLE 相信"数据会说话"、贝叶斯承认"先验也说话"、非参数则"让模型复杂度由数据决定"。这三类方法的工程选择往往由"数据量 + 领域知识 + 实时性要求"三者共同决定——这是一个多目标权衡问题。本节关注的"参数"是模型中具体控制分布形状的数值——它们与 2.4 节关注的"结构"是相互独立但又紧密耦合的两个层面:固定结构下学习参数,固定参数下学习结构,循环迭代这两者直至收敛——这是工程上常用的"结构 + 参数联合学习"的两阶段迭代策略。这种迭代策略在工程实现上对应 EM 算法的"广义"版本——后者把"隐变量"也视作一类可学习的"参数",从而把"结构 + 参数 + 隐变量"三者的联合学习统一在同一个框架内。这种"结构 + 参数 + 隐变量"三者联合学习的统一性是概率图模型相对其他机器学习框架的一个重要优势——它使得模型的不同组成部分可以用同一套算法(如 EM、变分推断)来估计。值得注意的是,参数学习与 2.2 节的推理也存在"对偶"关系——给定参数求后验是推理,给定数据求参数是学习。这两类操作在数学形式上紧密对应:推理是"积分掉隐变量",学习是"积分掉参数";两者的算法实现也常常可以共享。

2.3.1 Maximum Likelihood Parameter Learning

假设随机变量 C 表示某航班是否会导致空中相撞,我们希望估计分布 \(P(C)\)。因为 C 只取 0 或 1,所以只需估计参数 \(\theta = P(c^1)\)。我们想从数据 \(\mathcal{D}\) 中推断 \(\theta\)。假设我们有一个跨十年的历史数据库,并已知共有 n 次航班、m 次空中相撞。直觉上,\(\theta\) 基于数据 \(\mathcal{D}\) 的良好估计是 \(m/n\)。这一估计对应着最大似然估计(maximum likelihood estimate, MLE)

\(\hat{\theta} = \arg\max_{\theta} P(\mathcal{D} | \theta) \tag{2.51}\)

其思想是:选择使"在当前参数下观测到数据"这一事件的概率最大的 \(\theta\)。把"先假设 \(\theta\),再问数据出现的可能性"作为目标函数——"让数据看起来最像"的那个 \(\theta\) 就是我们的估计。

n 次航班中发生 m 次相撞的概率由二项分布给出:

\[ P(\mathcal{D} | \theta) = \frac{n!}{m! (n - m)!} \theta^m (1 - \theta)^{n-m} \propto \theta^m (1 - \theta)^{n-m} \tag{2.52,2.53} \]

其中 \(\propto\) 表示"忽略与 \(\theta\) 无关的常数因子"——这一简化是可行的,因为我们只关心使 \(P(\mathcal{D} | \theta)\) 最大的 \(\theta\),而与 \(\theta\) 无关的常数不会改变最优点。最大似然估计 \(\hat{\theta}\) 是使式 (2.53) 取最大值的 \(\theta\) 值。最大化式 (2.53) 等价于最大化似然函数的对数——通常称为对数似然(log-likelihood),常记作 \(\ell(\theta)\)。取对数的合理性来自对数函数在 (0, 1] 上的单调性:最大化原函数与其对数函数的最优解相同。

\[ \ell(\theta) \propto \ln \theta^m (1 - \theta)^{n-m} = m \ln \theta + (n - m) \ln(1 - \theta) \tag{2.54,2.55} \]

我们可以用求函数最大值的标准技巧:将 \(\ell\)\(\theta\) 的一阶导数设为 0 并解出 \(\theta\)。导数为

\[ \frac{\partial \ell(\theta)}{\partial \theta} = \frac{m}{\theta} - \frac{n - m}{1 - \theta} \tag{2.56} \]

设导数为 0 并解出 \(\hat{\theta}\)

\[ \frac{m}{\hat{\theta}} - \frac{n - m}{1 - \hat{\theta}} = 0 \tag{2.57} \]

经过几步代数即可验证 \(\hat{\theta} = m/n\)

对取 k 个值的变量 X 求最大似然估计也是直接的。若 \(m_{1:k}\) 是 k 个不同取值的观测计数,则 \(P(x_i | m_{1:k})\) 的最大似然估计为

\[ \hat{\theta}_i = \frac{m_i}{\sum_{j=1}^{k} m_j} \tag{2.58} \]

最大似然估计也可以应用于连续分布。假设我们有空速测量 \(v_{1:n}\),它们是图 2.12 中用于生成类条件分布的 n 条航迹的空速。尽管密度显然并不精确服从高斯分布,我们还是尝试以最大似然估计把高斯模型拟合到这些数据。均值 \(\mu\) 与方差 \(\sigma^2\) 的对数似然为

\[ \ell(\mu, \sigma^2) \propto -n \ln \sigma - \sum_i \frac{(v_i - \mu)^2}{2\sigma^2} \tag{2.59} \]

同样,我们对参数求导并令其为 0,解出最大似然估计:

\[ \frac{\partial \ell(\mu, \sigma^2)}{\partial \mu} = \sum_i \frac{(v_i - \hat{\mu})}{\hat{\sigma}^2} = 0 \tag{2.60} \]
\[ \frac{\partial \ell(\mu, \sigma^2)}{\partial \sigma^2} = -\frac{n}{\hat{\sigma}} + \sum_i \frac{(v_i - \hat{\mu})^2}{\hat{\sigma}^3} = 0 \tag{2.61} \]

经过一些代数整理,可得

\[ \hat{\mu} = \frac{\sum_i v_i}{n} \tag{2.62} \]
\[ \hat{\sigma}^2 = \frac{\sum_i (v_i - \hat{\mu})^2}{n} \tag{2.63} \]

图 2.17 展示了一个以最大似然估计 \(\hat{\mu} = 100.2\) 节、\(\hat{\sigma} = 31\) 节为参数的高斯分布。图中还把图 2.12 中的"真实"分布一同画出以供对比;该例中高斯是对真实分布一个相当合理的近似——尽管真实分布在两侧尾部呈现多峰性,但高斯模型在中心区域对真实密度的拟合程度已经足够用于工程决策。最大似然估计在数据量充足时具有"渐近无偏、渐近有效"等优良统计性质,但在数据稀少时容易给出"过自信"的极端值——这一不足正是 2.3.2 节要讨论的贝叶斯参数学习的根本动因。最大似然估计的另一个隐含限制是"模型必须正确指定"——如果真实分布不属于我们假设的分布族,那么无论数据再多,MLE 都只会收敛到分布族内"最接近"真实分布的那一点,而不是真实分布本身。这种"模型误设"的局限性在工程上提示我们:当有理由怀疑真实分布偏离预设形式时,应当考虑 2.3.3 节介绍的非参数方法,或者 2.3.2 节介绍的贝叶斯方法(后者通过先验引入额外的灵活性)。从优化角度看,最大似然估计对应"最小化 KL 散度 \(D_{KL}(P_{\text{data}} \| P_\theta)\) 的等价形式"——这一观察把 MLE 与"分布匹配"的目标联系起来:MLE 的目标就是"让模型分布尽可能接近真实数据分布",但只优化一个方向(\(P_\theta\) 向着 \(P_{\text{data}}\) 靠拢)。

2.3.2 Bayesian Parameter Learning

尽管最大似然估计对许多应用可能已经足够,但在数据量有限时它存在一些严重缺陷。例如,假设我们的航空安全数据库仅限于过去一周的事件,并发现没有空中相撞记录。若 \(\theta\) 是"航班导致空中相撞"的概率,则最大似然估计会是 \(\hat{\theta} = 0\)。认为"航班不存在相撞的可能性"并非合理结论,除非我们的先验假设——例如——要么所有航班都安全、要么所有航班都相撞。

参数学习的贝叶斯方法涉及对 \(\theta\) 的后验进行估计,并可视为在贝叶斯网络中所做的推理。例如,图 2.18 表示碰撞概率估计问题,其中观测变量 \(O_i\) 在第 i 次航班导致相撞时取 1、否则取 0。我们假设观测变量之间相互条件独立。我们必须指定 \(p(\theta)\)\(P(O_i | \theta)\)。若想用均匀先验,可以令密度 \(p(\theta) = 1\)。令 \(P(o_i^1 | \theta) = \theta\)

我们可以按贝叶斯网络中推理的标准方法进行。这里假设均匀先验:

\[ p(\theta | o_{1:n}) \propto p(\theta, o_{1:n}) = p(\theta) \prod_{i=1}^{n} P(o_i | \theta) = \prod_{i=1}^{n} P(o_i | \theta) = \prod_{i=1}^{n} \theta^{o_i} (1 - \theta)^{1 - o_i} = \theta^m (1 - \theta)^{n - m} \tag{2.64-2.68} \]

后验正比于 \(\theta^m (1 - \theta)^{n - m}\),其中 m 是数据中相撞的次数。为求归一化常数,对下式积分:

\[ \int_0^1 \theta^m (1 - \theta)^{n-m} d\theta = \frac{\Gamma(m + 1) \Gamma(n - m + 1)}{\Gamma(n + 2)} \tag{2.69} \]

其中 \(\Gamma\)伽马函数——是阶乘在实数上的一般化。若 n 为整数,则 \(\Gamma(n) = (n-1)!\)。把归一化考虑进来,有

\[ p(\theta | o_{1:n}) = \frac{\Gamma(n + 2)}{\Gamma(m + 1) \Gamma(n - m + 1)} \theta^m (1 - \theta)^{n - m} = \text{Beta}(\theta | m + 1, n - m + 1) \tag{2.70,2.71} \]

贝塔分布 \(\text{Beta}(\alpha, \beta)\) 由参数 \(\alpha\)\(\beta\) 定义;该分布的曲线如图 2.19 所示。\(\text{Beta}(1, 1)\) 对应于取值在 0 到 1 之间的均匀分布。

方便的是,若用贝塔分布作为二项分布参数的先验,则后验仍为贝塔分布。具体而言,若先验为 \(\text{Beta}(\alpha, \beta)\),且我们做一次观测 \(o_i\),则 \(o_i = 1\) 时得到后验 \(\text{Beta}(\alpha + 1, \beta)\)\(o_i = 0\) 时得到 \(\text{Beta}(\alpha, \beta + 1)\)。因此若我们以 \(\text{Beta}(\alpha, \beta)\) 为先验开始,且数据中有 n 次航班、m 次相撞,则后验为 \(\text{Beta}(\alpha + m, \beta + n - m)\)。先验中的 \(\alpha\)\(\beta\) 有时被称为伪计数(pseudocounts)——因为它们在形式上与后验中两类结果的观测计数作用相似,尽管伪计数不必取整。

原则上,先验的选择应当在不知道用于计算后验的数据的前提下进行。均匀先验在实践中往往效果很好,但若能获得专家知识,则可将其编码进先验。例如,假设我们有一枚略弯的硬币,想估计 \(\theta\)(硬币正面朝上的概率)。在收集任何翻硬币数据之前,我们的初始信念是 \(\theta\) 大致在 0.5 附近。比起用均匀先验 \(\text{Beta}(1, 1)\),我们或许会用 \(\text{Beta}(2, 2)\)(图 2.19),它给 0.5 附近的值赋予更多权重。若对"接近 0.5"这一估计更有信心,我们可以通过增大伪计数来减小先验的方差。\(\text{Beta}(10, 10)\)\(\text{Beta}(2, 2)\) 集中得多。然而总体上,先验的重要性会随着用于计算后验的数据量增加而下降。若观测了 n 次翻硬币且其中 m 次为正面,那么 \(\text{Beta}(1 + m, 1 + n - m)\)\(\text{Beta}(10 + m, 10 + n - m)\) 之间的差异在翻硬币数千次之后可以忽略。

狄利克雷分布(Dirichlet distribution)是贝塔分布的推广,可用于估计离散分布的参数。假设 X 是一个取 1 到 n 整数值的离散随机变量。我们把分布的参数定义为 \(\theta_{1:n}\),其中 \(P(x_i) = \theta_i\)。参数之和必须为 1,因此仅有前 n-1 个参数是独立的。狄利克雷分布可用于表示先验与后验分布,由 \(\alpha_{1:n}\) 参数化。其密度为

\[ \text{Dir}(\theta_{1:n} | \alpha_{1:n}) = \frac{\Gamma(\alpha_0)}{\prod_{i=1}^{n} \Gamma(\alpha_i)} \prod_{i=1}^{n} \theta_i^{\alpha_i - 1} \tag{2.72} \]

其中 \(\alpha_0\) 表示参数 \(\alpha_{1:n}\) 的求和。当 n = 2 时,显然式 (2.72) 等价于贝塔分布。

一种常见做法是采用均匀先验——即把所有的狄利克雷参数 \(\alpha_{1:n}\) 设为 1。对称狄利克雷分布(symmetric Dirichlet distribution)是指所有参数相同的分布。和贝塔分布一样,狄利克雷分布中的参数常被称为伪计数。

\(\theta_{1:n}\) 上的先验为 \(\text{Dir}(\alpha_{1:n})\),且 \(X = i\) 的观测有 \(m_i\) 次,则后验为

\[ p(\theta_{1:n} | \alpha_{1:n}, m_{1:n}) = \text{Dir}(\theta_{1:n} | \alpha_1 + m_1, \dots, \alpha_n + m_n) \tag{2.73} \]

如我们所见,对二元与离散随机变量来说,贝叶斯参数估计很直接——它只涉及对数据中各类结果的计数。对任何参数化分布,都可以用贝叶斯公式推断参数上的分布。视先验选择与参数化分布的形式而定,对参数空间上的后验进行计算也可能是可解析完成的。

2.3.3 Nonparametric Learning

前两节都假设概率模型的形式是固定的,所需学习的也只是一个固定的参数集合。一种替代方案基于非参数方法(nonparametric methods)——其参数个数会随数据量而变化。参数方法非参数方法的关键区别在于:前者假定模型的具体形式(如高斯、二项),并对该形式的参数进行估计;后者不预设固定形式,而是让模型的复杂度由数据本身决定。核密度估计(kernel density estimation, KDE)是最常见的非参数方法之一。

给定观测 \(o_{1:n}\),核密度估计把密度表示为

\[ p(x) = \frac{1}{n} \sum_{i=1}^{n} K(x - o_i) \tag{2.74} \]

其中 \(K\)核函数(kernel function),它积分为 1。核函数用来为接近观测数据点的值赋予更大的密度——其直观含义是把每个观测点"涂抹"成以它为中心的小分布,再把所有"涂抹"加在一起形成总体密度。核函数一般是对称的,即 \(K(x) = K(-x)\)——这一对称性确保密度估计不会"偏向"数据一侧。常用的核是零均值高斯分布;当使用高斯作为核时,其标准差常被称为带宽(bandwidth)。高斯核密度是平滑的,与之前讨论的分段均匀分布(图 2.1)形成对照。较大的带宽通常会带来更平滑的密度——这在工程上意味着更大的"抹平"效果,但也可能掩盖数据中真实的多峰结构。贝叶斯方法可以基于数据来选取合适的带宽——这是"用先验 + 数据决定超参数"这一思路的又一次应用。核密度估计的另一个优势在于其"非参数"特性:随着样本量增加,估计的精度会持续提升,不会因为模型预设的形式过于简单而遭遇"模型误设"(model misspecification)问题。核密度估计的局限在于:当数据维度较高(>5)时,每个观测点周围所需的"覆盖体积"急剧增大,需要的样本量随之指数级增长——这就是所谓的"维数灾难(curse of dimensionality)"。这也是为什么实际工程问题中常借助参数模型(如高斯混合)来控制维度——以"轻度模型误设"换取对样本量的可控要求。从概率论的角度看,KDE 与高斯混合都属于"无限维"或"高维"参数化模型的特例——前者把每个观测点都作为一个高斯成分,后者则用一个固定数量的高斯成分近似。在实际工程中,自适应带宽(adaptive bandwidth)是对固定带宽 KDE 的一种重要改进——它让每个观测点处的带宽根据局部数据密度自适应调整,从而在数据稀疏区更平滑、在数据稠密区更精细。这种局部自适应机制在工程上常通过"k 近邻距离"或"k 近邻逆距离"作为带宽的局部估计来实现——这也是非参数密度估计研究的一条主线。核密度估计在样本量趋于无穷时收敛到真实密度,这一理论性质在 1950 年代由 Parzen 与 Rosenblatt 分别独立证明——这为 KDE 在统计学习中的合法性奠定了基础。值得注意的是,KDE 在收敛速度上呈现出"非参数率"(n^{-4/(4+d)} 对 d 维数据)——这一速度比参数模型的对 n^{-1/2} 慢得多,因此"非参数"这一性质既是 KDE 的优势(无需预设形式),也是其劣势(需要更多数据才能达到同等精度)。在工程实践中,核密度估计常被用作"基线方法"——任何新的非参数方法都应至少在基准数据集上与 KDE 比较,以确认其是否真正带来了性能提升。

2.4 Structure Learning

前几节都假设贝叶斯网络的结构是事先已知的。本节讨论如何从数据中学习结构。学习贝叶斯网络结构的一种最大似然方法涉及:寻找使 \(P(G | \mathcal{D})\) 取最大值的图结构 \(G\),其中 \(\mathcal{D}\) 表示可用的数据。我们先说明如何基于 \(P(G | \mathcal{D})\) 计算一个贝叶斯网络评分,然后说明如何在该网络空间上搜索得分最高的网络。和贝叶斯网络中的推理一样,可以证明对一般的图与输入数据,学习贝叶斯网络的结构是 NP-hard 的——也就是说,即使我们愿意接受某种近似评分,也难以在多项式时间内找到全局最优的网络。这一复杂度结果由 Chickering 等人在 1990 年代给出,与 2.2.4 节的推理 NP-hard 结果并列——结构学习在工程上通常借助 2.4.2 节介绍的启发式搜索来应对这一负面结果。

结构学习比参数学习更具挑战性。参数学习中模型结构是固定的,搜索空间是连续的;结构学习中模型本身是搜索对象,搜索空间是离散的(所有可能的有向无环图)且巨大(n 节点时超指数级增长)。本节从 2.4.1 节起讨论贝叶斯评分(一种"既考虑模型复杂度、又考虑数据拟合程度"的评分函数),2.4.2 节讨论在有向无环图空间上的搜索策略,2.4.3–2.4.4 节讨论"哪些网络结构在数据上不可区分"这一深刻观察。结构学习的最终产物是一个有向无环图,它既能在数据上拟合得足够好,又能在结构上保持简洁——这正是贝叶斯评分所追求的两个目标的平衡。结构学习在工程上对应"当领域知识不足以手绘网络时,借助数据自动发现变量间的依赖结构"——例如生物信息学中基因调控网络的推断、故障诊断中故障传播图的学习,都属于这一范畴。结构学习与因果发现(causal discovery)密切相关——后者更进一步地追问"X 是否是 Y 的因";但纯结构学习只关心"哪些变量之间存在条件独立关系",而把因果方向留作"马尔可夫等价类"层面的不确定。本节余下四个小节构成本章的"结构学习工具箱":2.4.1 节是评分函数、2.4.2 节是搜索算法、2.4.3 节是不可区分性、2.4.4 节是状态空间优化。这四者构成一个完整的"评分 + 搜索 + 状态空间"三件套——任何结构学习算法都可以视为这三者的具体组合。这一工具箱在 2.4 节开头就给出,是因为结构学习涉及多个"看上去独立但实质耦合"的子问题——理解它们的耦合方式是工程上设计良好结构学习系统的关键。从工程实现的角度看,结构学习算法通常以"评分缓存 + 增量更新"作为性能优化核心——当一个搜索操作只改动图中少数几条边时,不必重新计算整个图的评分,只需对"受影响"的部分做局部更新——这一优化思路与 2.2.3 节变量消去的"增量计算"思想一脉相承。结构学习在工程上常被视为"机器学习中的组合优化"——它结合了图论、组合优化、信息论三个领域的工具,是机器学习应用数学工具的一个典型范例。结构学习的工程实践往往遵循"先粗后精"的两阶段流程:先用粗粒度搜索在较大的搜索空间内找到一个"差不多"的图结构,再用细粒度搜索在该结构附近做局部优化——这一两阶段策略在工程上比单一粒度的搜索更稳定、更高效。

2.4.1 Bayesian Structure Scoring

在讨论如何计算 \(P(G | \mathcal{D})\) 之前,我们需要先引入一些记法。我们假设贝叶斯网络中的 n 个变量 \(X_{1:n}\) 都是离散的,尽管一般而言这并非必须。用 \(r_i\) 表示 \(X_i\) 的可能取值数,用 \(q_i\) 表示 \(X_i\) 的父节点的可能取值数。若 \(X_i\) 无父节点,则 \(q_i = 1\)\(X_i\) 父节点的第 j 个赋值记作 \(\pi_{ij}\)

贝叶斯网络中共 \(\sum_{i=1}^{n} r_i q_i\) 个参数。每个参数写作 \(\theta_{ijk}\),且决定

\(P(X_i = k | \pi_{ij}) = \theta_{ijk} \tag{2.75}\)

尽管共有 \(\sum_{i=1}^{n} r_i q_i\) 个参数,但只有 \(\sum_{i=1}^{n} (r_i - 1) q_i\) 个是独立的。我们用 \(\theta\) 表示所有参数的集合。

\(m_{ijk}\) 表示在数据集 \(\mathcal{D}\) 中"在 \(\pi_{ij}\) 条件下 \(X_i = k\)"出现的次数。似然为

\[ P(\mathcal{D} | \theta, G) = \prod_{i=1}^{n} \prod_{j=1}^{q_i} \prod_{k=1}^{r_i} \theta_{ijk}^{m_{ijk}} \tag{2.76} \]

贝叶斯网络参数 \(\theta\) 上的先验可以分解。若有 \(\theta_{ij} = (\theta_{ij1}, \dots, \theta_{ij r_i})\),则

\[ p(\theta | G) = \prod_{i=1}^{n} \prod_{j=1}^{q_i} p(\theta_{ij}) \tag{2.77} \]

在某些弱假设下,先验 \(p(\theta_{ij})\) 可被证明服从式 (2.72) 中的狄利克雷分布。给定父节点第 j 个赋值的条件下,\(X_i\) 上的分布为 \(\text{Dir}(\alpha_{ij1}, \dots, \alpha_{ij r_i})\)

我们借助贝叶斯公式与全概率公式计算 \(P(G | \mathcal{D})\)

\[ P(G | \mathcal{D}) \propto P(G) P(\mathcal{D} | G) = P(G) \int P(\mathcal{D} | \theta, G) p(\theta | G) d\theta \tag{2.78,2.79} \]

将式 (2.76) 与式 (2.77) 的乘积对 \(\theta\) 积分后,可得

\[ P(G | \mathcal{D}) = P(G) \prod_{i=1}^{n} \prod_{j=1}^{q_i} \frac{\Gamma(\alpha_{ij0})}{\Gamma(\alpha_{ij0} + m_{ij0})} \prod_{k=1}^{r_i} \frac{\Gamma(\alpha_{ijk} + m_{ijk})}{\Gamma(\alpha_{ijk})} \tag{2.80} \]

其中

\[ \alpha_{ij0} = \sum_{k=1}^{r_i} \alpha_{ijk} \tag{2.81} \]
\[ m_{ij0} = \sum_{k=1}^{r_i} m_{ijk} \tag{2.82} \]

寻找使式 (2.79) 最大的 \(G\) 与寻找使贝叶斯评分(Bayesian score) 最大的 \(G\) 是一致的:

\[ \ln P(G | \mathcal{D}) = \ln P(G) + \sum_{i=1}^{n} \sum_{j=1}^{q_i} \left[ \ln \frac{\Gamma(\alpha_{ij0})}{\Gamma(\alpha_{ij0} + m_{ij0})} + \sum_{k=1}^{r_i} \ln \frac{\Gamma(\alpha_{ijk} + m_{ijk})}{\Gamma(\alpha_{ijk})} \right] \tag{2.83} \]

贝叶斯评分在数值上更便于计算,因为"对小数取对数再相加"比"把它们乘起来"要容易。许多软件库都能以合理精度直接计算伽马函数的对数。另一个细节是:贝叶斯评分中"独立参数 + 数据"的对数项组合可以写成"自信息"的形式——\(\log P(D | G) = \sum_{i,j,k} m_{ijk} \log \theta_{ijk} + \text{constant}\)——这一形式在信息论视角下有清晰含义:评分本质上是"在最佳参数下数据相对于模型的对数似然"减去"参数个数 × 数据量"的复杂度罚项。这种罚项在形式上与统计学习中的 AIC、BIC 等模型选择准则同构——贝叶斯评分可以视为一种基于先验选择的正则化形式。具体来说,AIC ≈ -2 log P + 2k(其中 k 是参数个数),BIC ≈ -2 log P + k log n(其中 n 是样本数);贝叶斯评分的复杂度罚项则更精细——它对每个参数都基于其先验做调整,因此对模型的不同部分施加不同的罚。

文献中已经探索了多种不同的图先验,但实践中常用的是均匀先验——即把所有的图先验视为相等——此时 \(\ln P(G)\) 可从式 (2.83) 的贝叶斯评分计算中略去。贝叶斯评分的一个有用性质是:即使在均匀图先验下,它也能在"模型复杂度"与"可用数据量"之间达成最优平衡。这种"自动权衡"特性使得贝叶斯评分成为结构学习中最常用的评分函数之一。

为说明贝叶斯评分如何平衡模型复杂度,考虑图 2.20 中的简单贝叶斯网络。A 的取值对 B 的取值有弱影响,C 与其他变量相互独立。我们从该"真实"模型中采样以生成数据 \(\mathcal{D}\),然后尝试学习模型结构。三个变量共有 \(2^3 = 8\) 种网络结构(实际上若考虑边是否存在则为 \(2^{3 \times 2} / \text{无环约束}\)),但我们将聚焦于以下几种模型的评分:

  • 真实模型:\(1 + 2 + 1 = 4\) 个独立参数;
  • 全连接模型 \(A \to B, A \to C, B \to C\)\(1 + 2 + 4 = 7\) 个独立参数;
  • 完全不连接模型:\(1 + 1 + 1 = 3\) 个独立参数。

图 2.21 展示了随数据量增加时全连接模型与不连接模型的贝叶斯评分相对于真实模型的变化情况。在图中我们用真实模型的评分做了减法,因此 0 以上的值表示"该模型相对真实模型在给定数据下是更好的表示"。图示表明:当样本少于 5000 时,不连接模型优于真实模型;全连接模型始终不比真实模型更好,但当样本量达到约 10000 时,它开始优于不连接模型——因为此时已有足够数据来充分估计其 7 个独立参数。

可能贝叶斯网络结构的数量是超指数级增长的。10 个节点就有约 \(4.2 \times 10^{18}\) 种有向无环图;20 个节点则有 \(2.4 \times 10^{72}\) 种。除了节点数很少的贝叶斯网络,我们无法枚举所有可能结构来寻找最高评分的网络。因此我们必须依赖某种搜索策略。幸运的是,搜索是一个通用问题,多年来已有大量不同的通用搜索算法被研究过。

最常见的搜索策略之一是K2 算法(之所以这样命名,是因为它是一个被称为 Kutató 的系统的演化版本)。K2 搜索(算法 2.8)以多项式时间运行,但不能保证找到全局最优的网络结构。它可以使用任何评分函数 \(f\),但常常与贝叶斯评分一起使用——原因是后者能在"模型复杂度"与"可用数据量"之间取得平衡。K2 从一张没有任何有向边的图开始,然后按一个假定的变量顺序迭代,以"使评分提升最多"的方式贪婪地给各节点添加父节点。实践中常对任一节点的父节点数加上一个上界以减少所需计算。原 K2 算法假设狄利克雷先验参数 \(\alpha_{ijk} = 1\) 对所有 i、j、k 成立,但原则上可以使用任何先验。

另一种通用搜索策略是局部搜索(local search),有时也称为爬山法(hill climbing)梯度上升(gradient ascent),其梗概如算法 2.9 所示。我们从一张初始图 \(G_0\) 出发,向"评分最高的邻居"移动。图的邻居(neighborhood)由那些"只通过一个基本图操作即可到达"的图构成。基本图操作包括:

  • 若 A 与 B 之间没有边,则引入 \(A \to B\)
  • \(A \to B\),则移除从 A 指向 B 的边;
  • \(A \to B\),则反转边的方向得到 \(A \leftarrow B\)

当然,并非所有操作都可从一张特定图出发进行;那些会向图中引入环的操作是无效的。搜索会持续到当前图对其所有邻居都不再有更高评分。

算法 2.8 有向无环图空间上的 K2 搜索 1. 函数 K2Search(\(f\)) 2. \(X_{1:n} \leftarrow\) 节点顺序 3. \(G_0 \leftarrow\) 包含节点 \(X_{1:n}\) 且无边的图 4. 对 \(i \leftarrow 1\)\(n\) 5. 重复 6. \(G \leftarrow G_0\) 7. 给 \(G_0\) 中的节点 \(X_i\) 添加一个父节点,使得 \(f(G_0)\) 取最大 8. 直到 \(f(G_0) \le f(G)\) 9. 返回 \(G\)

算法 2.9 有向无环图空间上的局部搜索 1. 函数 LocalDirectedGraphSearch(\(f, G_0\)) 2. \(G' \leftarrow G_0\) 3. 重复 4. \(G \leftarrow G'\) 5. \(G' \leftarrow\) \(G\) 的邻居 6. \(G' \leftarrow \arg\max_{G' \in G'} f(G')\) 7. 直到 \(f(G') \le f(G)\) 8. 返回 \(G'\)

局部搜索可能会陷入局部极值,从而无法找到全局最优的网络结构。文献中已提出多种应对局部极值的策略,其中包括:

  • 随机重启(randomized restart):一旦找到一个局部极值,就在搜索空间中随机选择一点重新开始搜索。
  • 模拟退火(simulated annealing):并非总是移动到适应度最大的邻居,而是按某种随机化探索策略访问适应度较低的邻居。随着搜索推进,探索中的随机性按某种调度逐渐减小。这种方法被称为"模拟退火",灵感来自冶金学中的退火过程。
  • 禁忌搜索(tabu search):维护一张包含最近访问过的搜索空间点的禁忌表;搜索算法会避开禁忌表中的邻居。
  • 遗传算法(genetic algorithms):从一个由搜索空间中的随机点(以二进制串表示)组成的初始种群开始。在有向图空间中,串中的一位表示两个节点之间箭头的存在与否。种群中个体按其评分成比例地繁殖。被选中繁殖的个体的串通过遗传交叉(genetic crossover)随机重组。遗传交叉在两个随机选中的个体上选定一个交叉点,然后交换该点之后的串。种群中也通过随机翻转串中的某些位来引入变异(mutation)。这一演化过程持续到搜索空间中找到一个令人满意的点。
  • 模因算法(memetic algorithms):有时也被称为"遗传局部搜索",是遗传算法与局部搜索的组合。遗传重组后,对个体施加局部搜索。

某些搜索策略在特定数据集上可能优于其他策略;但总体上寻找全局最优仍是 NP-hard 的。不过许多应用并不要求找到全局最优的网络结构——一个局部最优的结构通常已经可以接受。

2.4.3 Markov Equivalence Classes

如前所述,贝叶斯网络的结构编码了一组条件独立假设。在学习网络结构时需要观察一个重要现象:不同的两个图可以编码相同的独立假设。举一个简单的例子:考虑两变量网络 \(A \to B\)。该网络所包含的独立假设与 \(A \leftarrow B\) 是相同的——单凭"两节点之间是否有边"是无法确定条件独立性的,必须考察整个图上的 d-分隔结构。仅凭数据无法确定 A 与 B 之间边的方向。这一观察在数据驱动的因果发现中至关重要:它意味着仅靠"观测数据"无法区分因果方向——若不引入额外假设(如"因果马尔可夫条件"加"忠实性"),数据本身无法告诉我们"X 导致 Y"还是"Y 导致 X"。这也是为什么因果推断文献常需要借助"干预实验"或"时间先后"等额外信息来打破这一对称性。这一不可区分性同时也划定了"纯数据驱动结构学习"的能力上限——在只有观测数据的情形下,我们能恢复的至多是"马尔可夫等价类",而不能进一步细化到具体的因果方向。这一上限对工程实践有重要意义:当数据来源只能是被动观测时,结构学习算法不应声称"已经找到了因果方向",而应诚实地说"我们找到了一个马尔可夫等价类"。

我们说两个图马尔可夫等价(Markov equivalent)当且仅当它们编码相同的条件独立假设集。可以证明:两个图马尔可夫等价当且仅当(1)它们的无向边集相同——即忽略方向后,边的存在性一致;(2)它们的不道德 v-结构(immoral v-structures)相同。一个不道德 v-结构是指形如 \(X \to Y \leftarrow Z\) 的 v-结构,且 X 与 Z 不直接相连——这意味着 v-结构是"未被 X 与 Z 之间的边'短路'掉的汇合结构"。马尔可夫等价类是一个包含所有彼此马尔可夫等价的有向无环图的集合。该定理由 Verma 与 Pearl 在 1990 年代证明,是因果发现与结构学习文献中引用最广的结果之一。

一般而言,同一马尔可夫等价类中的两个结构可能被赋予不同的评分。然而,若在使用贝叶斯评分的同时使用满足 \(\kappa = \sum_{jk} \alpha_{ijk}\) 对所有 i 都恒定的狄利克雷先验,那么同一等价类中的两个马尔可夫等价结构会被赋予相同的评分。这种先验被称为 BDe 先验;其特例是 BDeu 先验,它令 \(\alpha_{ijk} = \kappa / (q_i r_i)\)。尽管常用的均匀先验 \(\alpha_{ijk} = 1\) 不总能保证同一等价类中的结构被赋予相同评分,但所赋评分通常彼此相当接近。一种"为同一等价类中的所有结构赋予相同评分"的评分函数被称为评分等价(score equivalent)。评分等价这一性质在工程上很有价值:当我们只关心"哪个等价类最匹配数据"时,不必再为"等价类内部的具体有向图选择"担忧——只要评分等价,我们就可以在"等价类代表"上做搜索,而不必在每个成员上重复计算。

一个马尔可夫等价类可以用部分有向图(partially directed graph) 表示——有时也称为本质图(essential graph)有向无环图模式(DAG pattern)。部分有向图可以同时包含有向边与无向边:其中有向边表示"该边方向在所有等价类成员中保持一致",无向边表示"该边方向在不同成员中可能不同"。图 2.22a 展示了一个表示某马尔可夫等价类的部分有向图。有向无环图 \(G\) 是该部分有向图 \(G'\) 所编码的马尔可夫等价类的成员,当且仅当 \(G\)(1)忽略方向后具有与 \(G'\) 相同的边集、(2)与 \(G'\) 具有相同的 v-结构。图 2.22b、2.22c 是该等价类的成员示例;图 2.22d 不是其成员——这正是因为它引入了一个不在 \(G'\) 中的新 v-结构。

比起搜索有向无环图空间,我们可以搜索由部分有向图表示的马尔可夫等价类空间。马尔可夫等价类的空间比有向无环图的空间小,因此搜索可以更高效。2.4.2 节中介绍的任何通用搜索策略都可被采用。若使用某种形式的局部搜索,我们需要定义"邻居"——即用于刻画"图中邻居"概念的局部图操作,例如:

  • 若 A 与 B 之间没有边,则添加 \(A - B\)\(A \to B\)——前者表示"我们暂不能确定方向",后者表示"我们确定方向";
  • \(A - B\)\(A \to B\),则移除 A 与 B 之间的边;
  • \(A \to B\),则反转边的方向得到 \(A \leftarrow B\)
  • \(A - B - C\),则添加 \(A \to B \leftarrow C\)

这四条操作共同构成了"在部分有向图空间上的局部搜索"的邻域定义。贝叶斯评分是针对有向无环图定义的。要为部分有向图评分,我们需要先生成其马尔可夫等价类的一个成员,并计算其评分。从部分有分图生成一个成员,涉及把无向边转化为有向边,且转化的方式不引入新的 v-结构——这保证了生成出的有向图仍属于原等价类。搜索部分有向图空间的好处不仅在于状态空间更小,更在于"在数据上不可区分的网络"会被自然归入同一等价类,从而避免了"在等价类内部做无效搜索"的时间浪费——这是 2.4.3 节观察的直接工程应用。Chickering 在 2002 年的工作系统化了"在部分有向图空间上的搜索",并证明了在某些条件下该方法能保证收敛到全局最优。实践中,"先生成等价类、再在等价类上选代表"的两步策略比"直接搜索有向无环图空间"要快一个量级——这一加速在大型贝叶斯网络(如 50 节点以上)的结构学习中具有决定性意义。这一做法也是"算法按问题结构定制"的一个范例——它表明,理解问题的数学结构(这里即"马尔可夫等价类")往往比单纯优化搜索算法更能带来性能提升。值得补充的是,部分有向图搜索的实现中常需要"成员生成算法"——给定一个部分有向图 P,输出一个满足"v-结构保持"的有向无环图 G;该子问题的复杂度本身也是 NP-hard,但实际中通过启发式规则(如按拓扑序为每条无向边选择方向)可以快速得到可行解。

2.5 Summary

  • 不确定性可能源于信息的不完全,或源于在预测未来事件时面临的实践或理论限制。两种来源在形式上同构,都用"概率分布"这一语言来刻画。
  • 恰当地刻画不确定性,对于构建稳健的决策系统十分重要——决策理论的一个核心结论是:在不确定下的最优决策必须基于对不确定性的"完整"描述(即整个分布),而不是仅基于点估计。
  • 贝叶斯网络以紧凑方式表示变量上的分布——参数规模从指数级降至与图结构相关的多项式级。
  • 网络的结构编码了条件独立的假设——这正是"贝叶斯网络能用更少参数表示联合分布"的原因所在。
  • 贝叶斯网络是一种灵活的表示,可用于编码多种多样的模型——既包括纯离散模型、纯连续模型,也包括 2.1.6 节介绍的混合模型;既包括静态模型,也包括 2.1.7 节介绍的时序模型。
  • 若能利用网络结构,概率推理可以做到高效——变量消去、连接树、置信传播等算法都把"网络结构"作为削减计算量的关键信息。
  • 贝叶斯方法与最大似然方法都可以用于推断模型参数与结构——前者引入先验以缓解数据稀少时的过拟合,后者则是数据充足时的渐近最优选择。

本章五个主要小节构成一条"由静到动、由已知到学习"的递进:2.1 节讨论表示(贝叶斯网络)、2.2 节讨论推理(精确与近似)、2.3 节讨论参数学习(最大似然与贝叶斯)、2.4 节讨论结构学习(评分与搜索);这四者合起来构成一个完整的"概率模型工具箱"。第 3 章将在此基础上加入"决策节点"与"效用节点",把贝叶斯网络推广为决策网络。

需要强调的是,本章所讨论的"不确定性下的推理"与"不确定性下的决策"是两类不同的问题——前者关心"在给定证据下后验是什么",后者关心"在给定目标函数下最优行动是什么"。本章已经为后者准备好了所有"概率"侧的工具;第 3 章将补全"效用"这一侧。读者若熟悉机器学习,可以把本章的"贝叶斯网络 + 推理 + 学习"对应于"概率图模型"的核心理论,把第 3 章开始的"决策网络 + 效用"对应于"决策理论"的工程化路径——两者在"概率"这一侧完全共享,但在"目标函数"这一侧分道扬镳。本章在 2.2.4 节给出的"贝叶斯推理 NP-hard"这一负面结论并不意味着贝叶斯网络在工程上无用——它只是告诉我们"在最坏情况下需要借助近似"。这一观察贯穿本章:精确方法的失败并不是终章,而是促使我们设计更聪明的近似策略的起点。值得指出的是,本章展示的概率模型工具箱有"模块化"的优势——表示、推理、参数学习、结构学习这四个组件可以独立替换,从而在工程上支持灵活的系统组合。例如,可以"用 GMM 做连续变量表示 + 用置信传播做推理 + 用 MLE 做参数学习 + 用局部搜索做结构学习",这种组合在实践中是常见的。这种模块化设计也意味着,本章每一节都可以单独深入——例如"近似推理"本身就是一个独立的研究方向;但作为入门,本章给出的全景图足以让读者判断"在面对具体问题时应当使用哪种工具"。

2.6 Further Reading

关于概率模型最深入的系统论述之一是 Koller 与 Friedman 合著的 Probabilistic Graphical Models: Principles and Techniques [2]。Barber 在 Bayesian Reasoning and Machine Learning [3] 中也对概率模型及其应用做了综述。Russell 与 Norvig 在其广为使用的人工智能教材 Artificial Intelligence: A Modern Approach [4] 中贯穿使用了贝叶斯网络。

概率论的基础在 Jaynes 的 Probability Theory: The Logic of Science [1] 中讨论。Fishburn 综述了主观概率的公理化 [5],Dupré 与 Tipler 给出了一种更晚近的公理化 [6]。Bertsekas 与 Tsitsiklis 的教材提供了对概率论的全面介绍 [7]。

数本教材讨论了贝叶斯网络以及其他概率图模型(如马尔可夫随机场与因子图)中的推理 [2]–[4], [8]–[10]。这些书讨论了包括置信传播与连接树算法在内的多种推理方法——后者曾在 2.2 节中提及。Kim 与 Pearl 给出了多叉树(polytrees)中精确推理的消息传递算法 [11]。Cooper 证明了贝叶斯网络中的推理是 NP-hard [12]。

对贝叶斯网络结构与参数学习的综述可参见教材 Probabilistic Graphical Models: Principles and Techniques [2] 与 Learning Bayesian Networks [13]。如正文所述,学习最优网络结构是 NP-hard [14], [15]。Cooper 与 Herskovits 发展了算法 2.8 所介绍的 K2 搜索算法 [16]。部分有向图空间上的搜索由 Chickering 讨论 [17]。Heckerman、Geiger 与 Chickering 证明:使用 BDe 先验时,任何两个马尔可夫等价结构具有相同的贝叶斯评分 [18]。BDeu 先验最早由 Buntine 提出 [19]。

本章个人批注

作为全书理论篇的第二章,本章承担了一个非常重要的职能——把"不确定性"从第 1 章的口语化直觉转化为可计算的形式化语言。我个人读下来的体会是:第 1 章是"我们为什么要面对不确定性下的决策",而本章则回答了"我们要用什么样的数学对象来承载这种不确定性"。这一转换看似平淡,实则奠定了后续所有算法(第 3 章的决策网络、第 4 章的 MDP、第 6 章的 POMDP)的共同基础。

2.1 节从信念度的公理化切入,把"概率"放在"对命题的偏好"这个更原始的概念之上,这是贝叶斯学派的标志——也是本书一以贯之的立场。Jaynes 的那条"普遍可比性 + 传递性 → 实值函数"的论证路径在 2.1.1 节被压缩成两行,但它的重要性远超篇幅本身。2.1.2 节关于高斯分布的讨论则相对工程化——特别值得注意的是作者用 JFK 机场终端区高度分布这一具体例子指出"高斯并不总是合适",从而顺势引出 GMM 与离散化两种替代方案。这种"用一个真实反例逼出方法选择"的写法很值得借鉴。

2.1.4–2.1.5 节是本章最具结构性意义的两节:贝叶斯网络以有向无环图作为联合分布的紧凑表示,而图形的边结构又编码了条件独立假设。"链式法则"是前者的实现工具,"d-分隔"则是后者的判别工具。2.1.5 节末尾关于"解释消退"的讨论是一个常被忽视但极具诊断价值的概念——当 B 与 S 共享子节点 E 时,观测到 E = 1 后 B 与 S 不再独立,这种"两个本来独立的因变量在公共结果被观察到后变得条件依赖"的现象在医学诊断、故障诊断、风险归因等领域都有非常具体的应用。

2.2 节的推理部分是本章最具计算色彩的一节。2.2.3 节关于变量消去的介绍让我想到几个值得在后续章节关注的问题:一是"消去顺序"为什么 NP-hard——直觉上是因为中间表的大小可能爆炸;二是在线性动态系统中"为什么高斯分布在滤波后仍然是高斯"——这是卡尔曼滤波可行性的根本。2.2.4 节把"贝叶斯推理是 NP-hard"这件事与 3SAT 显式地连接起来,对我而言是一个意外的收获:以前只是从文献里知道这个结论,但作者把 3SAT 实例显式编码为贝叶斯网络这一构造非常直接,让这个抽象结论变得具体可读。

2.2.5 节对近似推理的讨论围绕"采样"这一主线展开。直接采样、似然加权、吉布斯采样这三个方法并不是平行的——它们是对"似然 vs. 后验"这一张力的不同回应。直接采样忠实于后验但效率差;似然加权人为把样本拉到证据一致但引入权重偏差;吉布斯采样则用 MCMC 兼顾两者但需要老化与稀疏。化学探测的例子(99.9% 的似然加权样本权重为 0.001)非常生动地揭示了似然加权在极端先验下失效的情形。2.3–2.4 节则把视角从"给定结构与参数做推理"切换到"从数据中学习结构与参数"——这两节构成本章的应用面。2.3.1 节用二项分布的最大似然估计作为入门,引出"零计数问题"(数据少时 MLE 会给出 0 或 1 的极端值),从而为 2.3.2 节的贝叶斯参数估计做铺垫。伪计数(pseudocount)这一概念让我想到它在工业实践中常被用作"专家经验编码器"——把不可量化的领域知识转成先验的 \(\alpha\)\(\beta\) 增量,这是一条非常实用的工程经验。

2.4 节的结构学习是本章中最"算法味"最浓的部分。K2 搜索、局部搜索及其改进(随机重启、模拟退火、禁忌搜索、遗传算法、模因算法)一气呵成,构成了一组在不同问题规模下可替换的搜索策略工具箱。马尔可夫等价类(2.4.3)与部分有向图(2.4.4)是本章最具结构洞见的两个小节——它们都指向同一个事实:"数据只能约束到条件独立的程度,而不能约束到边方向的程度"。这一观察在因果推断(causal inference)领域是一条核心原则,在本书后续章节中应该还会反复出现。

整体而言,本章是一份"概率模型导论"——它把信念度公理化、贝叶斯网络、推理算法、参数与结构学习、采样方法、马尔可夫等价类这六个话题压缩在 50 页内,每一节都给出了该话题的最小可用版本并指明了 Further Reading。这种"以一本书的篇幅做地图"的工作方式是 Kochenderfer 教材的标志性风格。

与上下章的衔接(一段话)

本章是全书 12 章中的第 2 章,也是理论篇的第二章——紧接在第 1 章"问题域与方法谱系"之后,为后续 5 章理论章节奠定概率模型这一共同基础。从位置上看:第 1 章给出了"为什么要做不确定性下的决策"以及"五种设计方法";本章则把其中"概率视角"这一支(对应 1.3.2 监督学习小节中提到的"概率化表征")彻底展开。第 3 章 Decision Problems 会在本章的概率模型之上加入"效用(utility)"与"决策节点"的概念,把贝叶斯网络推广为决策网络(decision networks)——这是从"表示不确定性"到"在不确定性下做最优决策"的关键一跳。第 4 章 Sequential Problems 则会引入"行动的结果是概率性的"这一时间维度,把本章的静态贝叶斯网络推广为动态贝叶斯网络(即本章 2.1.7 节末尾已经初步讨论的动态贝叶斯网络的更系统版本)。换言之,本章是后续四章理论章节(决策问题、序贯问题、模型未知、状态观测不完备、多智能体协作)的共同前置——理解了本章的概率表示与推理,第 3 章的决策网络、第 4 章的 MDP、第 6 章的 POMDP 都可以视为在本章的图结构上加入新的节点类型(决策、状态转移、观测)或新的变量类型(效用、奖励、信念)。从作者的整体布局看,本章在理论篇中的位置类似于"先把积木搭好,再让后续章节用这些积木搭出更大的结构"。