"在统计学中，蒙特卡洛方法是一种基于随机采样的计算方法。"
"Monte Carlo methods are a broad class of computational algorithms that rely on repeated random sampling."
—— 统计学经典原理
🐍Python工具箱
numpy.random.normal()正态分布随机数
scipy.stats.multivariate_normal()多元正态分布
numpy.linalg.cholesky()Cholesky分解
matplotlib.pyplot.scatter()散点图
matplotlib.pyplot.plot()轨迹线图
numpy.random.choice()随机选择
本章概览
🎯 核心内容
- • 为什么需要MCMC方法
- • 马尔可夫链基础概念
- • Metropolis-Hastings算法
- • Gibbs采样方法
- • 收敛性诊断
📐 关键概念
- • 马尔可夫链（Markov Chain）
- • 平稳分布（Stationary Distribution）
- • 接受概率（Acceptance Probability）
- • 提议分布（Proposal Distribution）
- • 老化期（Burn-in Period）
22.1
为什么需要MCMC？
Why Do We Need MCMC?
在前面的章节中，我们学习了贝叶斯推断的基本原理。 然而，在实际应用中，我们常常遇到一个关键困难：后验分布无法直接计算。 本节将探讨这个问题，并介绍MCMC方法作为解决方案。
🚫 贝叶斯推断的计算瓶颈
回顾贝叶斯定理的核心公式：
fΘ|X(θ|x) = fX|Θ(x|θ) · fΘ(θ) / fX(x)
后验分布 = (似然 × 先验) / 证据
✅ 容易计算
似然 fX|Θ(x|θ)
通常是已知的概率模型
通常是已知的概率模型
✅ 容易指定
先验 fΘ(θ)
根据经验或无信息先验设定
根据经验或无信息先验设定
❌ 计算困难
证据 fX(x)
需要对θ积分，通常无解析解
需要对θ积分，通常无解析解
💥 证据因子的积分问题
fX(x) = ∫ fX|Θ(x|θ) · fΘ(θ) dθ
困难原因：
• θ可能是高维向量（d维）
• 积分区域复杂（如θ∈ℝd）
• 被积函数没有简单形式
• 解析积分不可行
• θ可能是高维向量（d维）
• 积分区域复杂（如θ∈ℝd）
• 被积函数没有简单形式
• 解析积分不可行
📋 具体例子：何时会遇到计算困难？
1️⃣高维参数空间
场景：多元回归模型
参数：θ = (β₀, β₁, ..., βd, σ²)
困难：需要计算 d+2 维积分
例如：d=10时，需要12维积分！
参数：θ = (β₀, β₁, ..., βd, σ²)
困难：需要计算 d+2 维积分
例如：d=10时，需要12维积分！
2️⃣复杂似然函数
场景：混合模型（Mixture Model）
似然：Σk πk fk(x|θk)
困难：求和与积分交织
问题：无法交换顺序简化
似然：Σk πk fk(x|θk)
困难：求和与积分交织
问题：无法交换顺序简化
3️⃣非共轭先验
场景：使用了非共轭先验
例如：Logistic回归的先验
问题：后验不是已知分布族
结果：无法用Beta、Gamma等已知分布表示
例如：Logistic回归的先验
问题：后验不是已知分布族
结果：无法用Beta、Gamma等已知分布表示
4️⃣分层模型
场景：多层次贝叶斯模型
结构：θ₃ → θ₂ → θ₁ → x
积分：需要多重积分边缘化
复杂度：层数增加，计算爆炸
结构：θ₃ → θ₂ → θ₁ → x
积分：需要多重积分边缘化
复杂度：层数增加，计算爆炸
💡 共同特点：在所有这些情况下，虽然我们知道后验∝似然×先验的正比关系， 但归一化常数（证据因子）无法解析计算，导致无法直接从后验分布采样。
🔄 之前的解决方案回顾
在前面的章节中，我们遇到过类似问题，使用了以下方法：
✅ 有效但有限的方法
1. 共轭先验
• Ch20：Beta先验+Bernoulli似然
• 后验也是Beta分布
• 利用Beta函数性质避免积分
• 局限：仅对特定模型组合有效
• 后验也是Beta分布
• 利用Beta函数性质避免积分
• 局限：仅对特定模型组合有效
2. MAP点估计
• Ch18-19：分类中的MAP决策
• 只需比较后验相对大小
• 不需要归一化常数
• 局限：丢失了分布的完整信息
• 只需比较后验相对大小
• 不需要归一化常数
• 局限：丢失了分布的完整信息
❌ 无法应对的情况
1. 需要完整后验分布
• 不确定性量化
• 置信区间估计
• 预测分布计算
• 决策理论应用
• 置信区间估计
• 预测分布计算
• 决策理论应用
2. 复杂模型
• 高维参数空间
• 非共轭先验
• 分层结构模型
• 混合分布
• 非共轭先验
• 分层结构模型
• 混合分布
🎯 核心需求
我们需要一种方法，能够：
1. 不需要计算证据因子（归一化常数）
2. 生成后验分布的样本
3. 适用于任意复杂模型
4. 保证采样的准确性
1. 不需要计算证据因子（归一化常数）
2. 生成后验分布的样本
3. 适用于任意复杂模型
4. 保证采样的准确性
💡 MCMC：问题的解决方案
MCMC (Markov Chain Monte Carlo) 方法巧妙地解决了这个问题。
核心思想
不直接计算后验分布，而是从后验分布中采样！
只要能生成足够多的样本 θ₁, θ₂, ..., θN ~ fΘ|X(θ|x)， 就可以用这些样本来近似后验分布的任何性质。
🎲 Monte Carlo部分
随机采样方法
• 通过大量样本估计
• 样本越多，估计越准
• 法则：大数定律、中心极限定理
• 回顾Ch15：蒙特卡洛基础
• 通过大量样本估计
• 样本越多，估计越准
• 法则：大数定律、中心极限定理
• 回顾Ch15：蒙特卡洛基础
⛓️ Markov Chain部分
构造特殊的随机游走
• 当前样本依赖上一个
• 形成马尔可夫链
• 平稳分布=目标后验
• 长期运行收敛到后验
• 当前样本依赖上一个
• 形成马尔可夫链
• 平稳分布=目标后验
• 长期运行收敛到后验
✨ 神奇之处
1. 只需要后验∝似然×先验
不需要知道归一化常数！
不需要知道归一化常数！
2. 适用于任意复杂模型
高维、非共轭、分层模型都可以
高维、非共轭、分层模型都可以
3. 理论保证收敛
在一定条件下，必然收敛到真实后验
在一定条件下，必然收敛到真实后验
🎨 直观类比：探索未知地形
想象你在一片浓雾中探索山地，想了解地形的高度分布：
❌ 传统方法
方法：画出整个地形图
需要：测量每个点的高度
（计算所有θ的后验概率）
（计算所有θ的后验概率）
问题：
• 点太多，测不过来
• 高维空间更困难
• 点太多，测不过来
• 高维空间更困难
✅ MCMC方法
方法：随机游走采样
策略：
• 从一个点开始
• 向高处走的概率大
• 向低处走的概率小
• 记录访问过的点
• 从一个点开始
• 向高处走的概率大
• 向低处走的概率小
• 记录访问过的点
结果：
高海拔区域访问次数多
→ 样本分布反映地形
高海拔区域访问次数多
→ 样本分布反映地形
🗺️ 类比对应关系
地形高度 ↔ 后验概率
位置 ↔ 参数θ
随机游走 ↔ 马尔可夫链
访问记录 ↔ MCMC样本
🔧 有了样本，我们能做什么？
假设我们通过MCMC获得了N个后验样本：θ₁, θ₂, ..., θN ~ fΘ|X(θ|x)
📊
1. 点估计
后验均值: E[θ|x] ≈ (1/N) Σ θᵢ
后验中位数: median(θ₁, ..., θN)
MAP估计: argmax f(θ|x) ≈ 样本密度最大处
📏
2. 不确定性量化
后验标准差：std(θ₁, ..., θN)
置信区间：2.5%和97.5%分位数 → 95% CI
最高密度区间：HPD (Highest Posterior Density)
📈
3. 后验分布可视化
直方图：样本的频率分布
核密度估计：平滑的概率密度曲线
轨迹图：观察采样过程和收敛性
🔮
4. 预测分布
f(xnew|x) = ∫ f(xnew|θ) f(θ|x) dθ
近似：(1/N) Σ f(xnew|θᵢ)
对每个样本θᵢ生成预测，再平均
🎯
5. 假设检验
后验概率：P(θ > 0 | x) ≈ #{{θᵢ > 0} / N
贝叶斯因子：比较不同模型
决策分析：基于后验的最优决策
🎯 22.1 关键要点
计算瓶颈
• 证据因子积分困难
• 高维参数空间
• 非共轭先验
• 复杂模型结构
• 高维参数空间
• 非共轭先验
• 复杂模型结构
MCMC方案
• 不计算归一化常数
• 通过采样近似后验
• 马尔可夫链+蒙特卡洛
• 理论保证收敛
• 通过采样近似后验
• 马尔可夫链+蒙特卡洛
• 理论保证收敛
样本用途
• 点估计（均值/中位数）
• 置信区间
• 分布可视化
• 预测和决策
• 置信区间
• 分布可视化
• 预测和决策
💡 核心洞察
• 核心困难：复杂模型的后验分布无法直接计算
• 解决思路：用采样代替解析计算
• 关键创新：只需后验正比关系，不需归一化
• 实际价值：使得复杂贝叶斯模型成为可能
• 下一步：学习如何构造这样的马尔可夫链（22.2-22.3）
• 解决思路：用采样代替解析计算
• 关键创新：只需后验正比关系，不需归一化
• 实际价值：使得复杂贝叶斯模型成为可能
• 下一步：学习如何构造这样的马尔可夫链（22.2-22.3）
✦
22.2
Metropolis-Hastings算法
The Metropolis-Hastings Algorithm
在22.1节中，我们了解了为什么需要MCMC方法。 本节将介绍MCMC的第一个经典算法：Metropolis-Hastings (MH) 算法。 这是一个优雅而强大的采样方法，至今仍在广泛使用。
🎯 算法核心思想
目标：从目标分布 π(θ)（通常是后验分布 fΘ|X(θ|x)）中采样
基本策略
1. 提议（Proposal）：
从当前位置θt出发，用提议分布q(θ'|θt)生成候选点θ'
从当前位置θt出发，用提议分布q(θ'|θt)生成候选点θ'
2. 评估（Evaluation）：
计算接受概率α，决定是否接受候选点
计算接受概率α，决定是否接受候选点
3. 决策（Decision）：
以概率α接受θ'作为θt+1，否则保留θt
以概率α接受θ'作为θt+1，否则保留θt
✅ 关键优势
• 只需要π(θ)的正比形式
• 不需要归一化常数
• 自动调整接受率
• 理论保证收敛
• 不需要归一化常数
• 自动调整接受率
• 理论保证收敛
🎨 直观理解
• 像在山上随机游走
• 倾向于往高处走
• 但不是完全贪婪
• 保持足够的探索性
• 倾向于往高处走
• 但不是完全贪婪
• 保持足够的探索性
📝 算法详细步骤
0
初始化
• 选择初始值 θ₀（可以是随机的）
• 设置迭代次数 N
• 选择提议分布 q(θ'|θt)
• 设置迭代次数 N
• 选择提议分布 q(θ'|θt)
1
生成候选点
θ' ~ q(θ'|θt)
从提议分布中采样一个候选值θ'
2
计算接受概率
α = min{1, [π(θ') · q(θt|θ')] / [π(θt) · q(θ'|θt)]}
注意：π(θ)只需正比形式，归一化常数会约掉！
3
接受/拒绝决策
• 生成随机数 u ~ Uniform(0, 1)
• 如果 u < α：接受，令 θt+1 = θ'
• 否则：拒绝，令 θt+1 = θt（保持不变）
• 如果 u < α：接受，令 θt+1 = θ'
• 否则：拒绝，令 θt+1 = θt（保持不变）
4
重复迭代
令 t = t + 1，返回步骤1，直到 t = N
🎯 最终结果
得到样本序列：θ₀, θ₁, θ₂, ..., θN
去除前B个"老化"样本后，剩余样本近似服从目标分布π(θ)
去除前B个"老化"样本后，剩余样本近似服从目标分布π(θ)
🔍 接受概率深入解析
接受概率α是MH算法的核心，让我们详细理解它：
α = min{1, R}
其中 R = [π(θ') / π(θt)] × [q(θt|θ') / q(θ'|θt)]
公式(1)：MH接受概率
第一项：目标分布比
π(θ') / π(θt)
含义：
• 新点相对于旧点的"好坏"
• 如果θ'更好（π(θ')更大），比值> 1
• 更倾向于接受更好的点
• 新点相对于旧点的"好坏"
• 如果θ'更好（π(θ')更大），比值> 1
• 更倾向于接受更好的点
第二项：提议分布比（校正项）
q(θt|θ') / q(θ'|θt)
含义：
• 校正提议分布的不对称性
• 如果对称（q(a|b)=q(b|a)），比值=1
• 保证细致平衡条件
• 校正提议分布的不对称性
• 如果对称（q(a|b)=q(b|a)），比值=1
• 保证细致平衡条件
情况1：R ≥ 1（新点更好）
• α = min{1, R} = 1
• 总是接受新点θ'
• 符合直觉：往"山上"走
• 总是接受新点θ'
• 符合直觉：往"山上"走
情况2：R < 1（新点较差）
• α = R < 1
• 以概率R接受新点θ'
• 不完全贪婪：允许往"山下"走
• 保持探索能力，避免困在局部
• 以概率R接受新点θ'
• 不完全贪婪：允许往"山下"走
• 保持探索能力，避免困在局部
💡 为什么这样设计？
1. 细致平衡条件（Detailed Balance）：
π(θ) · P(θ→θ') = π(θ') · P(θ'→θ)
保证马尔可夫链的平稳分布是π(θ)
π(θ) · P(θ→θ') = π(θ') · P(θ'→θ)
保证马尔可夫链的平稳分布是π(θ)
2. 平衡探索与利用：
既要探索高概率区域（利用），也要保持足够的随机性（探索）
既要探索高概率区域（利用），也要保持足够的随机性（探索）
🎲 提议分布的选择
提议分布q(θ'|θt)的选择对算法性能有重要影响。
🚶
1. 随机游走（Random Walk）
θ' = θt + ε, ε ~ N(0, σ²I)
特点：
• 对称提议：q(θ'|θt) = q(θt|θ')
• 接受概率简化：α = min{1, π(θ')/π(θt)}
• 最常用的选择
• 对称提议：q(θ'|θt) = q(θt|θ')
• 接受概率简化：α = min{1, π(θ')/π(θt)}
• 最常用的选择
参数选择：
• σ太小：步长小，探索慢，但接受率高
• σ太大：步长大，探索快，但接受率低
• 目标接受率：约23%-50%（经验法则）
• σ太小：步长小，探索慢，但接受率高
• σ太大：步长大，探索快，但接受率低
• 目标接受率：约23%-50%（经验法则）
🎯
2. 独立采样器（Independence Sampler）
θ' ~ q(θ'), 不依赖θt
特点：
• 候选点独立于当前位置
• 需要好的提议分布（接近目标分布）
• 不常用，除非有先验知识
• 候选点独立于当前位置
• 需要好的提议分布（接近目标分布）
• 不常用，除非有先验知识
🔧
3. 自适应提议（Adaptive Proposals）
思想：
根据已有样本，动态调整提议分布
根据已有样本，动态调整提议分布
例子：
• 用前期样本估计协方差矩阵Σ
• 用 θ' = θt + ε, ε ~ N(0, cΣ) 作为提议
• c是调整因子（如c=2.38²/d）
• 用前期样本估计协方差矩阵Σ
• 用 θ' = θt + ε, ε ~ N(0, cΣ) 作为提议
• c是调整因子（如c=2.38²/d）
优势：自动适应参数空间的形状
💡 实践建议
• 默认选择：对称随机游走（简单有效）
• 调参策略：通过试验调整步长σ，目标接受率30%-50%
• 多维情况：考虑参数间的相关性，使用协方差矩阵
• 高级技巧：自适应方法（但需注意理论条件）
• 调参策略：通过试验调整步长σ，目标接受率30%-50%
• 多维情况：考虑参数间的相关性，使用协方差矩阵
• 高级技巧：自适应方法（但需注意理论条件）
📊 简单例子：一维正态分布
让我们用一个简单例子来具体演示MH算法的工作过程。
目标
从标准正态分布 N(0, 1) 中采样
π(θ) = (1/√(2π)) exp(-θ²/2)
注意：这里我们故意用MCMC来采样一个简单分布，
目的是验证算法的正确性（实际中N(0,1)可以直接采样）
目的是验证算法的正确性（实际中N(0,1)可以直接采样）
设置
初始值：θ₀ = 0
提议分布：θ' = θt + ε, ε ~ N(0, 1)
迭代次数：N = 10,000
执行步骤（示例迭代）
迭代1：
• 当前：θ₀ = 0
• 提议：θ' = 0 + 0.5 = 0.5（假设ε=0.5）
• 计算：π(0.5)/π(0) = exp(-0.5²/2)/exp(0) = 0.882
• 对称提议，α = min{1, 0.882} = 0.882
• 假设u=0.6 < 0.882，接受！θ₁ = 0.5
• 当前：θ₀ = 0
• 提议：θ' = 0 + 0.5 = 0.5（假设ε=0.5）
• 计算：π(0.5)/π(0) = exp(-0.5²/2)/exp(0) = 0.882
• 对称提议，α = min{1, 0.882} = 0.882
• 假设u=0.6 < 0.882，接受！θ₁ = 0.5
迭代2：
• 当前：θ₁ = 0.5
• 提议：θ' = 0.5 + 1.2 = 1.7（假设ε=1.2）
• 计算：π(1.7)/π(0.5) = exp(-1.7²/2)/exp(-0.5²/2) = 0.086
• α = min{1, 0.086} = 0.086
• 假设u=0.3 > 0.086，拒绝！θ₂ = 0.5（保持）
• 当前：θ₁ = 0.5
• 提议：θ' = 0.5 + 1.2 = 1.7（假设ε=1.2）
• 计算：π(1.7)/π(0.5) = exp(-1.7²/2)/exp(-0.5²/2) = 0.086
• α = min{1, 0.086} = 0.086
• 假设u=0.3 > 0.086，拒绝！θ₂ = 0.5（保持）
预期结果
• 样本的直方图应该近似正态分布形状
• 样本均值 ≈ 0
• 样本标准差 ≈ 1
• 接受率约40%-60%（取决于步长）
• 样本均值 ≈ 0
• 样本标准差 ≈ 1
• 接受率约40%-60%（取决于步长）
⚙️ 实践中的考虑因素
1. 老化期（Burn-in Period）
问题：初始样本可能不代表目标分布
解决：丢弃前B个样本（如B=N/10）
如何判断：观察轨迹图，看是否"稳定下来"
2. 样本相关性（Autocorrelation）
问题：连续样本不独立，θt和θt+1相关
影响：有效样本量小于名义样本量
应对：
• 增加总样本量
• 稀疏化（Thinning）：每k个取一个
• 改进提议分布（减少相关性）
• 增加总样本量
• 稀疏化（Thinning）：每k个取一个
• 改进提议分布（减少相关性）
3. 收敛诊断
方法1：轨迹图（Trace Plot）
观察样本序列是否在合理范围内稳定游走
观察样本序列是否在合理范围内稳定游走
方法2：Gelman-Rubin统计量
运行多条链，比较链间和链内方差
运行多条链，比较链间和链内方差
方法3：有效样本量（ESS）
考虑自相关性后的实际独立样本数
考虑自相关性后的实际独立样本数
4. 多模态问题
挑战：目标分布有多个峰（模态）
困难：单条链可能困在一个模态中
解决：
• 多条链（不同初始值）
• 平行回火（Parallel Tempering）
• 使用更激进的提议分布
• 多条链（不同初始值）
• 平行回火（Parallel Tempering）
• 使用更激进的提议分布
⚖️ MH算法的优势与局限
✅ 优势
1. 通用性强
适用于几乎任何目标分布，
只需知道正比形式
只需知道正比形式
2. 理论完善
在一定条件下，
理论保证收敛到目标分布
理论保证收敛到目标分布
3. 实现简单
算法逻辑清晰，
代码实现容易
代码实现容易
4. 灵活性高
提议分布可以灵活选择，
适应不同问题
适应不同问题
⚠️ 局限
1. 调参困难
提议分布的选择对性能影响大，
需要经验和试错
需要经验和试错
2. 收敛慢
高维空间中，
可能需要大量迭代
可能需要大量迭代
3. 样本相关
连续样本高度相关，
有效样本量减少
有效样本量减少
4. 诊断困难
难以确定是否真正收敛，
特别是多模态情况
特别是多模态情况
💡 实践建议
MH算法是MCMC的"瑞士军刀"——通用但需要技巧。
对于特定问题（如有条件独立结构），
可以考虑更专门的方法，如Gibbs采样（22.3节）。
对于特定问题（如有条件独立结构），
可以考虑更专门的方法，如Gibbs采样（22.3节）。
🎯 22.2 关键要点
算法流程
• 提议候选点
• 计算接受概率
• 接受/拒绝决策
• 迭代重复
• 计算接受概率
• 接受/拒绝决策
• 迭代重复
核心公式
• α = min{1, R}
• R = (π比) × (q比)
• 对称提议简化
• 细致平衡条件
• R = (π比) × (q比)
• 对称提议简化
• 细致平衡条件
实践要点
• 老化期处理
• 提议分布选择
• 收敛性诊断
• 接受率监控
• 提议分布选择
• 收敛性诊断
• 接受率监控
💡 核心洞察
• 关键创新：接受概率的巧妙设计，只需正比关系
• 平衡艺术：探索（接受较差点）vs 利用（倾向好的点）
• 提议重要：提议分布的选择直接影响效率
• 实用技巧：随机游走+对称提议是最常用组合
• 下一步：学习Gibbs采样，处理高维条件独立问题（22.3）
• 平衡艺术：探索（接受较差点）vs 利用（倾向好的点）
• 提议重要：提议分布的选择直接影响效率
• 实用技巧：随机游走+对称提议是最常用组合
• 下一步：学习Gibbs采样，处理高维条件独立问题（22.3）
✦
22.3
Gibbs采样方法
Gibbs Sampling Algorithm
在22.2节中，我们学习了通用的Metropolis-Hastings算法。 本节将介绍另一个重要的MCMC算法：Gibbs采样。 它是MH算法的一个特例，但在高维条件独立问题中特别高效。
💡 核心思想
关键洞察：高维联合分布难以直接采样，但条件分布往往简单得多！
Gibbs采样策略
问题：从联合分布 p(θ₁, θ₂, ..., θd) 采样困难
观察：条件分布 p(θᵢ | θ₁, ..., θᵢ₋₁, θᵢ₊₁, ..., θd) 容易采样
策略：每次只更新一个变量，固定其他变量
从该变量的条件分布中采样
从该变量的条件分布中采样
✅ 巨大优势
• 不需要提议分布
• 接受率始终为1
• 实现简单直接
• 特别适合条件独立结构
• 接受率始终为1
• 实现简单直接
• 特别适合条件独立结构
🎯 关键条件
• 需要知道所有条件分布
• 条件分布必须易于采样
• 适用于分层/图模型
• 变量有条件独立性
• 条件分布必须易于采样
• 适用于分层/图模型
• 变量有条件独立性
📝 Gibbs采样算法步骤
假设目标分布为 p(θ₁, θ₂, ..., θd)，有d个变量
0
初始化
• 选择初始值 θ⁽⁰⁾ = (θ₁⁽⁰⁾, θ₂⁽⁰⁾, ..., θd⁽⁰⁾)
• 设置迭代次数 N
• 设置迭代次数 N
1
更新第1个变量
θ₁⁽ᵗ⁺¹⁾ ~ p(θ₁ | θ₂⁽ᵗ⁾, θ₃⁽ᵗ⁾, ..., θd⁽ᵗ⁾)
从第1个变量的条件分布采样，固定其他变量
2
更新第2个变量
θ₂⁽ᵗ⁺¹⁾ ~ p(θ₂ | θ₁⁽ᵗ⁺¹⁾, θ₃⁽ᵗ⁾, ..., θd⁽ᵗ⁾)
注意：使用刚更新的θ₁⁽ᵗ⁺¹⁾（而非θ₁⁽ᵗ⁾）
⋮
依次更新其他变量
继续更新 θ₃, θ₄, ..., θd
每次都使用最新更新的值
每次都使用最新更新的值
d
更新第d个变量
θd⁽ᵗ⁺¹⁾ ~ p(θd | θ₁⁽ᵗ⁺¹⁾, θ₂⁽ᵗ⁺¹⁾, ..., θd-1⁽ᵗ⁺¹⁾)
完成一轮完整的扫描（sweep）
↻
重复迭代
令 t = t + 1，返回步骤1，直到 t = N
🎯 关键特点
• 每个变量的更新自动被接受（接受率=100%）
• 使用最新的值（Gauss-Seidel风格）
• 一次完整扫描称为一个迭代（iteration）
• 最终得到样本序列：θ⁽⁰⁾, θ⁽¹⁾, θ⁽²⁾, ..., θ⁽ᴺ⁾
• 使用最新的值（Gauss-Seidel风格）
• 一次完整扫描称为一个迭代（iteration）
• 最终得到样本序列：θ⁽⁰⁾, θ⁽¹⁾, θ⁽²⁾, ..., θ⁽ᴺ⁾
🔄 Gibbs采样 vs Metropolis-Hastings
Gibbs采样实际上是MH算法的特例，让我们看看它们的关系：
| 特征 | Metropolis-Hastings | Gibbs采样 |
|---|---|---|
| 提议分布 | 需要选择 q(θ'|θₜ) 需要调参 | 使用条件分布 无需选择 |
| 接受率 | α = min{1, R} 可能拒绝候选点 | α = 1 总是接受 |
| 更新方式 | 同时更新所有变量 θ⁽ᵗ⁾ → θ⁽ᵗ⁺¹⁾ | 逐个更新变量 θᵢ⁽ᵗ⁾ → θᵢ⁽ᵗ⁺¹⁾ |
| 适用范围 | 通用 任何目标分布 | 需要条件分布 适合特定结构 |
| 效率 | 取决于提议分布 可能较低 | 条件独立时高 强相关时低 |
| 实现难度 | 需要调参和测试 中等 | 如果知道条件分布 非常简单 |
💡 理论联系
Gibbs采样是MH算法的特例，其中提议分布选为条件分布：
q(θ' | θₜ) = p(θᵢ' | θ₋ᵢ⁽ᵗ⁾)
此时接受概率α恰好等于1（这可以通过细致平衡条件证明）
q(θ' | θₜ) = p(θᵢ' | θ₋ᵢ⁽ᵗ⁾)
此时接受概率α恰好等于1（这可以通过细致平衡条件证明）
📊 经典例子：二元正态分布
这是Gibbs采样最经典的例子，帮助理解算法的工作机制。
目标分布
从二元正态分布采样：(θ₁, θ₂) ~ N(μ, Σ)
μ = [0, 0]ᵀ, Σ = [[1, ρ], [ρ, 1]]
ρ是相关系数（如ρ=0.8表示强正相关）
条件分布（关键！）
条件分布1：
θ₁ | θ₂ ~ N(ρθ₂, 1-ρ²)
θ₁ | θ₂ ~ N(ρθ₂, 1-ρ²)
条件分布2：
θ₂ | θ₁ ~ N(ρθ₁, 1-ρ²)
θ₂ | θ₁ ~ N(ρθ₁, 1-ρ²)
神奇之处：
二元联合分布复杂，但每个条件分布都是简单的一元正态！
二元联合分布复杂，但每个条件分布都是简单的一元正态！
Gibbs采样步骤
初始化：θ⁽⁰⁾ = (0, 0)
迭代t：
1. 采样 θ₁⁽ᵗ⁺¹⁾ ~ N(ρθ₂⁽ᵗ⁾, 1-ρ²)
2. 采样 θ₂⁽ᵗ⁺¹⁾ ~ N(ρθ₁⁽ᵗ⁺¹⁾, 1-ρ²)
1. 采样 θ₁⁽ᵗ⁺¹⁾ ~ N(ρθ₂⁽ᵗ⁾, 1-ρ²)
2. 采样 θ₂⁽ᵗ⁺¹⁾ ~ N(ρθ₁⁽ᵗ⁺¹⁾, 1-ρ²)
重复N次，得到样本 (θ₁⁽ᵗ⁾, θ₂⁽ᵗ⁾), t=1,2,...,N
几何直观
• 样本点沿着坐标轴方向移动
• 第1步：水平移动（更新θ₁）
• 第2步：垂直移动（更新θ₂）
• 形成"之字形"轨迹
• 第1步：水平移动（更新θ₁）
• 第2步：垂直移动（更新θ₂）
• 形成"之字形"轨迹
相关性影响：
• ρ接近0：快速收敛（变量近独立）
• ρ接近1：收敛慢（强相关，之字形明显）
• ρ接近0：快速收敛（变量近独立）
• ρ接近1：收敛慢（强相关，之字形明显）
🎯 何时使用Gibbs采样？
✅ 理想场景
🔗
1. 分层模型
• 贝叶斯分层模型
• 每层参数条件独立
• 条件分布已知
例：随机效应模型
• 每层参数条件独立
• 条件分布已知
例：随机效应模型
📊
2. 图模型
• 马尔可夫随机场
• 贝叶斯网络
• 变量有条件独立性
例：图像去噪
• 贝叶斯网络
• 变量有条件独立性
例：图像去噪
🧩
3. 缺失数据
• EM算法的贝叶斯版本
• 数据增强方法
• 隐变量模型
例：混合模型
• 数据增强方法
• 隐变量模型
例：混合模型
🎲
4. 共轭先验
• 条件后验是标准分布
• 易于采样
• 计算高效
例：正态-逆Gamma模型
• 易于采样
• 计算高效
例：正态-逆Gamma模型
⚠️ 不适合场景
🔴
1. 强相关变量
• 参数高度相关
• 之字形轨迹严重
• 收敛极慢
问题：混合缓慢
• 之字形轨迹严重
• 收敛极慢
问题：混合缓慢
❓
2. 条件分布复杂
• 条件分布无解析形式
• 无法直接采样
• 需要嵌套MCMC
问题：实现复杂
• 无法直接采样
• 需要嵌套MCMC
问题：实现复杂
🚫
3. 约束参数空间
• 参数有复杂约束
• 条件采样困难
• 需要投影或变换
问题：可能违反约束
• 条件采样困难
• 需要投影或变换
问题：可能违反约束
⏱️
4. 离散参数
• 大量离散状态
• 条件分布计算昂贵
• 更新低效
问题：每步计算成本高
• 条件分布计算昂贵
• 更新低效
问题：每步计算成本高
💡 选择建议
如果模型有明显的条件独立结构且条件分布易于采样，
Gibbs采样通常是首选。否则，考虑MH算法或混合使用。
Gibbs采样通常是首选。否则，考虑MH算法或混合使用。
🔧 分块Gibbs采样
改进Gibbs采样效率的一个重要技巧。
核心思想
标准Gibbs：每次更新一个变量
θ₁, θ₂, θ₃, ... 逐个更新
θ₁, θ₂, θ₃, ... 逐个更新
分块Gibbs：同时更新一组变量
(θ₁, θ₂), (θ₃, θ₄, θ₅), ... 按块更新
(θ₁, θ₂), (θ₃, θ₄, θ₅), ... 按块更新
✅ 优势
1. 减少相关性
同时更新相关变量，减少之字形轨迹
同时更新相关变量，减少之字形轨迹
2. 加速收敛
每次迭代信息量更大
每次迭代信息量更大
3. 更好的混合
参数空间探索更充分
参数空间探索更充分
🎯 分块策略
1. 自然分块
按照模型结构分组
按照模型结构分组
2. 相关性分块
把强相关参数放一组
把强相关参数放一组
3. 层次分块
分层模型按层分块
分层模型按层分块
💡 实例
回归模型：θ = (β, σ²)
• 标准Gibbs：依次更新β₁, β₂, ..., βp, σ²
• 分块Gibbs：先更新整个β向量，再更新σ²
• 效果：通常收敛更快，特别是β维度高时
• 标准Gibbs：依次更新β₁, β₂, ..., βp, σ²
• 分块Gibbs：先更新整个β向量，再更新σ²
• 效果：通常收敛更快，特别是β维度高时
💡 实践技巧与陷阱
✅ 技巧1：随机扫描顺序
问题：固定顺序可能引入偏差
改进：每次迭代随机打乱变量更新顺序
效果：理论上保证遍历性，实践中更稳健
✅ 技巧2：过度松弛
思想：在Gibbs采样中引入"动量"
θᵢ⁽ᵗ⁺¹⁾ = 2θ̃ᵢ - θᵢ⁽ᵗ⁾
θ̃ᵢ 是从条件分布采样的值，"过冲"减少相关性
⚠️ 陷阱1：不可约性问题
问题：条件分布可能无法到达某些状态
例子：参数空间有多个不连通的区域
检查：确保链是不可约的（irreducible）
⚠️ 陷阱2：慢混合
问题：强相关导致收敛极慢
表现：之字形轨迹，自相关系数衰减慢
应对：
• 使用分块Gibbs
• 重参数化减少相关性
• 考虑其他MCMC方法
• 使用分块Gibbs
• 重参数化减少相关性
• 考虑其他MCMC方法
⚖️ Gibbs采样的优势与局限
✅ 优势
1. 无需提议分布
自动使用条件分布，
无调参困扰
无调参困扰
2. 接受率100%
每个提议都被接受，
无样本浪费
无样本浪费
3. 实现简单
逻辑清晰直观，
代码易于编写
代码易于编写
4. 理论完善
收敛性质清楚，
有丰富理论支撑
有丰富理论支撑
5. 自然并行化
某些情况下
可并行更新多个变量
可并行更新多个变量
⚠️ 局限
1. 适用范围受限
需要知道所有条件分布，
不如MH通用
不如MH通用
2. 强相关时慢
参数高度相关时，
收敛极慢
收敛极慢
3. 条件分布复杂
有时条件分布
难以采样
难以采样
4. 可能不可约
需要验证
马尔可夫链性质
马尔可夫链性质
5. 诊断仍困难
收敛判断
依然有挑战
依然有挑战
🎯 总结建议
Gibbs采样是特定问题的"利器"——当模型结构合适时，
它比MH更简单高效；但遇到强相关或复杂条件分布时，
需要考虑混合策略或其他高级MCMC方法。
它比MH更简单高效；但遇到强相关或复杂条件分布时，
需要考虑混合策略或其他高级MCMC方法。
🎯 22.3 关键要点
核心思想
• 逐个更新变量
• 使用条件分布
• 接受率100%
• MH的特例
• 使用条件分布
• 接受率100%
• MH的特例
适用场景
• 分层模型
• 图模型
• 条件独立结构
• 共轭先验
• 图模型
• 条件独立结构
• 共轭先验
改进技巧
• 分块采样
• 随机顺序
• 过度松弛
• 重参数化
• 随机顺序
• 过度松弛
• 重参数化
💡 核心洞察
• 精妙简化：用简单的条件分布采样代替复杂的联合分布
• 自动接受：巧妙的设计使得接受率自动为1
• 结构优势：最适合有明确条件独立结构的模型
• 相关陷阱：强相关参数是Gibbs采样的阿喀琉斯之踵
• 实践智慧：分块采样和随机顺序是重要改进技巧
• 自动接受：巧妙的设计使得接受率自动为1
• 结构优势：最适合有明确条件独立结构的模型
• 相关陷阱：强相关参数是Gibbs采样的阿喀琉斯之踵
• 实践智慧：分块采样和随机顺序是重要改进技巧