蒙特卡洛方法与MCMC

蒙特卡洛方法用随机采样来解决确定性问题——从积分计算到参数估计；MCMC（马尔可夫链蒙特卡洛）则让我们能从任意复杂分布中高效采样。本页梳理从最朴素的拒绝采样到 HMC 的完整谱系。前置阅读：数值优化与数学基础、概率论基础、EM算法与变分推断。

核心直觉

蒙特卡洛的核心思想极其简单：用频率逼近概率，用样本均值逼近积分。

蒙特卡洛积分：你想知道一个湖的面积，但形状不规则。往地图上随机撒 1000 个点，数落在湖里的比例，乘以地图总面积——这就是蒙特卡洛估计。点撒得越多，估计越准（误差按 $1/\sqrt{N}$ 速度下降）。
直接采样：如果你有一个均匀分布的随机数生成器（如骰子），直接用它模拟。但很多分布无法直接采样——比如一个形状奇怪的概率密度函数。
拒绝采样：你想按一个曲线采点，但没有对应工具。那就找一个大盒子把曲线包住，在盒子里均匀撒点，落在曲线下方的留下，曲线上方的拒绝。简单但效率低——如果曲线很窄，大量采样被浪费。
MCMC：聪明的做法：不独立撒点，而是像散步一样在空间里”漫步”。每一步根据当前位置决定下一步去哪，往概率高的地方多走走，往概率低的地方少走走。走得足够久之后，你走过的足迹就是按目标分布采的样。

类比：直接采样是”闭眼撒网”；拒绝采样是”大海捞针”；MCMC 是”跟着人流走”——人流密的地方（高概率区）你就多待，自然就按密度分布采样了。

为什么对 AI 重要？ 蒙特卡洛方法贯穿了现代 AI 的多个核心环节：贝叶斯神经网络的参数推断（用 MCMC 采样后验）、强化学习中的策略评估（用重要采样做 off-policy 评估）、扩散模型的去噪采样（逆向随机微分方程的数值求解）、以及大语言模型的对齐训练（RLHF 中的方差减小采样）。理解蒙特卡洛，就拿到了理解这些前沿技术的底层钥匙。

原理详解

蒙特卡洛积分

我们常需要计算形如下式的期望（即函数 $f$ 在概率分布 $p$ 下的平均值）：

I = \mathbb{E}_{x \sim p}[f(x)] = \int f(x)\, p(x)\, dx

当 $p(x)$ 是高维分布、或者 $f(x)$ 形状复杂时，解析积分几乎不可能。蒙特卡洛方法绕开解析计算：从 $p(x)$ 中采 $N$ 个独立样本 $x_1, x_2, \dots, x_N$ ，用样本均值估计期望：

\hat{I}_N = \frac{1}{N} \sum_{i=1}^{N} f(x_i)

为什么这样做有效？ 根据大数定律，当 $N \to \infty$ 时，样本均值 $\hat{I}_N$ 几乎必然收敛到真实期望 $I$ 。更进一步，根据中心极限定理，估计误差近似服从正态分布：

\hat{I}_N - I \;\approx\; \mathcal{N}\!\left(0,\; \frac{\sigma_f^2}{N}\right)

其中 $\sigma_f^2 = \mathrm{Var}_{x \sim p}[f(x)]$ 是 $f(x)$ 在分布 $p$ 下的方差。因此估计误差的标准差为 $\sigma_f / \sqrt{N}$ 。

维度无关性：这是蒙特卡洛相对网格积分的决定性优势。在一维空间，网格积分用 $N$ 个点能达到 $O(1/N)$ 的精度；但在 $d$ 维空间，同样 $N$ 个点每维只有 $N^{1/d}$ 个网格点，精度退化到 $O(N^{-1/d})$ ——在 100 维空间（深度学习的常见规模）， $N=10^6$ 时每维也只有约 2 个点。而蒙特卡洛的 $O(1/\sqrt{N})$ 收敛速度与维度无关。这也是为什么高维积分（如贝叶斯后验、物理配分函数）几乎只能用蒙特卡洛。

重要采样（Importance Sampling）

当 $f(x)$ 在 $p(x)$ 下出现概率很小但贡献很大时（“稀有事”），直接采样效率极低——绝大多数样本 $f(x_i) \approx 0$ ，只有极少数样本贡献巨大，估计方差 $\sigma_f^2$ 爆炸。重要采样改从一个提议分布 $q(x)$ 中采样（要求 $q(x)>0$ 只要 $p(x)f(x) \neq 0$ ）：

I = \int f(x)\, p(x)\, dx = \int \frac{f(x)\, p(x)}{q(x)}\, q(x)\, dx = \mathbb{E}_{x \sim q}\!\left[\frac{f(x)\, p(x)}{q(x)}\right]

因此从 $q$ 中采 $N$ 个样本后：

\hat{I}_N = \frac{1}{N} \sum_{i=1}^{N} f(x_i)\, \frac{p(x_i)}{q(x_i)} = \frac{1}{N} \sum_{i=1}^{N} f(x_i)\, w_i

权重 $w_i = p(x_i)/q(x_i)$ 修正了”从 $q$ 采但目标是 $p$ “的分布偏差。关键是选好 $q(x)$ ，让它在 $f(x)$ 大的地方采样多。理论上最优提议分布为 $q^*(x) \propto |f(x)|\, p(x)$ ，此时方差为零——但这要求知道我们正想计算的量，实践中只能近似。

强化学习中的关键应用：off-policy 评估（用旧策略收集的数据评估新策略的期望回报）本质就是重要采样。当新旧策略差异大时，权重 $w_i$ 的方差会爆炸，这是离线强化学习的核心难题。详见下文”典型应用”。

马尔可夫链蒙特卡洛（MCMC）

当目标分布 $p(x)$ 只能算到未归一化形式 $\tilde{p}(x)$ （即 $p(x) = \tilde{p}(x) / Z$ ，归一化常数 $Z = \int \tilde{p}(x)\, dx$ 未知且高维下难以计算）时，前面方法都失效：直接采样需要 $Z$ ，拒绝采样需要知道分布的上界。

MCMC 的思路是构造一条马尔可夫链（Markov Chain，一种”下一步只依赖当前状态”的随机过程），使其平稳分布（stationary distribution，即链运行足够久后样本所服从的分布）恰为目标分布 $p(x)$ 。关键工具是细致平衡条件（Detailed Balance）：

p(x)\, T(x \to x') = p(x')\, T(x' \to x)

其中 $T(x \to x')$ 是从状态 $x$ 转移到状态 $x'$ 的概率。直观理解：从 $x$ 跳到 $x'$ 的”流量”等于从 $x'$ 跳回 $x$ 的”流量”，系统达到动态平衡。满足细致平衡就能保证链的平稳分布是 $p(x)$ 。

为什么只要未归一化形式就够？ 注意在细致平衡中，如果我们构造的转移核 $T$ 依赖于比值 $p(x')/p(x) = \tilde{p}(x')/\tilde{p}(x)$ ，那么未知的归一化常数 $Z$ 会在比值中约掉。这就是 MCMC 能避开 $Z$ 的关键。

Metropolis-Hastings 算法

最通用的 MCMC 方法。用一个（通常对称的）提议分布 $q(x' \mid x)$ （如以当前点为中心的高斯）生成候选样本，再按接受率（acceptance probability）决定是否接受：

从 $q(x' \mid x_{\text{current}})$ 采一个候选 $x'$
计算接受率：

\alpha = \min\!\left(1,\; \frac{\tilde{p}(x')\, q(x_{\text{current}} \mid x')}{\tilde{p}(x_{\text{current}})\, q(x' \mid x_{\text{current}})}\right)

以概率 $\alpha$ 接受（ $x_{\text{next}} = x'$ ），否则拒绝（ $x_{\text{next}} = x_{\text{current}}$ ，原地不动）

当提议分布对称时（如高斯 $q(x' \mid x) = q(x \mid x')$ ）， $q$ 项约掉，接受率简化为 $\alpha = \min(1, \tilde{p}(x') / \tilde{p}(x_{\text{current}}))$ ——新样本比当前样本概率高就一定接受，低就按比例概率接受。这个”倾向于往高处走”的简单规则，恰好满足细致平衡。

吉布斯采样（Gibbs Sampling）

当联合分布 $p(x_1, x_2, \dots, x_d)$ 难以直接采样，但每个条件分布 $p(x_j \mid x_{-j})$ （给定其他所有维度时第 $j$ 维的分布）容易采样时（如图模型、主题模型 LDA），逐维度交替采样：

\begin{aligned} x_1^{(t+1)} &\sim p\!\left(x_1 \mid x_2^{(t)}, x_3^{(t)}, \dots, x_d^{(t)}\right) \\ x_2^{(t+1)} &\sim p\!\left(x_2 \mid x_1^{(t+1)}, x_3^{(t)}, \dots, x_d^{(t)}\right) \\ &\vdots \\ x_d^{(t+1)} &\sim p\!\left(x_d \mid x_1^{(t+1)}, x_2^{(t+1)}, \dots, x_{d-1}^{(t+1)}\right) \end{aligned}

每次只采一个维度，用其余维度的当前值（注意用刚更新的最新值）。它是 Metropolis-Hastings 的特例（接受率恒为 1，不会拒绝）。

哈密顿蒙特卡洛（HMC）

Metropolis-Hastings 在高维空间中效率极低：随机提议方向几乎不可能对准高概率区的狭长”山脊”，导致接受率极低（这就是”维度灾难”）。HMC 引入物理直觉：把参数空间想象成一个地形，目标分布的对数 $\log \tilde{p}(x)$ 是势能函数 $U(x) = -\log \tilde{p}(x)$ ，再给每个参数引入一个动量变量 $v$ ，让粒子在地形上按哈密顿动力学滑动：

H(x, v) = U(x) + \frac{1}{2} v^\top M^{-1} v

其中 $M$ 是质量矩阵（通常取单位阵）。模拟时在动量空间采 $v \sim \mathcal{N}(0, M)$ ，然后用蛙跳法（leapfrog integration）按梯度 $\nabla U(x)$ 更新若干步，最后用 Metropolis 接受/拒绝修正数值误差。HMC 能利用梯度信息快速穿越高维空间，大幅减少样本间的自相关——在贝叶斯深度学习中比 MH 效率高几个数量级。

HMC 需要计算 $\nabla \log \tilde{p}(x)$ （相当于反向传播），因此与深度学习框架天然兼容。No-U-Turn Sampler（NUTS） 是 HMC 的自适应版本，自动调步长和路径长度（通过”走到开始回头就停”的启发式），是 PyMC / Stan 的默认采样器。

方差减小技术

朴素蒙特卡洛虽然无偏，但方差可能很大。以下技巧能在不引入偏差的前提下显著降低方差（即同样样本量下精度更高）。

对偶变量（Antithetic Variates）

利用 $f(x)$ 与 $f(-x)$ 或 $f(1-x)$ 的负相关性：每采一个 $u \sim \mathrm{Uniform}(0,1)$ ，同时使用 $u$ 和 $1-u$ ，取平均 $[f(u)+f(1-u)]/2$ 作为一次估计。当 $f$ 单调时两者负相关，方差可降低 50% 以上。

控制变量（Control Variates）

找一个期望 $\mu_C$ 已知、且与 $f(x)$ 正相关的函数 $C(x)$ ，用调整后的估计量：

\hat{I}_{\text{cv}} = \frac{1}{N}\sum_i \left[f(x_i) - \beta\bigl(C(x_i) - \mu_C\bigr)\right]

最优系数 $\beta^* = \mathrm{Cov}(f, C)/\mathrm{Var}(C)$ 。直觉：用 $C$ 的”可预测偏差”来校正 $f$ 的波动。扩散模型推理中的 ODE-based 校正和强化学习中的 baseline 减小方差都是这个思路。

分层采样与准蒙特卡洛（QMC）

把采样空间均匀划分为若干”层”，每层内独立采样（分层采样）；或直接用低差异序列（low-discrepancy sequence，如 Sobol、Halton 序列）代替随机数——这就是准蒙特卡洛（Quasi-Monte Carlo, QMC）。QMC 的收敛速度可提升到 $O((\log N)^d / N)$ ，远好于 $O(1/\sqrt{N})$ 。2025 年以来 QMC 在扩散模型采样中的应用已成为活跃研究方向，详见文末”2025-2026 年进展”。

原理图解

蒙特卡洛方法谱系

MCMC 采样过程

HMC 与 Metropolis-Hastings 的对比

动手实践

NumPy 蒙特卡洛估计圆周率

最经典的入门例子：在边长为 2 的正方形内切一个单位圆，向正方形内均匀撒点，落入圆内的比例趋近于面积比 $\pi/4$ 。

import numpy as np

# 在正方形内撒点，落入内切圆的比例 ≈ π/4
N = 1_000_000
x, y = np.random.uniform(-1, 1, N), np.random.uniform(-1, 1, N)  # 正方形内随机撒点
inside = (x**2 + y**2 < 1).sum()     # 落入单位圆的数量
pi_estimate = 4 * inside / N          # 圆面积/正方形面积 = π/4
# 理论误差 ≈ 4 * sqrt(p(1-p)/N) ≈ 0.001，p = π/4 ≈ 0.785
print(f"π ≈ {pi_estimate:.4f}（误差 ≈ {abs(pi_estimate - np.pi):.4f}）")

用蒙特卡洛估计高维积分

下面展示蒙特卡洛在高维积分上的优势：估算单位超球（ $d$ 维球）内 $\int_{\|x\|\le 1} dx$ 的体积。解析值为 $\pi^{d/2} / \Gamma(d/2+1)$ ，但高维下用网格积分代价随 $d$ 指数增长，而蒙特卡洛代价几乎不变。

import numpy as np
from math import gamma, pi

def mc_sphere_volume(d: int, n_samples: int = 2_000_000) -> float:
    """用蒙特卡洛估计 d 维单位超球的体积。"""
    # 在 [-1,1]^d 的超立方体内均匀采样
    pts = np.random.uniform(-1, 1, size=(n_samples, d))
    inside = (np.sum(pts**2, axis=1) <= 1.0).mean()   # 落入超球的比例
    cube_volume = 2 ** d
    return inside * cube_volume

for d in [2, 5, 10, 20]:
    est = mc_sphere_volume(d)
    true = pi ** (d / 2) / gamma(d / 2 + 1)
    print(f"d={d:>2}: 蒙特卡洛 ≈ {est:10.4f}  真值 = {true:10.4f}")
# 注意：d 很大时真值趋于 0（"高维球体积消失"），蒙特卡洛仍能在固定样本量下给出估计

PyTorch 实现 Metropolis-Hastings 采样

用一个简单的双峰分布演示 MH 算法的工作过程。这个例子帮助理解 MCMC 如何仅用未归一化概率就能采样。

import torch

def unnormalized_pdf(x: torch.Tensor) -> torch.Tensor:
    """双峰混合高斯：未归一化密度 p̃(x) = N(x; -3, 1) + N(x; 3, 1)"""
    comp1 = torch.exp(-0.5 * ((x - 3) / 1) ** 2)
    comp2 = torch.exp(-0.5 * ((x + 3) / 1) ** 2)
    return comp1 + comp2

def metropolis_hastings(n_iter: int = 20_000, burn_in: int = 2_000, seed: int = 0):
    torch.manual_seed(seed)
    x_current = torch.tensor([0.0])                  # 初始点
    samples = torch.zeros(n_iter)
    proposal_std = 2.0                                # 提议高斯的标准差

    for i in range(n_iter):
        # 1) 对称提议：以当前点为中心采高斯候选
        x_proposal = x_current + proposal_std * torch.randn(1)
        # 2) 计算接受率（对称提议 q 项约掉）
        alpha = min(1.0, (unnormalized_pdf(x_proposal) / unnormalized_pdf(x_current)).item())
        # 3) 按概率 α 接受
        if torch.rand(1).item() < alpha:
            x_current = x_proposal
        samples[i] = x_current.item()

    return samples[burn_in:]                          # 丢弃预热期

samples = metropolis_hastings()
print(f"均值 = {samples.mean():.3f}（理论 ≈ 0）")
print(f"标准差 = {samples.std():.3f}（理论 ≈ √(1+9) ≈ 3.16）")

PyMC 贝叶斯推断（MCMC 采样）

实战中极少手写 MCMC，而是用概率编程框架。下面用 PyMC 的 NUTS 采样器估计线性回归的后验：

import pymc as pm
import numpy as np

# 生成模拟数据：线性关系 y = 2x + 1 + 噪声
np.random.seed(42)
x = np.linspace(0, 10, 50)
y = 2 * x + 1 + np.random.randn(50)

# 定义贝叶斯模型，用 NUTS（HMC 的自适应版本）采样后验
with pm.Model():
    a = pm.Normal("a", mu=0, sigma=10)       # 斜率先验
    b = pm.Normal("b", mu=0, sigma=10)       # 截距先验
    sigma = pm.HalfNormal("sigma", sigma=1)  # 噪声先验
    likelihood = pm.Normal("y", mu=a * x + b, sigma=sigma, observed=y)
    trace = pm.sample(1000, tune=500)        # MCMC 采样
print(pm.summary(trace)[["mean", "sd"]].round(2))  # a≈2.0, b≈1.0

重要采样在强化学习中的简化示意

off-policy 评估：用行为策略 $\mu$ 收集的轨迹评估目标策略 $\pi$ 的期望回报。每条轨迹的权重是其上每步 $\pi/\mu$ 比值的乘积。

import torch
import torch.nn.functional as F

def off_policy_estimate(
    rewards: list[float],           # 实际收到的奖励序列
    log_prob_mu: torch.Tensor,      # 行为策略下各动作的对数概率（旧策略）
    log_prob_pi: torch.Tensor,      # 目标策略下各动作的对数概率（新策略）
    gamma: float = 0.99,            # 折扣因子
) -> float:
    """用重要采样估计目标策略 π 的期望折扣回报。"""
    T = len(rewards)
    discounts = gamma ** torch.arange(T)
    returns = sum(r * d for r, d in zip(rewards, discounts))
    # 重要性权重 = Π π(a|s)/μ(a|s) = exp(Σ [log π - log μ])
    log_ratio = (log_prob_pi - log_prob_mu).sum()
    weight = torch.exp(log_ratio).item()
    return weight * returns

实践中逐轨迹重要采样的方差极大，现代算法（PPO、GRPO）会做截断（clipping）或改用基于 baseline 的优势函数来稳定训练。详见损失函数与正则化相关章节。

实践要点

收敛诊断：MCMC 需要”预热期”（burn-in）丢弃初期未收敛的样本。用 R-hat（Gelman-Rubin 统计量）检查收敛——它比较多条独立链的组间/组内方差，R-hat < 1.01 表示多条链已收敛到同一分布。
样本自相关：MCMC 样本不是独立的（马尔可夫链有记忆），有效样本量（Effective Sample Size, ESS）远小于采样总数。用 ESS 而非原始样本数评估质量；ESS/N 比值越接近 1 说明链效率越高。
维度灾难：拒绝采样在高维下几乎全部被拒（效率趋于零），高维问题必须用 MCMC 或 HMC。即便如此，HMC 在数千维以上也会遇到调参困难。
HMC 需要梯度：HMC 利用对数概率的梯度信息（类似反向传播）加速采样，但要求目标分布可导。PyMC 和 Stan 已内置 NUTS 自动化调参。
变分推断 vs MCMC 的取舍：VI 快但近似（有偏），MCMC 精确但慢（无偏）。大规模数据用 VI，需要精确后验用 MCMC。详见EM算法与变分推断。
重要采样的方差：提议分布 $q$ 选得不好会导致方差爆炸——少数样本权重极大，估计被它们主导。实践中需要监控有效样本量 $\mathrm{ESS} = (\sum w_i)^2 / \sum w_i^2$ 。
步长与稳定性：HMC 的蛙跳步长太大 → 接受率低、能量不守恒；太小 → 每步前进慢、自相关高。NUTS 用对偶平均（dual averaging）在预热期自动把接受率调到约 0.8。

典型应用

贝叶斯推断：PyMC / Stan 用 MCMC（NUTS）估计模型参数的后验分布，广泛应用于统计建模、A/B 测试、临床试验。详见贝叶斯推断基础、模型评估。
强化学习中的策略评估：重要采样用于 off-policy 评估——用旧策略的数据评估新策略的效果，是离线强化学习的核心工具；RLHF 中也用重要采样做方差减小。
扩散模型采样：DDPM（去噪扩散概率模型）的逆向去噪过程本质是逆向随机微分方程（Reverse SDE）的数值求解，属于广义蒙特卡洛方法。每个去噪步都涉及一次条件采样。详见扩散模型。
AlphaGo / AlphaZero 的蒙特卡洛树搜索：用随机模拟（rollout）评估棋局胜率，结合 UCB 选择，在围棋等大状态空间中实现超越人类的决策——这是蒙特卡洛思想在博弈 AI 中的巅峰应用。
物理模拟：分子动力学、量子蒙特卡洛、伊辛模型模拟等用 MCMC 模拟粒子系统的统计行为，材料科学和药物设计的基础工具。
金融衍生品定价：期权定价需要对底层资产价格的随机过程（如几何布朗运动）积分，蒙特卡洛是最常用的数值方法；拟蒙特卡洛已被广泛用于降低定价方差。

典型类库与工具

类库	语言	说明
PyMC	Python	概率编程框架，NUTS 采样器 + 变分推断（ADVI），贝叶斯建模首选
Stan / CmdStanPy	C++/Python	高性能贝叶斯推断引擎，NUTS 采样器，学术统计广泛使用
NumPyro	Python	基于 JAX 的概率编程，GPU 加速 MCMC，Google 团队维护
emcee	Python	仿射不变 MCMC 采样器，天文学和物理领域流行
TensorFlow Probability	Python	Google 概率编程库，支持 HMC 和变分推断
Pyro / NumPyro	Python	基于 PyTorch / JAX 的概率编程，支持 HMC 和 SVI
BlackJAX	Python	基于 JAX 的高性能 MCMC 采样器库，模块化设计，支持 NUTS / HMC / 扩散采样

术语表

术语	英文	解释
蒙特卡洛方法	Monte Carlo Method	用随机采样来近似计算积分、期望等确定量的方法族
拒绝采样	Rejection Sampling	用包络函数采点，落在目标分布下方的接受，其余拒绝
重要采样	Importance Sampling	从提议分布采样，用权重 $p/q$ 修正分布偏差
马尔可夫链	Markov Chain	下一步只依赖当前状态的随机过程
平稳分布	Stationary Distribution	马尔可夫链运行足够久后样本的分布
细致平衡	Detailed Balance	转移概率满足的可逆条件，保证平稳分布为目标分布
Metropolis-Hastings	MH	最通用的 MCMC 算法，提议-接受-拒绝机制
吉布斯采样	Gibbs Sampling	逐维度从条件分布采样的 MCMC 方法，MH 的特例
哈密顿蒙特卡洛	HMC	引入动量变量和梯度信息加速采样的 MCMC 方法
NUTS	No-U-Turn Sampler	HMC 的自适应版本，自动选路径长度，PyMC / Stan 默认
预热期	Burn-in	MCMC 初期未收敛的样本，丢弃不用
R-hat	Gelman-Rubin 统计量	比较多条链的组内/组间方差判断是否收敛
有效样本量	Effective Sample Size (ESS)	衡量自相关后等效的独立样本数
准蒙特卡洛	Quasi-Monte Carlo (QMC)	用低差异序列代替随机数，收敛更快
对偶变量 / 控制变量	Antithetic / Control Variates	利用相关性降低估计方差的经典技巧

2025-2026 年进展

蒙特卡洛方法在 2025-2026 年迎来了一波与生成式 AI 深度融合的新进展，主要沿着三条线展开。

1. 扩散模型采样与 ODE/SDE 数值方法的统一

扩散模型的逆向去噪过程既可以视为随机微分方程（SDE）的蒙特卡洛模拟，也可以视为其对应的确定性常微分方程（ODE）的数值积分。2025 年的研究热点是：

确定性 vs 随机性的权衡：DDPM 的随机版本（每步加噪）有更强的误差纠正能力但方差大；DDIM 的确定性版本方差为零但有偏。新一代采样器（如 EDM2、DPM-Solver++(3M)）通过高阶 ODE 求解器把采样步数压到 10-20 步，同时保持图像质量。
准蒙特卡洛（QMC）进入扩散采样：用 Sobol / Halton 低差异序列替代伪随机数来驱动采样噪声，实验显示在相同步数下 FID 指标可降低 10-20%，且不增加任何计算成本。这把经典统计学的方差减小技术直接接入了生成模型推理。
一致性模型（Consistency Models）：OpenAI、DeepMind 等团队在 2024-2025 年推动的”一致性蒸馏”，把多步扩散采样蒸馏为 1-2 步生成，本质上是把蒙特卡洛模拟的训练信号压缩为单次前向传播。

2. MCMC 与可微编程 / JAX 生态的融合

传统 MCMC（Stan / PyMC）以 Python 前端 + C++ 后端为主，2025 年的趋势是全面转向基于 JAX 的可微、可 JIT 编译、可 GPU 并行的生态：

BlackJAX 与 NumPyro 的成熟：把 HMC / NUTS / NUTSSampler 写成纯函数式 JAX 代码，能直接在 TPU/GPU 上做大规模并行采样。相比传统框架，在贝叶斯神经网络（数千维参数）上的吞吐量可提升 10-100 倍。
可微 MCMC：由于 JAX 的自动微分天然可用，研究者可以把 MCMC 采样步骤嵌入更大的可微计算图（如 MCMC-as-Variational-Inference），用梯度下降联合优化采样器超参数和模型参数。
并行回火（Parallel Tempering）的复兴：在多峰后验（如混合模型、Bayesian 网络结构学习）中，标准 NUTS 会困在局部峰里。2025 年的工作用 JAX 的 vmap 把多条不同”温度”的链并行跑，通过周期性交换样本实现跨峰跳跃，在大规模贝叶斯推断中恢复了可靠性。

3. LLM 时代的重要采样与方差减小

随着 RLHF / RLAIF 成为大模型对齐的标准流程，重要采样重新成为研究焦点：

GRPO（Group Relative Policy Optimization）：DeepSeek 在 2024 年底提出、2025 年被广泛采用的 GRPO 算法，放弃了逐轨迹重要采样的显式权重，改为对一组采样轨迹用组内均值做 baseline（即控制变量思想），大幅降低了策略梯度的方差，同时避免了 PPO 的价值网络开销。这是控制变量技术在 LLM 训练中的大规模成功实践。
自一致性（Self-Consistency）与多样本聚合：在 LLM 推理（数学题、代码生成）中，对同一问题采样多条推理链再投票/校验，本质就是蒙特卡洛估计”最可能正确答案”的期望。详见自一致性解码。
蒙特卡洛树搜索（MCTS）回归 LLM 推理：受 AlphaGo 启发，2025 年 o1 / DeepSeek-R1 等推理模型在训练和测试时重新引入了基于采样的搜索（MCTS / beam search 的混合），用蒙特卡洛模拟估计不同推理路径的”价值”。详见Beam Search。

4. 准蒙特卡洛在科学机器学习中的扩展

在物理信息神经网络（PINN）、算子学习等科学计算场景，高维积分仍是瓶颈。2025 年的工作把 QMC 序列与神经网络训练结合：

用 Sobol 序列替代小批量的随机采样，使 PINN 在 PDE 求解中的损失景观更平滑，收敛更稳定。
在神经辐射场（NeRF）/ 3D Gaussian Splatting 的体渲染积分中，QMC 采样减少了渲染噪声，提升了新视角合成的质量。