📚 课程目录第1章课程导入与MCMC全景

1. 课程导入与MCMC全景

在统计建模与机器学习中，我们经常遇到这样一类问题：模型参数的后验分布形式已知，但归一化常数无法解析计算，或者参数维度很高，导致我们无法直接从这个分布中抽取样本。马尔可夫链蒙特卡洛（Markov Chain Monte Carlo，MCMC）正是为解决这类"知道分布形状、却抽不出样本"的困境而生的一类算法。本课程将围绕MCMC的实践展开，从贝叶斯推断的采样需求出发，逐步深入到Metropolis-Hastings、Gibbs采样、收敛诊断与真实案例的完整实现。

为什么需要采样：贝叶斯推断的核心诉求

贝叶斯推断的基本公式是后验正比于似然乘以先验：

p(θ | D) = p(D | θ) · p(θ) / p(D)
其中 p(D) = ∫ p(D | θ) · p(θ) dθ

分母 p(D) 称为边缘似然或证据，它需要对整个参数空间积分。对于简单共轭模型，这个积分有闭式解；但一旦模型包含多层结构、非线性链接函数或高维参数，积分往往不可解。此时后验分布只剩下一个"未归一化的密度函数"——我们能计算任意 θ 处的 p(D|θ)·p(θ)，却不知道它的归一化常数。

贝叶斯推断关心的量几乎都是后验期望：参数均值、可信区间、预测分布、风险函数等，它们都可以写成 E[f(θ)] = ∫ f(θ)p(θ|D)dθ 的形式。如果能够从后验分布中抽取大量样本 θ₁, θ₂, …, θₙ，那么这些期望就可以用样本均值 (1/N)Σf(θᵢ) 来近似。这就是采样方法的价值：把积分问题转化为样本平均问题，绕开归一化常数。

MCMC的核心思想：用马尔可夫链"走"出目标分布

普通蒙特卡洛方法（如逆变换法、拒绝采样）要求能够独立地从目标分布抽样，在高维情形下效率极低甚至不可行。MCMC换了一个思路：构造一条马尔可夫链，使其平稳分布恰好等于目标后验分布。链从某个初始点出发，按转移核逐步游走，经过足够长的"预热期"后，链上的状态就可以近似看作来自目标分布的样本。

这里有两个关键概念需要先建立直觉。其一是细致平稳条件：若转移核 P 满足 π(x)P(x→y) = π(y)P(y→x)，则 π 是该链的平稳分布。Metropolis-Hastings正是通过"提议—接受/拒绝"的机制来构造满足该条件的转移核。其二是遍历性：链需要不可约且正常返，才能保证无论从哪个初始点出发，长期运行后都会收敛到平稳分布。实践中我们无法严格验证遍历性，只能通过多链、不同初值、收敛诊断等手段来增加信心。

需要强调的是，MCMC样本不是独立的。相邻样本之间存在自相关，有效样本量（ESS）通常远小于实际迭代次数。这直接影响后续的误差估计与收敛判断，也是本课程反复强调诊断环节的原因。

典型应用场景

MCMC并非只服务于贝叶斯统计。只要一个问题可以写成"对某个复杂分布求期望"，MCMC就有用武之地。常见的场景包括：

贝叶斯参数估计：线性回归、Logistic回归、泊松回归中，对回归系数与超参数的后验抽样，进而给出系数可信区间与预测分布。
层级/混合模型：随机效应、多层结构模型中参数维度随分组数增长，MCMC（尤其是Gibbs与HMC）是主流工具。
缺失数据与潜变量模型：把缺失值或潜变量当作待采样参数，通过数据增广（data augmentation）在观测与潜变量之间交替采样。
统计物理与组合优化：模拟退火、Ising模型配分函数估计等，本质都是构造马尔可夫链。
贝叶斯神经网络与概率编程：PyMC、Stan、NumPyro等工具底层均以MCMC或变分推断为核心。

一个最小可运行的示例：Metropolis-Hastings采样

为了让"提议—接受"的机制落地，下面给出一个对一维目标分布采样的最小实现。假设目标分布是未归一化的 p(θ) ∝ exp(−θ²/2)·(1+0.3·sin(3θ))，它没有简单的闭式采样方法，但我们可以计算任意点的密度值。

import numpy as np

def target(theta):
    return np.exp(-theta**2 / 2) * (1 + 0.3 * np.sin(3 * theta))

def metropolis_hastings(n_iter=20000, step=1.2, seed=42):
    rng = np.random.default_rng(seed)
    theta = 0.0
    samples = np.empty(n_iter)
    accepted = 0
    for i in range(n_iter):
        proposal = theta + rng.normal(0, step)
        ratio = target(proposal) / target(theta)
        if rng.random() < min(1.0, ratio):
            theta = proposal
            accepted += 1
        samples[i] = theta
    return samples, accepted / n_iter

samples, acc_rate = metropolis_hastings()
print("接受率:", round(acc_rate, 3))
print("后验均值估计:", round(samples[5000:].mean(), 3))

这段代码中有几个实践要点值得注意。第一，step 是提议分布的标准差，它直接决定接受率：步长过小则接受率高但链移动缓慢、自相关强；步长过大则大量提议被拒绝、链停滞。一维情形下接受率落在0.4～0.5附近通常较为理想，高维时理论最优接受率约为0.234。第二，前5000次迭代作为预热期丢弃，是为了让链从初值0.0"忘掉"起点、进入平稳分布区域。第三，samples[5000:].mean() 就是后验均值的蒙特卡洛估计，其精度取决于有效样本量而非总迭代次数。

判断这段示例是否"跑对了"，可以看三件事：接受率是否在合理区间（约0.2～0.5）；迹图是否围绕某个水平带稳定波动而非持续漂移；用不同随机种子或不同初值重复运行，后验均值估计是否接近。这三条正是后续收敛诊断章节的雏形。

本课程的实践路线

本课程共30章，整体遵循"原理—算法—诊断—案例"的递进结构。前段建立贝叶斯采样需求与马尔可夫链基础；中段逐一实现Metropolis-Hastings、Gibbs采样、以及哈密顿蒙特卡洛（HMC）与NUTS等现代采样器，并对比它们在不同维度与相关性下的表现；后段集中处理收敛诊断（迹图、Gelman-Rubin R̂、有效样本量ESS、自相关函数）、调参策略（步长自适应、预热、重参数化）以及真实数据案例（如贝叶斯线性回归、层级模型、缺失数据插补）。每一章都会给出可运行的代码片段、示例参数与结果判读方法，帮助你把"能跑"推进到"可信"。

学习本课程不需要预先精通测度论，但需要熟悉概率分布、条件概率与基础Python编程。建议边读边运行代码，并刻意修改步长、迭代次数、初值等参数，观察链的行为如何变化——MCMC的直觉，很大程度上来自亲手调试这些参数时看到的迹图与接受率。