与蒙特卡洛直接相关的算法和数据结构详解
与蒙特卡洛直接相关的算法和数据结构详解
1. 蒙特卡洛方法到底在做什么
1.1 一句话定义
蒙特卡洛方法是一类以随机采样为核心的计算方法:当一个量难以直接求出时,生成许多样本,用样本的统计结果近似真实答案。
它不是某一个固定算法,而是一套共同思想。它可以用来:
- 估计积分、概率和期望;
- 从复杂概率分布中采样;
- 追踪随时间变化的隐藏状态;
- 在巨大决策树中寻找较优动作;
- 用多个精度层级降低仿真成本。
1.2 最核心的期望公式
设随机变量 \(X\) 的概率密度为 \(p(x)\),我们关心函数 \(f(X)\) 的平均值:
怎么读
“真实答案 \(\mu\),等于在分布 \(p\) 下,对 \(f(X)\) 求平均。”
符号含义
| 符号 | 含义 |
|---|---|
| \(X\) | 随机变量,例如一次随机生成的状态 |
| \(p(x)\) | \(X\) 在位置 \(x\) 附近出现的概率密度 |
| \(f(x)\) | 我们真正关心的量,例如收益、损失或函数值 |
| \(\mathbb{E}_{p}\) | 按照分布 \(p\) 取平均 |
| \(\mu\) | 想得到的真实期望 |
直觉解释
积分 \(\int f(x)p(x)\,dx\) 可以理解成:
- 每个位置 \(x\) 有一个结果 \(f(x)\);
- 这个位置出现的可能性由 \(p(x)\) 决定;
- 用出现概率作为权重,把所有结果加权平均。
当这个积分难以解析计算时,可以直接模拟“按 \(p\) 出现”的样本。
1.3 样本平均为什么能近似真实平均
从 \(p(x)\) 中独立采样:
然后计算:
怎么读
“把 \(N\) 次模拟得到的函数值相加,再除以 \(N\)。”
上标“帽子”表示它是估计值,不是真实值。下标 \(N\) 表示它由 \(N\) 个样本计算得到。
为什么成立
在期望和方差满足常见有限性条件时,大数定律保证:
这不是说有限样本时一定准确,而是说样本越多,偏离真实值很远的可能性通常越小。
1.4 误差为什么下降得慢
若 \(\operatorname{Var}[f(X)]=\sigma^2<\infty\),样本平均的方差为:
因此标准误差为:
怎么读
“样本数增加 \(N\) 倍,误差只缩小到原来的 \(1/\sqrt{N}\)。”
例如:
- 样本数乘以 \(4\),标准误差约减半;
- 样本数乘以 \(100\),标准误差约缩小到十分之一;
- 想多得到一位十进制精度,往往需要大约 \(100\) 倍样本。
这就是蒙特卡洛方法常常需要方差缩减、重要性采样、QMC 或 MLMC 的原因。
1.5 一套统一的观察框架
几乎所有直接相关的蒙特卡洛算法,都可以用以下五个问题理解:
- 样本是什么? 状态、参数、粒子、路径还是树节点?
- 样本从哪里来? 独立分布、提议分布、状态转移还是搜索策略?
- 样本是否带权? 每个样本贡献相同,还是有不同权重?
- 样本如何汇总? 求均值、加权平均、反向传播还是分层相加?
- 如何判断结果可信? 标准误差、ESS、收敛诊断还是访问次数?
2. 直接相关算法的整体地图
| 分支 | 核心思想 | 典型算法 | 核心数据结构 |
|---|---|---|---|
| 普通蒙特卡洛 | 独立采样后求平均 | Monte Carlo Integration | 样本集合、在线统计 |
| 采样改进 | 改变采样方式以提高效率 | 重要性采样、拒绝采样、分层采样 | 带权样本集合、接受样本缓冲区 |
| MCMC | 构造以目标分布为平稳分布的链 | MH、Gibbs、HMC、NUTS | 马尔可夫链轨迹 |
| SMC | 用带权粒子逐步逼近分布 | 粒子滤波、序贯重要性采样 | 加权粒子集合 |
| MCTS | 用随机模拟估计动作长期价值 | UCT、PUCT | 蒙特卡洛搜索树或图 |
| QMC | 用低差异序列均匀覆盖空间 | Sobol、Halton | 低差异序列生成器 |
| MLMC | 用多层差值替代单一高精度估计 | 多层蒙特卡洛 | 分层样本统计 |
最核心的四个分支可以记成:
3. 普通蒙特卡洛:用样本平均代替真实平均
3.1 蒙特卡洛积分
要计算一维积分:
令 \(X\) 在区间 \([a,b]\) 上均匀分布。均匀分布的密度为 \(1/(b-a)\),于是:
把等式重新排列:
所以蒙特卡洛估计量为:
怎么读
“在区间内随机取点,计算这些点上的函数值平均数,再乘区间长度。”
小例子:计算 \(\int_0^1x^2\,dx\)
真实答案是 \(1/3\)。随机生成 \(N\) 个 \([0,1]\) 内的点:
假设只抽到五个点:\(0.1,0.4,0.6,0.8,0.9\),则:
五个样本很少,所以结果与 \(1/3\) 有明显差异。增加样本后通常会更稳定。
3.2 估计概率
事件 \(A\) 的概率可以写成指示函数的期望:
其中:
因此概率估计就是“命中次数除以总次数”:
这解释了估计圆周率时为什么可以统计落在圆内的点的比例。
3.3 置信区间
当样本量足够大、方差有限且样本近似独立时,可用中心极限定理构造近似置信区间:
符号含义
| 符号 | 含义 |
|---|---|
| \(\hat\mu_N\) | 样本均值 |
| \(s\) | 样本标准差 |
| \(N\) | 样本数 |
| \(z_{1-\alpha/2}\) | 标准正态分布分位数;95% 区间时约为 \(1.96\) |
注意
这个公式不适合被机械使用。重尾分布、样本高度相关、样本量过小或方差不存在时,正态近似可能很差。MCMC 样本还必须考虑自相关,不能直接把总迭代次数当成独立样本数。
3.4 直接蒙特卡洛伪代码
输入:样本数 N、采样分布 p、目标函数 f
输出:期望估计与标准误差
初始化在线统计量 statistics
for i = 1 to N:
x = Sample(p)
y = f(x)
statistics.update(y)
estimate = statistics.mean
standard_error = statistics.standard_deviation / sqrt(N)
return estimate, standard_error
3.5 核心数据结构:样本集合与在线统计
SampleSet
├── samples[] 可选:原始样本
├── values[] 可选:f(samples[i])
├── sample_count
├── running_mean
├── running_variance
└── random_seed
若只需要均值和方差,没有必要保存全部样本。可以使用 Welford 在线统计结构:
OnlineStatistics
├── count
├── mean
└── M2 所有样本相对当前均值的平方偏差累计量
加入新值 \(x_n\) 时:
最终样本方差为:
直觉解释
Welford 算法不需要保存历史数据。每来一个新样本,就小幅修正均值,并记录这个新样本对总离散程度的贡献。它比“先求平方和再相减”的写法更稳定。
3.6 方差缩减
普通蒙特卡洛的关键瓶颈不是偏差,而往往是方差。常见方法包括:
- 对偶变量:成对使用负相关样本,例如 \(U\) 与 \(1-U\);
- 控制变量:利用期望已知、且与目标量相关的变量修正估计;
- 分层采样:把空间切成多个区域,避免样本偶然挤在某一处;
- 条件蒙特卡洛:先解析消掉一部分随机性;
- 拉丁超立方采样:让每个坐标方向的覆盖更均匀。
这些方法的共同目标不是“让样本更多”,而是“让每个样本提供更多有效信息”。
4. 重要性采样与拒绝采样
4.1 重要性采样:换一个更会抽重点的分布
目标仍然是:
如果直接从 \(p\) 采样很困难,或真正重要的区域在 \(p\) 下很少出现,可以引入更容易采样的分布 \(q\):
怎么读
“虽然样本来自 \(q\),但用 \(p/q\) 进行补偿,仍然可以计算在 \(p\) 下的平均。”
权重定义为:
估计量为:
为什么要乘 \(p/q\)
假设某一区域在 \(q\) 中被抽得比在 \(p\) 中更频繁,那么该区域的每个样本就应该少算一点;反之则多算一点。
- \(p(x)>q(x)\):该点在 \(q\) 下抽得不够多,权重 \(p/q>1\);
- \(p(x)<q(x)\):该点在 \(q\) 下抽得太多,权重 \(p/q<1\)。
可以把它理解为“换了抽样账本之后的会计修正”。
必要条件
只要某处 \(f(x)p(x)\) 可能非零,\(q(x)\) 就不能为零。否则那个区域永远抽不到,再大的权重也无法修正。
形式上要求:
4.2 自归一化重要性采样
有时只知道目标密度的未归一化形式:
其中未知常数 \(Z\) 无法计算。此时使用未归一化权重:
再归一化:
估计量为:
直觉解释
未知常数 \(Z\) 同时出现在所有权重中。归一化时,分子分母中的 \(Z\) 会相互抵消,因此不需要知道它。
严谨性提醒
自归一化估计在有限样本下通常有偏,但在常见正则条件下是一致的,即样本数增加时趋近真实值。
4.3 有效样本量 ESS
若一个样本几乎占据全部权重,其他样本就几乎没有贡献。常用近似指标是:
怎么看这个公式
- 若所有权重相等,\(\bar w_i=1/N\):
这表示 \(N\) 个样本都充分有效。
- 若一个权重为 \(1\),其余为 \(0\):
这表示虽然存了 \(N\) 个样本,但实际上只有一个样本在决定结果。
因此:
ESS 是很有用的退化诊断,但它是经验指标,不应被误解成严格等价于真正的独立样本数。
4.4 核心数据结构:带权样本集合
WeightedSampleSet
├── samples[N][D]
├── log_weights[N]
├── normalized_weights[N]
├── cumulative_weights[N]
├── effective_sample_size
└── weighted_statistics
实践中优先保存 log_weights,因为大量小概率相乘很容易下溢为零。
4.5 拒绝采样:先多提候选,再按比例接受
选择容易采样的提议密度 \(q(x)\) 和常数 \(M\),满足:
每次先采样:
再按以下条件接受:
几何直觉
曲线 \(Mq(x)\) 像一张完全盖住 \(p(x)\) 的“外罩”。先在外罩下抽一个点,只有落在目标曲线 \(p(x)\) 下方时才接受。
若 \(p\) 和 \(q\) 都已归一化,则平均接受率为:
所以 \(M\) 越接近 \(1\) 越好;\(M\) 很大意味着大量候选会被浪费。
输入:目标密度 p、提议密度 q、包络常数 M
重复:
x ~ q
u ~ Uniform(0, 1)
if u <= p(x) / (M q(x)):
输出 x
主要限制
在高维空间中,找到紧密包络通常很困难。即使每一维只多出一点空隙,整体体积差也可能随维度迅速放大,导致接受率极低。
5. MCMC:用马尔可夫链从复杂分布采样
5.1 为什么需要 MCMC
贝叶斯推断中常遇到:
这里 \(\pi(x)\) 是目标后验分布。符号 \(\propto\) 表示“成比例”:我们能计算分子形状,但归一化常数可能很难求。
普通蒙特卡洛要求能直接从 \(\pi\) 采样。MCMC 的办法是:构造一条随机移动的链,让它长期停留频率恰好符合 \(\pi\)。
5.2 马尔可夫链是什么
马尔可夫性质写成:
怎么读
“知道现在的状态后,预测下一步不再需要更早的历史。”
这不代表历史完全没有影响,而是历史影响已经被当前状态概括。
5.3 Metropolis–Hastings 接受率
当前状态为 \(x\),先从提议分布生成候选:
接受概率为:
这是 MCMC 中最容易觉得晦涩的公式,可以拆成两部分:
第一部分:目标分布比值
它比较候选状态和当前状态谁更符合目标分布。
- 若候选更可能,比例大于 \(1\),通常直接接受;
- 若候选更不可能,仍可能以一定概率接受,避免链被局部高点困住。
第二部分:提议方向修正
如果从 \(x\) 到 \(x'\) 很容易,而从 \(x'\) 回到 \(x\) 很难,那么单靠目标概率比会产生方向偏差。这一项负责修正提议机制的不对称。
对称提议的简化
若:
提议修正相互抵消,得到 Metropolis 接受率:
一个数值例子
若候选状态的目标密度只有当前状态的一半,而且提议对称:
这意味着候选虽然更差,但仍有一半概率被接受。允许“偶尔下坡”正是 MCMC 能穿过低概率区域的原因之一。
5.4 为什么这个接受率能工作:详细平衡
一个常见的充分条件是详细平衡:
其中 \(K(x,x')\) 是链从 \(x\) 转移到 \(x'\) 的概率或密度。
直觉解释
在长期稳定状态下,从 \(x\) 流向 \(x'\) 的概率质量,与从 \(x'\) 流回 \(x\) 的概率质量相等。于是总体分布不会再发生变化,\(\pi\) 成为平稳分布。
详细平衡是常用的充分条件,但不是平稳性的必要条件。并非所有正确的 MCMC 方法都必须满足详细平衡。
5.5 Metropolis–Hastings 伪代码
输入:初始状态 x0、目标未归一化密度 π、提议分布 q、迭代次数 T
x = x0
for t = 1 to T:
x_candidate ~ q(. | x)
log_ratio = log π(x_candidate) - log π(x)
+ log q(x | x_candidate) - log q(x_candidate | x)
log_u = log Uniform(0, 1)
if log_u <= min(0, log_ratio):
x = x_candidate
accepted[t] = true
else:
accepted[t] = false
trace[t] = x
return trace
使用对数比值可以避免直接计算极小概率。
5.6 Gibbs 采样
若 \(x=(x_1,\ldots,x_d)\),Gibbs 采样依次更新每个分量:
怎么读
更新第 \(j\) 个变量时:
- 已经在本轮更新过的变量,使用新值;
- 尚未更新的变量,使用上一轮的旧值;
- 只随机生成当前这个变量。
直觉解释
直接在高维联合分布里抽样可能很难,但固定其他变量后,每一个条件分布可能很简单。Gibbs 采样通过“逐坐标刷新”完成整体采样。
5.7 HMC:把随机游走变成有方向的运动
HMC 给参数 \(x\) 引入辅助动量 \(r\),定义哈密顿量:
其中:
直觉解释
- \(U(x)\) 是“势能”:目标概率越高,\(-\log\pi(x)\) 越低;
- \(K(r)\) 是“动能”:控制移动方向和速度;
- 梯度帮助轨迹沿着高概率区域远距离移动,而不是像随机游走一样小步乱晃。
数值积分会产生能量误差,所以 HMC 末尾仍使用 Metropolis 接受步骤进行修正。
NUTS 在 HMC 基础上自动选择轨迹长度,避免走得太短或已经折返还继续计算。
5.8 核心数据结构:MCMC 轨迹
MCMCTrace
├── samples[T][D]
├── log_probabilities[T]
├── accepted[T]
├── proposal_or_hmc_states
├── warmup_length
├── chain_id
└── diagnostics
├── acceptance_rate
├── autocorrelation
├── effective_sample_size
└── convergence_statistics
概念上的链是:
但程序中通常使用连续数组或张量,而不是链表。这里的“链”指概率依赖关系,不是编程语言中的链式存储。
5.9 自相关与 MCMC 有效样本量
MCMC 相邻样本相关,因此 \(T\) 次迭代不等于 \(T\) 个独立样本。一个常见关系是:
其中 \(\rho_k\) 是间隔 \(k\) 的自相关。
怎么读
分母是“相关性惩罚”。若自相关都接近零,ESS 接近 \(T\);若正相关很强,ESS 会远小于 \(T\)。
实际软件会使用截断、谱估计或多链方法稳定估计 ESS,因此不应直接把无限和式机械实现。
5.10 MCMC 的关键风险
- 预热不足:链仍受初始状态强烈影响;
- 混合缓慢:链在参数空间移动太慢;
- 多峰困境:链长期停留在一个模态;
- 步长不合适:过小则移动慢,过大则接受率低;
- 只看轨迹图不够:还应结合多链诊断、ESS 和模型检查。
6. SMC 与粒子滤波:用一群带权粒子追踪分布
6.1 状态空间模型
动态系统通常包含看不见的隐藏状态 \(x_t\) 和可观察数据 \(y_t\):
怎么读
- 第一式:下一时刻状态由上一时刻状态随机演化;
- 第二式:观测由当前隐藏状态随机生成。
目标是在看到 \(y_1,\ldots,y_t\) 后估计:
这里 \(y_{1:t}\) 表示从第 \(1\) 时刻到第 \(t\) 时刻的全部观测。
6.2 用粒子表示一个分布
粒子滤波使用带权点集合近似后验分布:
这个公式怎么理解
它看起来抽象,其实表达的是:
“我们不再保存一条连续的概率密度曲线,而是用 \(N\) 个带权代表点来近似它。”
| 符号 | 含义 |
|---|---|
| \(x_t^{(i)}\) | 第 \(i\) 个粒子在时刻 \(t\) 的状态 |
| \(w_t^{(i)}\) | 这个粒子的可信权重,通常满足权重和为 \(1\) |
| \(\delta(\cdot)\) | 狄拉克点质量,表示概率集中在该粒子位置 |
不要把 \(\delta\) 理解成普通函数。这里它只是一种紧凑写法,用于表达“离散加权点集所代表的经验分布”。
6.3 Bootstrap 粒子滤波的权重更新
先传播粒子:
然后用当前观测计算未归一化权重:
归一化:
直觉解释
- 状态转移模型负责“预测下一步可能在哪里”;
- 观测似然负责“根据新证据给预测打分”;
- 归一化后,所有粒子的权重重新组成一个概率分布。
6.4 为什么需要重采样
多次更新后,常出现少数粒子权重很大,大多数权重接近零。此时计算资源被浪费在几乎没有贡献的粒子上。
使用:
判断退化程度。当 ESS 低于阈值时进行重采样。
重采样的直觉
- 高权重粒子被复制多次;
- 低权重粒子大概率被淘汰;
- 重采样后通常把权重重置为 \(1/N\)。
重采样不会凭空创造新状态,它只重新分配计算资源。因此重采样过于频繁也可能损失粒子多样性。
6.5 粒子滤波伪代码
初始化:x[i] ~ initial_distribution,w[i] = 1/N
for t = 1 to T:
for i = 1 to N:
x_new[i] ~ Transition(. | x[i])
log_weight[i] = log w[i] + log Likelihood(y[t] | x_new[i])
weights = NormalizeLogWeights(log_weight)
ESS = 1 / sum(weights[i]^2)
if ESS < threshold:
ancestor_indices = Resample(weights)
x = x_new[ancestor_indices]
w = [1/N, ..., 1/N]
else:
x = x_new
w = weights
estimate[t] = WeightedMean(x, w)
6.6 核心数据结构:加权粒子集合
ParticleSet
├── states[N][D]
├── log_weights[N]
├── normalized_weights[N]
├── ancestor_indices[N]
├── particle_ids[N]
├── effective_sample_size
├── time_step
└── resampling_statistics
祖先索引非常重要。它记录重采样后的每个粒子来自上一步的哪个粒子,可用于路径回溯、平滑和谱系分析。
高性能实现通常采用“结构分离数组”:状态、权重、祖先索引分别连续存放,而不是创建 \(N\) 个包含许多字段的小对象。这样更利于缓存和 GPU 批处理。
6.7 常见重采样方法
| 方法 | 核心做法 | 特点 |
|---|---|---|
| 多项式重采样 | 每次独立按权重抽祖先 | 简单,但随机方差较大 |
| 系统重采样 | 用一个随机偏移生成等间隔抽样点 | 快,常用,方差较低 |
| 分层重采样 | 每个等分区间各抽一个随机点 | 降低随机波动 |
| 残差重采样 | 先确定性复制整数部分,再抽剩余部分 | 减少不必要随机性 |
7. MCTS:用随机模拟指导树搜索
7.1 MCTS 解决什么问题
当动作分支很多、搜索深度很大、无法穷举全部决策时,MCTS 通过重复模拟,把计算预算集中到更有希望的分支。
一次迭代分为:
- 选择:沿当前树选择值得继续探索的节点;
- 扩展:创建一个尚未尝试的子节点;
- 模拟:从新节点出发,用随机或启发式策略获得结果;
- 反向传播:把结果沿路径更新回去。
7.2 节点统计量
对节点或动作边,常保存:
| 符号 | 含义 |
|---|---|
| \(N_j\) | 被访问的次数 |
| \(W_j\) | 所有模拟收益之和 |
| \(Q_j\) | 平均模拟收益 |
直觉解释
\(Q_j\) 是“目前看起来有多好”,\(N_j\) 是“这个判断建立在多少次试验上”。一个高收益但只试过一次的节点,可信度通常不如一个试过很多次且收益稳定的节点。
7.3 UCT:利用与探索的平衡
UCT 评分为:
逐项解释
利用项
优先选择目前平均收益高的节点。
探索项
- 父节点被访问得越多,仍未充分尝试的子节点越值得补试;
- 子节点访问次数 \(N_j\) 越小,探索奖励越大;
- \(C\) 控制探索强度。
为什么要取对数和平方根
这来自多臂老虎 JI的上置信界思想:探索奖励应随总试验次数缓慢增长,并随某个选项自身试验次数增加而衰减。
未访问节点怎么办
当 \(N_j=0\) 时公式会除以零。实现中通常把未访问节点的分数视为正无穷,或在进入 UCT 比较前优先访问每个动作至少一次。
7.4 PUCT:加入策略先验
PUCT 常写为:
符号含义
| 符号 | 含义 |
|---|---|
| \(Q(s,a)\) | 在状态 \(s\) 选择动作 \(a\) 的平均价值 |
| \(P(s,a)\) | 策略模型给出的先验概率 |
| \(N(s)\) | 状态 \(s\) 的总访问次数 |
| \(N(s,a)\) | 动作 \(a\) 的访问次数 |
| \(c_{\mathrm{puct}}\) | 先验探索强度 |
直觉解释
UCT 主要根据“试出来的结果”探索;PUCT 还会参考“策略模型预先认为哪些动作值得尝试”。但随着某动作访问次数增加,先验奖励逐渐衰减,真实模拟统计会越来越重要。
7.5 核心数据结构:蒙特卡洛搜索树
MCTSNode
├── state_or_state_reference
├── parent
├── children 动作到子节点的映射
├── unexpanded_actions
├── visit_count N
├── total_value W
├── mean_value Q = W / N
├── prior_probability P,可选
├── terminal_flag
└── virtual_loss 并行搜索时可选
更严谨的设计常把动作统计放在边上:
MCTSEdge
├── action
├── child_node
├── visit_count
├── total_value
├── mean_value
└── prior_probability
因为 \(Q(s,a)\)、\(N(s,a)\) 本质上属于“在状态 \(s\) 选择动作 \(a\)”这条边,而不只是子状态本身。
7.6 置换表:树可能变成图
同一个状态可能由不同动作序列到达。可使用:
TranspositionTable
state_hash -> node_or_shared_statistics
合并重复状态。这样可以减少重复计算,但底层结构会从严格树变成有向图。实现时要处理循环、引用计数和不同父路径的价值视角。
7.7 MCTS 伪代码
root = CreateNode(initial_state)
repeat simulation_budget times:
node = root
path = []
while node is fully expanded and not terminal:
edge = SelectByUCTorPUCT(node)
path.append(edge)
node = edge.child
if node is not terminal:
edge = ExpandOneAction(node)
path.append(edge)
node = edge.child
reward = SimulateOrEvaluate(node.state)
for edge in reverse(path):
edge.visit_count += 1
edge.total_value += reward_from_edge_player_view
严谨性提醒:收益视角
两人零和博弈中,向父层回传时可能需要改变收益符号。若不明确“价值属于当前玩家、根玩家还是固定玩家”,很容易写出表面运行正常但统计方向错误的实现。
8. QMC:用更均匀的点代替纯随机点
8.1 为什么纯随机点可能浪费样本
纯随机点可能偶然聚成一团,也可能留下空洞。QMC 使用低差异序列,使前 \(N\) 个点尽量均匀覆盖积分区域。
常见序列:
- Sobol 序列;
- Halton 序列;
- Faure 序列;
- Niederreiter 序列。
8.2 QMC 估计量
在单位超立方体 \([0,1]^d\) 上:
QMC 使用确定性点列 \(x_1,\ldots,x_N\):
公式和普通蒙特卡洛的样本平均几乎相同,区别在于点的产生方式:
- MC:点是独立伪随机样本;
- QMC:点是按低差异规则构造的有序序列。
8.3 “低差异”是什么意思
可以把单位区域切成许多小盒子。理想情况下,每个盒子里的点数应该接近“盒子体积 × 总点数”。差异度衡量实际点数与理想点数之间的最大偏差。
一个经典理论界是 Koksma–Hlawka 不等式:
不必被符号吓到
它只表达一件事:
积分误差不超过“函数有多难积分”乘以“点分布有多不均匀”。
- \(V_{\mathrm{HK}}(f)\):函数的某种总变差,越大表示函数越不平滑;
- \(D_N^*\):星差异度,越小表示点覆盖越均匀。
这个界在高维时可能不够紧,且对函数有条件要求,因此不能简单宣称 QMC 总是优于 MC。
8.4 随机化 QMC
纯 QMC 是确定性的,难以直接用重复实验估计误差。随机化 QMC 对低差异序列进行扰乱或随机平移,使每次运行具有随机性,同时尽量保留均匀覆盖。
这样可以多次独立扰乱,使用不同运行结果之间的样本方差估计误差。
8.5 核心数据结构:低差异序列生成器
LowDiscrepancySequence
├── dimension
├── index
├── direction_numbers
├── base_parameters
├── scrambling_state
└── current_point[D]
它不是简单的随机点数组。第 \(i\) 个点通常由索引、维度参数和方向数增量生成,序列顺序本身具有意义。
9. MLMC:把高精度估计拆成多层差值
9.1 为什么需要多层
很多仿真有多个精度层级:
- 粗网格:便宜,但偏差大;
- 细网格:准确,但一次模拟很贵。
若只在最高精度层做普通蒙特卡洛,成本可能非常高。
9.2 望远镜分解
设 \(P_0,P_1,\ldots,P_L\) 表示从粗到细的多个近似。恒等式:
中间项会相互抵消。例如前三层:
对两边取期望:
直觉解释
不直接反复估计昂贵的 \(P_L\),而是:
- 大量估计便宜的粗层 \(P_0\);
- 分别估计相邻精度之间的小修正;
- 把所有修正加起来。
若相邻层使用耦合随机数,使 \(P_\ell\) 与 \(P_{\ell-1}\) 高度相关,那么差值 \(P_\ell-P_{\ell-1}\) 的方差通常很小,因此高精度层只需要少量样本。
9.3 MLMC 估计量
令:
每层用 \(N_\ell\) 个样本估计:
总估计为:
关键设计原则
- 便宜且方差大的低层:多采样;
- 昂贵且差值方差小的高层:少采样;
- 同一差值样本内,粗细模拟应尽量共享随机性。
9.4 核心数据结构:分层统计
MLMCState
├── levels[0...L]
│ ├── level_index
│ ├── sample_count
│ ├── sum_difference
│ ├── sum_squared_difference
│ ├── estimated_mean
│ ├── estimated_variance
│ └── cost_per_sample
├── target_mean_square_error
└── total_estimate
MLMC 通常不保存全部原始路径,只保存每层的累计统计、成本和样本数。
10. 直接相关的核心数据结构
严格来说,没有一种像“栈”或“红黑树”那样统一定义的“蒙特卡洛数据结构”。更准确的说法是:不同蒙特卡洛算法拥有与其采样、权重更新和统计汇总直接绑定的状态结构。
10.1 随机样本集合 SampleSet
适用:普通 MC、方差缩减、部分 QMC。
SampleSet
├── samples
├── evaluated_values
├── count
├── mean
├── variance
└── metadata
关键操作:
- 添加样本;
- 更新均值和方差;
- 合并多个批次;
- 计算置信区间;
- 分块写盘。
10.2 带权样本集合 WeightedSampleSet
适用:重要性采样、一般加权估计。
WeightedSampleSet
├── samples
├── log_weights
├── normalized_weights
├── cumulative_weights
├── effective_sample_size
└── weighted_moments
加权均值:
10.3 马尔可夫链轨迹 MCMCTrace
适用:MH、Gibbs、HMC、NUTS。
MCMCTrace
├── ordered_samples
├── acceptance_flags
├── log_probabilities
├── sampler_states
├── warmup_boundary
└── diagnostics
核心特点:
- 样本有时间顺序;
- 相邻样本相关;
- 常同时运行多条链;
- 支持丢弃预热段、切片读取和诊断统计。
10.4 加权粒子集合 ParticleSet
适用:SMC、粒子滤波。
ParticleSet
├── states
├── log_weights
├── normalized_weights
├── ancestor_indices
├── time_step
└── ESS
核心特点:
- 整个集合表示一个经验分布;
- 权重随新观测变化;
- 重采样会产生复制和淘汰;
- 祖先索引用于重建历史路径。
10.5 蒙特卡洛搜索树 MonteCarloSearchTree
适用:MCTS、UCT、PUCT。
MonteCarloSearchTree
├── root
├── node_pool
├── edge_pool
├── state_index_or_transposition_table
└── search_statistics
节点和边中的核心字段:
- 状态或状态引用;
- 动作;
- 父子关系;
- 访问次数;
- 累计收益;
- 平均价值;
- 策略先验;
- 终局标记。
10.6 低差异序列状态 LowDiscrepancySequence
适用:QMC、随机化 QMC。
LowDiscrepancySequence
├── dimension
├── current_index
├── direction_numbers_or_bases
├── scrambling_state
└── current_point
它保存的是“如何生成下一个均匀覆盖点”的状态,而不一定保存全部历史点。
10.7 分层样本统计 MLMCState
适用:MLMC。
MLMCState
├── level_statistics[]
├── per_level_cost
├── per_level_variance
├── allocated_sample_counts
└── total_estimate
它的重点是按层管理样本预算,而不是组织单个样本对象。
11. 加权抽样的辅助数据结构
重要性采样和粒子重采样经常需要按离散权重:
随机选择索引 \(i\)。
11.1 前缀和数组
构造:
生成 \(u\in[0,C_N)\),找到第一个满足 \(C_i>u\) 的索引。
CumulativeWeights
├── prefix_sums[N]
└── total_weight
- 构建:\(O(N)\);
- 单次二分采样:\(O(\log N)\);
- 修改一个权重后重新构建:\(O(N)\)。
11.2 Alias Table
把每个离散概率拆成“主选项 + 备用选项”,使采样只需:
- 均匀选择一个桶;
- 再做一次伯努利判断。
AliasTable
├── probability[N]
└── alias[N]
- 构建:\(O(N)\);
- 单次采样:\(O(1)\);
- 权重频繁变化时重建成本较高。
适合固定分布下进行大量重复抽样。
11.3 Fenwick 树与线段树
树节点保存某个范围内的权重和。随机采样时,从根开始根据左子树权重决定向左还是向右。
- 构建:\(O(N)\);
- 单点更新:\(O(\log N)\);
- 单次加权采样:\(O(\log N)\)。
适合权重持续动态变化的场景。
11.4 复杂度比较
| 结构 | 构建 | 单次采样 | 单点更新 | 适用场景 |
|---|---|---|---|---|
| 线性扫描 | \(O(1)\) | \(O(N)\) | \(O(1)\) | 样本很少 |
| 前缀和数组 | \(O(N)\) | \(O(\log N)\) | \(O(N)\) | 权重基本固定 |
| Alias Table | \(O(N)\) | \(O(1)\) | 通常需重建 | 固定权重、大量采样 |
| Fenwick 树 | \(O(N)\) | \(O(\log N)\) | \(O(\log N)\) | 权重频繁变化 |
| 线段树 | \(O(N)\) | \(O(\log N)\) | \(O(\log N)\) | 动态采样及范围统计 |
这些是蒙特卡洛算法的辅助结构,而不是独立的蒙特卡洛方法。
12. 算法与数据结构对应关系
| 算法 | 样本是否独立 | 是否带权 | 核心数据结构 | 最重要的诊断量 |
|---|---|---|---|---|
| 普通 MC | 通常独立 | 通常等权 | 样本集合、在线统计 | 标准误差 |
| 重要性采样 | 通常独立 | 是 | 带权样本集合 | 权重方差、ESS |
| 拒绝采样 | 接受后通常独立 | 否 | 接受样本缓冲区 | 接受率 |
| MH / Gibbs | 否 | 最终通常等权 | MCMC 轨迹 | 自相关、ESS、收敛诊断 |
| HMC / NUTS | 否 | 最终通常等权 | 扩展 MCMC 状态与轨迹 | 发散、ESS、能量诊断 |
| 粒子滤波 / SMC | 条件相关 | 是 | 加权粒子集合 | ESS、粒子多样性 |
| MCTS / UCT / PUCT | 不是概率分布样本的传统含义 | 用访问统计隐式加权 | 搜索树或图 | 访问次数、价值稳定性 |
| QMC | 确定性或随机化序列 | 等权 | 低差异序列生成器 | 多次扰乱间误差 |
| MLMC | 层内通常独立,粗细样本耦合 | 按层组合 | 分层样本统计 | 每层方差与成本 |
13. 复杂度、数值稳定性与工程实现
13.1 普通蒙特卡洛
设单次采样与函数评价成本为 \(C_f\),样本数为 \(N\):
- 时间:\(O(NC_f)\);
- 保存所有样本:空间 \(O(ND)\);
- 只保存在线统计:额外空间可降为 \(O(1)\)。
独立样本天然适合多线程、GPU 和分布式并行。
13.2 MCMC
设迭代数为 \(T\),每次更新成本为 \(C_q\):
- 时间:\(O(TC_q)\);
- 保存 \(D\) 维完整轨迹:空间 \(O(TD)\);
- 单链内部存在顺序依赖,但多链可以并行。
不要为了减少文件大小而盲目“稀疏抽样”。稀疏抽样一般不会增加统计信息,只会丢弃样本;只有存储或后处理成本成为主要瓶颈时才可能合理。
13.3 粒子滤波
时间步为 \(T\),粒子数为 \(N\),状态维度为 \(D\):
- 时间通常为 \(O(TN)\) 乘以单次传播和似然成本;
- 只保留当前粒子:空间 \(O(ND)\);
- 保存全历史:可能达到 \(O(TND)\)。
可使用祖先索引压缩历史路径,但若频繁回溯仍需注意内存和随机访问成本。
13.4 MCTS
模拟次数为 \(S\),平均深度为 \(H\):
- 主循环大致为 \(O(SH)\);
- 实际成本还取决于状态复制、合法动作生成和 rollout;
- 空间与已扩展节点、边及状态表示大小有关。
常见优化:节点池、状态增量更新、置换表、延迟展开、虚拟损失、批量神经网络推理。
13.5 对数域权重
直接计算:
可能因概率极小而下溢。改为:
归一化时使用 Log-Sum-Exp。设 \(a_i=\log w_i\),取:
则:
为什么有效
\(a_i-m\le0\),所以指数值不会爆炸;最大的那一项变为 \(e^0=1\),也不容易全部下溢为零。
归一化对数权重可写为:
13.6 随机数管理
可复现实现至少应记录:
- 随机数生成器算法;
- 主种子及子流生成方式;
- 线程、进程、GPU 设备的独立随机流;
- 算法版本与参数;
- 样本数、批大小和并行配置。
并行环境中不能简单让所有线程使用同一个种子,否则可能生成重复或强相关随机流。
13.7 内存布局
对大量粒子、链状态和树节点:
- 数值计算密集型数据优先使用连续数组;
- 热字段与冷字段可分离;
- 尽量避免每个样本单独动态分配对象;
- 大轨迹采用分块写盘或内存映射;
- MCTS 用节点池和整数索引常比大量指针更紧凑。
14. 容易混淆的概念与边界
14.1 概率数据结构不等于蒙特卡洛核心数据结构
Bloom Filter、Count-Min Sketch、HyperLogLog、MinHash 等允许概率误差,常被称为概率数据结构,某些理论语境也会称其为 Monte Carlo 型数据结构。
但它们不是 MC、MCMC、SMC 或 MCTS 计算流程的核心状态结构。本文所说的“直接相关数据结构”主要是:
- 样本集合;
- 带权样本集合;
- MCMC 轨迹;
- 粒子集合;
- 蒙特卡洛搜索树;
- 低差异序列状态;
- 分层样本统计。
14.2 随机化算法不一定是蒙特卡洛算法
随机化算法常分为:
- Monte Carlo 型:运行时间通常有界,但答案可能有小概率错误或估计误差;
- Las Vegas 型:答案保证正确,但运行时间带随机性。
Treap、Skip List、随机化快速排序虽然使用随机数,但不是本文讨论的蒙特卡洛采样体系的核心算法。
14.3 QMC 的名称边界
纯 QMC 使用确定性低差异序列,严格说并不依赖随机采样,因此有时被视为与蒙特卡洛并列的方法。随机化 QMC 则明确重新引入随机性。它们仍因目标、估计形式和应用领域高度一致而通常放在蒙特卡洛体系中讨论。
14.4 模拟退火与交叉熵方法
它们与蒙特卡洛密切相关,但核心目标是优化而非估计分布:
- 模拟退火通过温度控制接受较差解的概率;
- 交叉熵方法迭代更新采样分布,使样本集中到高质量区域。
它们可视为蒙特卡洛思想在随机优化中的扩展,而不是本文最核心的四类结构。
14.5 MCMC 的“链”不是链表
MCMC 中的链描述“下一状态依赖当前状态”的概率关系。底层最常用数组、张量或分块文件,而不是链表。
14.6 MCTS 的“树”不一定严格是树
使用置换表合并重复状态后,多个父节点可能指向同一状态,结构会变成有向图。因此应区分算法名称与实际内存拓扑。
15. 选型指南
| 问题特征 | 推荐方法 | 直接核心结构 | 首要风险 |
|---|---|---|---|
| 可直接独立采样,只需估计期望 | 普通 MC | 在线统计或样本集合 | 方差大、收敛慢 |
| 稀有事件或重点区域难抽到 | 重要性采样 | 带权样本集合 | 权重退化、支撑集不覆盖 |
| 有容易构造的包络分布,维度不高 | 拒绝采样 | 接受样本缓冲区 | 接受率过低 |
| 只知道未归一化目标密度 | MH | MCMC 轨迹 | 混合慢、多峰 |
| 条件分布容易采样 | Gibbs | 多维链轨迹 | 变量强相关时移动慢 |
| 高维连续变量且可求梯度 | HMC / NUTS | 位置、动量与轨迹 | 几何病态、数值发散 |
| 状态随时间变化并持续接收观测 | 粒子滤波 / SMC | 加权粒子集合 | 权重退化、粒子贫化 |
| 巨大决策树无法穷举 | MCTS / UCT / PUCT | 搜索树或图 | 模拟偏差、价值视角错误 |
| 中低有效维积分,函数较平滑 | QMC / 随机化 QMC | 低差异序列生成器 | 高有效维时优势下降 |
| 有粗细多个仿真层级 | MLMC | 分层样本统计 | 粗细耦合不佳、层分配错误 |
| 固定离散权重,需要海量抽样 | Alias Method | Alias Table | 权重更新昂贵 |
| 动态权重,需要持续抽样 | Fenwick 树或线段树 | 动态权重树 | 常数开销较高 |
16. 总结
16.1 算法层面
直接与蒙特卡洛相关的主要算法是:
它们的本质区别可以压缩成一句话:
- 普通 MC:生成独立样本后求平均;
- 重要性采样:从别处分布抽样,再用权重修正;
- 拒绝采样:从包络分布提候选,只保留合格样本;
- MCMC:通过有依赖的随机链间接获得目标分布样本;
- SMC:让一群带权粒子随时间传播、加权和重采样;
- MCTS:把随机模拟结果累计在搜索树节点或边上;
- QMC:用覆盖更均匀的序列计算样本平均;
- MLMC:把昂贵高精度答案拆成便宜粗层和小修正。
16.2 数据结构层面
直接绑定这些算法的核心结构是:
最重要的对应关系:
| 分支 | 主要对象 | 必须保存的核心信息 |
|---|---|---|
| 普通 MC | 样本或在线统计 | 数量、均值、方差 |
| 重要性采样 | 带权样本 | 对数权重、归一化权重、ESS |
| MCMC | 有序轨迹 | 状态、接受标记、概率与诊断 |
| SMC | 粒子群 | 状态、权重、祖先、时间步 |
| MCTS | 搜索节点和边 | 访问次数、累计价值、先验 |
| QMC | 序列生成器 | 维度、索引、方向数、扰乱状态 |
| MLMC | 多层统计 | 各层样本数、差值方差、成本 |
最终可以把“蒙特卡洛数据结构”理解为:
为了生成、保存、更新、重采样和汇总蒙特卡洛样本而设计的算法状态结构。
它不是一种单独的数据结构类别,而是一组与不同蒙特卡洛算法紧密绑定的表示方式。
未经作者同意请勿转载
本文来自博客园作者:aixueforever,原文链接:https://www.cnblogs.com/aslanvon/p/22188388

浙公网安备 33010602011771号