Bayes Nets -- Sampling(13)

上一讲我们学习了变量消除等算法，它们能精确地计算出概率查询的结果。但我们也知道，精确推理在最坏情况下是NP难的，对于非常庞大和复杂的贝叶斯网络，精确计算仍然是不可行的。

这一讲的核心思想是：当我们无法承担精确计算的巨大代价时，我们可以退而求其次，通过生成大量随机样本来“近似”地估算概率。这种方法被称为近似推理 (Approximate Inference)，而采样 (Sampling) 是实现它的核心技术。

先验采样 (Prior Sampling)和拒绝采样 (Rejection Sampling)

先验采样是最基础、最简单的采样方法, 它直接模拟贝叶斯网络的生成过程来创建一个完整的样本。

算法流程为:

• 按照贝叶斯网络的拓扑顺序（从父节点到子节点）遍历所有变量。

• 对于每个变量 X，根据其条件概率表 P(X|Parents(X))（此时其所有父节点的值都已被采样确定），随机采样一个值赋给 X。

• 重复此过程，直到所有变量都被赋值，就得到了一个完整的样本。

显而易见, 这种方法非常低效，尤其是在处理带证据的查询时。如果我们想计算 P(C|−t)，而 $P(-t) 的概率只有1%。那么，先验采样生成的样本中，大约有99%的样本都会因为 T=+t 而与证据不符，这些样本最终都会被丢弃。为了得到足够多的有效样本，我们需要生成海量的总样本，造成巨大的计算浪费。

而拒绝采样就是对先验采样的一个简单优化。它在生成样本的过程中，一旦发现某个变量的采样值与证据不符，就立刻停止并拒绝这个样本，然后重新开始一个新的。

它避免了为那些注定要被丢弃的样本浪费后续的计算时间。然而, 虽然有所改进，但它仍然没有解决根本问题。在证据本身是小概率事件的情况下，绝大多数的采样尝试仍然会在早期就被拒绝，效率依然很低。

似然加权 (Likelihood Weighting)

这是本次讲义介绍的第一个真正实用的采样算法，它保证了生成的每一个样本都与证据相符。

核心思想：我们不再随机地生成所有变量的值，而是强制将所有证据变量的值固定为观测到的值。但是，这种“强行设定”会扭曲样本的原始概率分布。为了修正这个偏差，我们为每个样本引入一个权重 (weight)。

证据变量指的是值固定的变量, 例如P(S ∣ C)中的C; 非证据变量是需要被采样的变量, 例如上述概率公式中的S。

算法流程（以计算 P(T|+c, +e) 为例）：

1. 初始化：将样本权重 w 设为 1.0。将证据变量 C 和 E 的值固定为 +c 和 +e。

2. 遍历变量（按拓扑序）：

   • 对于非证据变量（如 T 和 S）：像先验采样一样，根据其父节点的值进行随机采样。

   • 对于证据变量（如 C 和 E）：不进行采样。而是，将当前样本的权重 w 乘以该证据变量在其父节点取相应值时的条件概率。

   • 例如，当处理到 C 时，我们将权重更新为 w = w * P(+c | tⱼ)，其中 tⱼ 是 T 在当前样本中的采样值。这个乘法操作的含义是：“在给定父节点 tⱼ 的情况下，我们观测到 +c 这个证据的可能性有多大？”这个可能性就被作为权重的一部分。

3. 最终计算：在收集了足够多的带权样本后，我们不再是简单地计数，而是计算加权计数。例如，要计算 P(+t|+c, +e)，我们将所有 T=+t 的样本的权重相加，然后除以所有样本的总权重。

原理推导

首先我们给出结论，似然加权方法之所以具有一致性，是因为它通过为样本分配权重，巧妙地使加权后的采样分布等同于真实的全联合概率分布。这意味着，通过对这些加权样本进行求和或计算平均值，我们能够得到对真实概率分布的无偏估计。

1. 采样分布 S_WS(z, e)

这个公式描述了似然加权算法生成一个特定样本的概率：

S_WS(z, e) = ∏_jP(z_j|parents(Z_j))

其中： - z：需要采样的变量集合 - e：固定的证据变量集合 - 显然，这两组变量的并集就是网络中的所有变量。

解释：在似然加权中，我们只对非证据变量进行采样。这个公式是所有被采样变量(z_j)的条件概率乘积。例如，如果W是证据变量，那么C,S,R是被采样的变量。这个采样的分布就是P(C) ⋅ P(S|C) ⋅ P(R|C)。

2. 样本权重 w(z, e)

这个公式定义了每个生成样本的权重：

w(z, e) = ∏_kP(e_k|parents(E_k))

解释：这个权重是所有证据变量(e_k)的条件概率的乘积，每个证据变量的条件概率是基于其在当前样本中的父节点值。这个权重衡量了当前生成的样本与已知证据的一致程度。

3. 联合分布的一致性证明

这是证明似然加权正确性的核心：

S_WS(z, e) ⋅ w(z, e) = (∏_jP(z_j|parents(Z_j)))(∏_kP(e_k|parents(E_k)))

解释：当我们将采样分布和样本权重相乘时，我们得到的是所有变量（包括被采样的z_j和证据e_k）的条件概率乘积。这正是贝叶斯网络中全联合概率分布的定义：

P(z, e) = ∏_iP(X_i|parents(X_i))

当这两部分的结果相乘时，它们精确地、完整地重建了原始的全联合概率分布。这就是为什么加权后的采样分布能够无偏地代表真实的联合概率分布，从而保证了算法的一致性。因此，我们可以得出结论：

S_WS(z, e) ⋅ w(z, e) = P(z, e)

示例分析

• 网络结构：C→S, C→R, S→W, R→W

• 查询目标：计算 P(C=T∣S=T,W=T)，即在洒水器打开(S=True)且草地湿润(W=True)的情况下，多云(C=True)的概率

• 证据变量：e={S=T,W=T}

• 非证据变量：Z={C,R}

条件概率表如下:

P(C): | C | P(C) | |:—–:|:——-:| | T | 0.5 |

P(S∣C): | C | P(S=T∣C) | |:—–:|:——-:| | T | 0.1 | | F | 0.5 |

P(R∣C): | C | P(R=T∣C) | |:—–:|:——-:| | T | 0.8 | | F | 0.2 |

P(W∣S,R): | S | R | P(W=T∣S,R) | |:—:|:—:|:————-:| | T | T | 0.99 | | T | F | 0.9 |

下面展示似然加权生成一个样本的全过程:

1. 初始化

• 证据：S=T,W=T
• 样本权重：w=1.0

2. 处理变量 C

• C是非证据变量，需采样
• 根据P(C)采样，结果：C=F
• 当前样本：{C=F}，w=1.0

3. 处理变量 S

• S为证据变量(S=T)
• 更新权重：w = w × P(S=T∣C=F) = 1.0 × 0.5 = 0.5
• 当前样本：{C=F}，w=0.5

4. 处理变量 R

• R是非证据变量，需采样
• 根据P(R∣C=F)采样，结果：R=T
• 当前样本：{C=F,R=T}，w=0.5

5. 处理变量 W

• W为证据变量(W=T)
• 更新权重：w = w × P(W=T∣S=T,R=T) = 0.5 × 0.99 = 0.495
• 最终样本：{C=F,R=T}，w=0.495

下面验证”加权采样分布 = 真实联合概率”：

1. 采样概率SWS(z,e)：
   • SWS = P(C=F) × P(R=T∣C=F) = (1-0.5) × 0.2 = 0.1
2. 样本权重w(z,e)：
   • w = P(S=T∣C=F) × P(W=T∣S=T,R=T) = 0.5 × 0.99 = 0.495
3. 两者相乘：
   • SWS × w = 0.1 × 0.495 = 0.0495
4. 真实联合概率P(z,e)：
   • P(C=F,S=T,R=T,W=T) = P(C=F) × P(S=T∣C=F) × P(R=T∣C=F) × P(W=T∣S=T,R=T)
   • = 0.5 × 0.5 × 0.2 × 0.99 = 0.0495

可以看出验证的结果是两者完全相等.

重复上述过程多次后，使用以下公式计算：

P(C=T∣S=T,W=T) ≈

这个例子展示了似然加权如何通过采样和加权来模拟真实的联合概率分布，从而计算出所需的条件概率。

缺点与不足

似然加权最大的不足之处是上游变量的值不受下游证据的影响。由于采样过程是前向的（从上游到下游），在对某个变量进行采样时，我们无法利用其子节点（下游）的证据信息。

这意味着如果证据变量位于网络中许多独立的下游节点，每个证据都会使最终的权重呈指数级下降。这可能导致权重变得极小，从而引发数值上的不稳定问题。

同时, 当大多数样本的权重都非常小时，最终的概率估计可能完全由一个或少数几个权重异常大的“幸运”样本决定。这会导致估计的方差很高，结果不稳定。

这个局限性意味着似然加权在进行诊断式推理（从结果推断原因）时效率很低。它生成样本的方式和证据是脱节的，只能在最后通过权重来“亡羊补牢”，这导致了大量的计算浪费。

举个栗子

例如, 让我们构建一个非常简单的贝叶斯网络来模拟交通事故。

• 上游变量 (原因): D (司机分心 - Distracted)，可能取值为 T (是) 或 F (否)。

• 下游变量 (结果): A (发生事故 - Accident)，可能取值为 T (是) 或 F (否)。

此时的网络结构为: D→A (司机分心会导致事故)

假设的概率:

• P(D=T)=0.1 (在正常情况下，司机分心的概率是10%)

• P(A=T∣D=T)=0.6 (如果司机分心，发生事故的概率是60%)

• P(A=T∣D=F)=0.01 (如果司机没分心，发生事故的概率是1%)

现在，我们有了一个下游证据：确实发生了一场事故 (A=T)。我们的目标是推断上游原因：司机当时分心的概率是多少？即 P(D=T∣A=T)。

从直觉上，既然事故已经发生了，那么司机分心的可能性应该会远高于平时的10%。

现在我们来看似然加权算法是如何一步步处理这个问题的，并观察其局限性。

第1步：采样上游变量 D : 算法开始生成第一个样本。它遇到的第一个变量是 D。算法会问：“我应该如何为 D 采样？”

由于采样过程是纯粹前向的，它完全看不到下游已经发生的事故证据 A=T。它唯一能依据的就是 D 的先验概率分布 P(D)。

• 有 10% 的概率，它会采样得到 D=T。

• 有 90% 的概率，它会采样得到 D=F。

这就是问题的核心所在！尽管我们知道事故已经发生，这应该让我们更倾向于相信司机分心了，但算法在采样 D 的时候，完全忽略了这个信息。它依然会花费 90% 的精力去生成“司机没有分心”这种与证据非常不符的样本。

第2步：计算权重

在采样完 D 之后，算法才会考虑证据 A=T，并用它来计算权重。

• 情况一 (90%的概率发生): 算法在上一步采样得到 D=F。此时的权重 w=P(A=T∣D=F)=0.01。这是一个极低的权重。

• 情况二 (10%的概率发生): 算法在上一步采样得到 D=T。此时的权重 w=P(A=T∣D=T)=0.6。这是一个相对高很多的权重。

通过上面的过程，我们可以清晰地看到这个局限性体现在哪里：

• 盲目采样：算法在生成样本的“构思”阶段是盲目的。它像一个不知道最终结局的编剧，按照通常的概率写故事的开头（采样上游变量）。

• 效率低下：算法花费了大量的计算资源（90%的样本）去探索那些几乎不可能导致已知证据（事故）发生的情况（司机没分心）。这些样本最终只会得到一个微不足道的权重，对最终结果的贡献极小。

• 结果被“幸运样本”主导：最终的概率估计，几乎完全依赖于那少数（10%）采样到了 D=T 的“幸运样本”。如果样本总数不够多，我们可能采不到足够多的幸运样本，导致最终结果的方差很大，非常不稳定。

吉布斯采样 (Gibbs Sampling)

这是一种完全不同的采样策略，属于马尔可夫链蒙特卡洛 (Markov Chain Monte Carlo - MCMC) 方法的一种。

为了理解吉布斯采样，首先要理解它所属的 MCMC 方法是什么意思: 蒙特卡洛 (Monte Carlo)指代任何依赖于重复随机采样的算法。而马尔可夫链 (Markov Chain)指的是一个状态序列，其中下一个状态的概率只取决于当前状态，而与之前的任何状态都无关

与似然加权那种“一气呵成”的前向采样不同，吉布斯采样的核心思想是迭代式地修正（Iterative Refinement）。它不试图一次性生成一个完美的样本，而是从一个随机的、不完美的完整状态开始, 通过轮流地、逐个地对变量进行重新采样来不断修正这个状态。在经过足够多的修正后，这个状态就会稳定下来，并开始在真实的目标概率分布附近“徘徊”。此时的状态就可以被当作一个有效的样本。

基本的算法流程如下：

1. 初始化：从一个随机的完整状态 x 开始（但证据变量的值是固定的）。

2. 迭代循环：

   • 选择一个非证据变量 Zᵢ。

   • 重新采样：将 Zᵢ 的当前值“清空”，然后根据所有其他变量的当前值，为 Zᵢ 重新采样一个新值。这个新值来自于条件概率分布 P(Z_i|mb(Z_i))，其中 mb(Z_i) 是 Zᵢ 的马尔可夫毯 (Markov Blanket)（即它的父节点、子节点、以及子节点的其他父节点）。

3. 收集样本：每完成一次（或若干次）这样的重采样，当前的状态 x 就被当作一个样本记录下来。

其优势在于, 吉布斯采样在每次重采样时，都会同时考虑“上游”和“下游”的变量（通过马尔可夫毯），因此信息可以在整个网络中双向传播。这使得它在处理似然加权遇到的“下游证据”问题时通常表现得更好。

示例分析

网络结构如下: C → S, C → R, S → W, R → W

查询目标: P(C = T|S = T, W = T)

• 证据变量 (E): S = T, W = T (值是固定的)

• 非证据变量 (Z): C, R (我们需要对它们进行迭代采样)

• 概率表 (CPT): 我们继续使用之前例子中的概率。

吉布斯采样 (Gibbs Sampling) 流程为:

步骤 1: 初始化 - 固定证据: 我们将 S 的值锁定为 T，将 W 的值锁定为 T。

• 随机初始化: 我们为所有非证据变量 C, R 随机赋一个初始值。
  • 比如，我们随机选择 C=F 和 R=F。现在，我们的马尔可夫链的初始状态是：C = F, R = F, S = T, W = T。

步骤 2: 迭代循环: 我们将轮流对非证据变量 C 和 R 进行重新采样。

首先是第 1 轮迭代

Part A: 重新采样变量 C

• “冻结”其他变量: 将其他所有变量的值固定在当前状态：R = F, S = T, W = T。

• 计算条件概率: 我们需要计算 P(C|R = F, S = T, W = T)。

P(C = T|...) ∝ P(C = T, R = F, S = T, W = T)

 = P(C = T) ⋅ P(S = T|C = T) ⋅ P(R = F|C = T) ⋅ P(W = T|S = T, R = F)

 = (0.5) ⋅ (0.1) ⋅ (0.2) ⋅ (0.9) = 0.009

P(C = F|...) ∝ P(C = F, R = F, S = T, W = T)

 = P(C = F) ⋅ P(S = T|C = F) ⋅ P(R = F|C = F) ⋅ P(W = T|S = T, R = F)

 = (0.5) ⋅ (0.5) ⋅ (0.8) ⋅ (0.9) = 0.18

归一化并采样:

P(C = T|...) ≈ 0.009/(0.009 + 0.18) ≈ 0.048 (约 5%)

P(C = F|...) ≈ 0.18/(0.009 + 0.18) ≈ 0.952 (约 95%)

根据这个概率，我们为 C 重新采样。有极大概率我们会采样到 C=F。假设我们这次采样的结果就是 C=F。此时更新状态: 系统的当前状态现在是 C = F, R = F, S = T, W = T (这次 C 的值恰好没变)。

Part B: 重新采样变量 R

• “冻结”其他变量: 将其他变量的值固定在当前最新的状态：C = F, S = T, W = T。

• 计算条件概率: 我们需要计算 P(R|C = F, S = T, W = T)。

P(R = T|...) ∝ P(R = T|C = F) ⋅ P(W = T|S = T, R = T)

 = (0.2) ⋅ (0.99) = 0.198

P(R = F|...) ∝ P(R = F|C = F) ⋅ P(W = T|S = T, R = F)

 = (0.8) ⋅ (0.9) = 0.72

归一化并采样:

P(R = T|...) ≈ 0.198/(0.198 + 0.72) ≈ 0.216 (约 22%)

P(R = F|...) ≈ 0.72/(0.198 + 0.72) ≈ 0.784 (约 78%)

根据这个概率，我们为 R 重新采样。假设这次我们采样的结果是 R=F。然后更新状态: 第1轮迭代结束，系统的最终状态是 C = F, R = F, S = T, W = T。

在之后的第 2 轮迭代 (及后续)中, 我们会重复上面的过程。在第2轮迭代开始时，我们会先重新采样 C，此时它将基于 R = F, S = T, W = T 来计算条件概率（与第一轮相同）。然后，我们会根据新采样的 C 值和固定的证据，再次为 R 重新采样。

随着迭代的进行，状态 C, R 会不断在各种可能的值之间跳转，但跳转的频率会符合它们在证据 S = T, W = T 下的真实后验概率。

步骤 3: 收集样本并得出结论

在经过足够长的”热身”期后，我们开始记录每一轮迭代结束时的状态。例如，我们可能收集到以下10,000个样本（只记录非证据变量）：

C = F, R = F

C = F, R = F

C = F, R = T

…

C = T, R = T

…

最后，我们要回答查询 P(C = T|S = T, W = T)。我们只需要在收集到的10,000个样本中，计算出 C=T 的样本所占的比例即可

理论证明, 只要这个过程重复足够多次，算法生成的样本分布最终会收敛到真实的后验概率分布 P(Q|e)。