与蒙特卡洛直接相关的算法和数据结构详解

与蒙特卡洛直接相关的算法和数据结构详解

1. 蒙特卡洛方法到底在做什么

1.1 一句话定义

蒙特卡洛方法是一类以随机采样为核心的计算方法:当一个量难以直接求出时,生成许多样本,用样本的统计结果近似真实答案。

它不是某一个固定算法,而是一套共同思想。它可以用来:

  • 估计积分、概率和期望;
  • 从复杂概率分布中采样;
  • 追踪随时间变化的隐藏状态;
  • 在巨大决策树中寻找较优动作;
  • 用多个精度层级降低仿真成本。

1.2 最核心的期望公式

设随机变量 \(X\) 的概率密度为 \(p(x)\),我们关心函数 \(f(X)\) 的平均值:

\[\mu=\mathbb{E}_{p}[f(X)] =\int f(x)p(x)\,dx. \]

怎么读

“真实答案 \(\mu\),等于在分布 \(p\) 下,对 \(f(X)\) 求平均。”

符号含义

符号 含义
\(X\) 随机变量,例如一次随机生成的状态
\(p(x)\) \(X\) 在位置 \(x\) 附近出现的概率密度
\(f(x)\) 我们真正关心的量,例如收益、损失或函数值
\(\mathbb{E}_{p}\) 按照分布 \(p\) 取平均
\(\mu\) 想得到的真实期望

直觉解释

积分 \(\int f(x)p(x)\,dx\) 可以理解成:

  1. 每个位置 \(x\) 有一个结果 \(f(x)\)
  2. 这个位置出现的可能性由 \(p(x)\) 决定;
  3. 用出现概率作为权重,把所有结果加权平均。

当这个积分难以解析计算时,可以直接模拟“按 \(p\) 出现”的样本。

1.3 样本平均为什么能近似真实平均

\(p(x)\) 中独立采样:

\[X_1,X_2,\ldots,X_N\overset{\mathrm{iid}}{\sim}p(x), \]

然后计算:

\[\hat{\mu}_N = \frac{1}{N}\sum_{i=1}^{N}f(X_i). \]

怎么读

“把 \(N\) 次模拟得到的函数值相加,再除以 \(N\)。”

上标“帽子”表示它是估计值,不是真实值。下标 \(N\) 表示它由 \(N\) 个样本计算得到。

为什么成立

在期望和方差满足常见有限性条件时,大数定律保证:

\[\hat{\mu}_N\longrightarrow\mu, \qquad N\longrightarrow\infty. \]

这不是说有限样本时一定准确,而是说样本越多,偏离真实值很远的可能性通常越小。

1.4 误差为什么下降得慢

\(\operatorname{Var}[f(X)]=\sigma^2<\infty\),样本平均的方差为:

\[\operatorname{Var}(\hat{\mu}_N)=\frac{\sigma^2}{N}. \]

因此标准误差为:

\[\operatorname{SE}(\hat{\mu}_N)=\frac{\sigma}{\sqrt{N}}. \]

怎么读

“样本数增加 \(N\) 倍,误差只缩小到原来的 \(1/\sqrt{N}\)。”

例如:

  • 样本数乘以 \(4\),标准误差约减半;
  • 样本数乘以 \(100\),标准误差约缩小到十分之一;
  • 想多得到一位十进制精度,往往需要大约 \(100\) 倍样本。

这就是蒙特卡洛方法常常需要方差缩减、重要性采样、QMC 或 MLMC 的原因。

1.5 一套统一的观察框架

几乎所有直接相关的蒙特卡洛算法,都可以用以下五个问题理解:

  1. 样本是什么? 状态、参数、粒子、路径还是树节点?
  2. 样本从哪里来? 独立分布、提议分布、状态转移还是搜索策略?
  3. 样本是否带权? 每个样本贡献相同,还是有不同权重?
  4. 样本如何汇总? 求均值、加权平均、反向传播还是分层相加?
  5. 如何判断结果可信? 标准误差、ESS、收敛诊断还是访问次数?

2. 直接相关算法的整体地图

分支 核心思想 典型算法 核心数据结构
普通蒙特卡洛 独立采样后求平均 Monte Carlo Integration 样本集合、在线统计
采样改进 改变采样方式以提高效率 重要性采样、拒绝采样、分层采样 带权样本集合、接受样本缓冲区
MCMC 构造以目标分布为平稳分布的链 MH、Gibbs、HMC、NUTS 马尔可夫链轨迹
SMC 用带权粒子逐步逼近分布 粒子滤波、序贯重要性采样 加权粒子集合
MCTS 用随机模拟估计动作长期价值 UCT、PUCT 蒙特卡洛搜索树或图
QMC 用低差异序列均匀覆盖空间 Sobol、Halton 低差异序列生成器
MLMC 用多层差值替代单一高精度估计 多层蒙特卡洛 分层样本统计

最核心的四个分支可以记成:

\[\text{MC:独立样本} \quad \text{MCMC:相关链样本} \quad \text{SMC:带权粒子} \quad \text{MCTS:带统计的搜索树}. \]


3. 普通蒙特卡洛:用样本平均代替真实平均

3.1 蒙特卡洛积分

要计算一维积分:

\[I=\int_a^b g(x)\,dx. \]

\(X\) 在区间 \([a,b]\) 上均匀分布。均匀分布的密度为 \(1/(b-a)\),于是:

\[\mathbb{E}[g(X)] = \int_a^b g(x)\frac{1}{b-a}\,dx = \frac{I}{b-a}. \]

把等式重新排列:

\[I=(b-a)\mathbb{E}[g(X)]. \]

所以蒙特卡洛估计量为:

\[\hat I_N = \frac{b-a}{N}\sum_{i=1}^{N}g(X_i), \qquad X_i\sim U(a,b). \]

怎么读

“在区间内随机取点,计算这些点上的函数值平均数,再乘区间长度。”

小例子:计算 \(\int_0^1x^2\,dx\)

真实答案是 \(1/3\)。随机生成 \(N\)\([0,1]\) 内的点:

\[\hat I_N=\frac{1}{N}\sum_{i=1}^{N}X_i^2. \]

假设只抽到五个点:\(0.1,0.4,0.6,0.8,0.9\),则:

\[\hat I_5 = \frac{0.1^2+0.4^2+0.6^2+0.8^2+0.9^2}{5} =0.396. \]

五个样本很少,所以结果与 \(1/3\) 有明显差异。增加样本后通常会更稳定。

3.2 估计概率

事件 \(A\) 的概率可以写成指示函数的期望:

\[\mathbb{P}(X\in A) = \mathbb{E}[\mathbf{1}_{A}(X)]. \]

其中:

\[\mathbf{1}_{A}(X) = \begin{cases} 1, & X\in A,\\ 0, & X\notin A. \end{cases} \]

因此概率估计就是“命中次数除以总次数”:

\[\widehat{\mathbb{P}}(X\in A) = \frac{1}{N}\sum_{i=1}^{N}\mathbf{1}_{A}(X_i). \]

这解释了估计圆周率时为什么可以统计落在圆内的点的比例。

3.3 置信区间

当样本量足够大、方差有限且样本近似独立时,可用中心极限定理构造近似置信区间:

\[\hat\mu_N \pm z_{1-\alpha/2}\frac{s}{\sqrt N}. \]

符号含义

符号 含义
\(\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\) 时:

\[\delta=x_n-\bar x_{n-1}, \]

\[\bar x_n=\bar x_{n-1}+\frac{\delta}{n}, \]

\[M_{2,n}=M_{2,n-1}+\delta(x_n-\bar x_n). \]

最终样本方差为:

\[s^2=\frac{M_{2,n}}{n-1}. \]

直觉解释

Welford 算法不需要保存历史数据。每来一个新样本,就小幅修正均值,并记录这个新样本对总离散程度的贡献。它比“先求平方和再相减”的写法更稳定。

3.6 方差缩减

普通蒙特卡洛的关键瓶颈不是偏差,而往往是方差。常见方法包括:

  • 对偶变量:成对使用负相关样本,例如 \(U\)\(1-U\)
  • 控制变量:利用期望已知、且与目标量相关的变量修正估计;
  • 分层采样:把空间切成多个区域,避免样本偶然挤在某一处;
  • 条件蒙特卡洛:先解析消掉一部分随机性;
  • 拉丁超立方采样:让每个坐标方向的覆盖更均匀。

这些方法的共同目标不是“让样本更多”,而是“让每个样本提供更多有效信息”。


4. 重要性采样与拒绝采样

4.1 重要性采样:换一个更会抽重点的分布

目标仍然是:

\[\mu=\mathbb{E}_p[f(X)] =\int f(x)p(x)\,dx. \]

如果直接从 \(p\) 采样很困难,或真正重要的区域在 \(p\) 下很少出现,可以引入更容易采样的分布 \(q\)

\[\mu = \int f(x)\frac{p(x)}{q(x)}q(x)\,dx = \mathbb{E}_q\left[f(X)\frac{p(X)}{q(X)}\right]. \]

怎么读

“虽然样本来自 \(q\),但用 \(p/q\) 进行补偿,仍然可以计算在 \(p\) 下的平均。”

权重定义为:

\[w_i=\frac{p(X_i)}{q(X_i)}, \qquad X_i\sim q. \]

估计量为:

\[\hat\mu_N = \frac{1}{N}\sum_{i=1}^{N}w_i f(X_i). \]

为什么要乘 \(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)\) 就不能为零。否则那个区域永远抽不到,再大的权重也无法修正。

形式上要求:

\[p(x)|f(x)|>0\quad\Longrightarrow\quad q(x)>0. \]

4.2 自归一化重要性采样

有时只知道目标密度的未归一化形式:

\[\tilde p(x)=Zp(x), \]

其中未知常数 \(Z\) 无法计算。此时使用未归一化权重:

\[\tilde w_i=\frac{\tilde p(X_i)}{q(X_i)}, \]

再归一化:

\[\bar w_i = \frac{\tilde w_i}{\sum_{j=1}^{N}\tilde w_j}. \]

估计量为:

\[\hat\mu_{\mathrm{SNIS}} = \sum_{i=1}^{N}\bar w_i f(X_i). \]

直觉解释

未知常数 \(Z\) 同时出现在所有权重中。归一化时,分子分母中的 \(Z\) 会相互抵消,因此不需要知道它。

严谨性提醒

自归一化估计在有限样本下通常有偏,但在常见正则条件下是一致的,即样本数增加时趋近真实值。

4.3 有效样本量 ESS

若一个样本几乎占据全部权重,其他样本就几乎没有贡献。常用近似指标是:

\[\operatorname{ESS} = \frac{1}{\sum_{i=1}^{N}\bar w_i^2}. \]

怎么看这个公式

  • 若所有权重相等,\(\bar w_i=1/N\)

\[\operatorname{ESS} = \frac{1}{N(1/N)^2} =N. \]

这表示 \(N\) 个样本都充分有效。

  • 若一个权重为 \(1\),其余为 \(0\)

\[\operatorname{ESS}=1. \]

这表示虽然存了 \(N\) 个样本,但实际上只有一个样本在决定结果。

因此:

\[1\le \operatorname{ESS}\le 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\),满足:

\[p(x)\le Mq(x) \qquad\text{对所有 }x\text{ 成立}. \]

每次先采样:

\[X\sim q, \qquad U\sim U(0,1), \]

再按以下条件接受:

\[U\le\frac{p(X)}{Mq(X)}. \]

几何直觉

曲线 \(Mq(x)\) 像一张完全盖住 \(p(x)\) 的“外罩”。先在外罩下抽一个点,只有落在目标曲线 \(p(x)\) 下方时才接受。

\(p\)\(q\) 都已归一化,则平均接受率为:

\[\mathbb{P}(\text{接受})=\frac{1}{M}. \]

所以 \(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) =p(x\mid y) \propto p(y\mid x)p(x). \]

这里 \(\pi(x)\) 是目标后验分布。符号 \(\propto\) 表示“成比例”:我们能计算分子形状,但归一化常数可能很难求。

普通蒙特卡洛要求能直接从 \(\pi\) 采样。MCMC 的办法是:构造一条随机移动的链,让它长期停留频率恰好符合 \(\pi\)

5.2 马尔可夫链是什么

马尔可夫性质写成:

\[\mathbb{P}(X_{t+1}\mid X_t,X_{t-1},\ldots,X_0) = \mathbb{P}(X_{t+1}\mid X_t). \]

怎么读

“知道现在的状态后,预测下一步不再需要更早的历史。”

这不代表历史完全没有影响,而是历史影响已经被当前状态概括。

5.3 Metropolis–Hastings 接受率

当前状态为 \(x\),先从提议分布生成候选:

\[x'\sim q(x'\mid x). \]

接受概率为:

\[\alpha(x,x') = \min\left( 1, \frac{\pi(x')q(x\mid x')} {\pi(x)q(x'\mid x)} \right). \]

这是 MCMC 中最容易觉得晦涩的公式,可以拆成两部分:

\[\frac{\pi(x')}{\pi(x)} \times \frac{q(x\mid x')}{q(x'\mid x)}. \]

第一部分:目标分布比值

\[\frac{\pi(x')}{\pi(x)} \]

它比较候选状态和当前状态谁更符合目标分布。

  • 若候选更可能,比例大于 \(1\),通常直接接受;
  • 若候选更不可能,仍可能以一定概率接受,避免链被局部高点困住。

第二部分:提议方向修正

\[\frac{q(x\mid x')}{q(x'\mid x)} \]

如果从 \(x\)\(x'\) 很容易,而从 \(x'\) 回到 \(x\) 很难,那么单靠目标概率比会产生方向偏差。这一项负责修正提议机制的不对称。

对称提议的简化

若:

\[q(x'\mid x)=q(x\mid x'), \]

提议修正相互抵消,得到 Metropolis 接受率:

\[\alpha(x,x') = \min\left(1,\frac{\pi(x')}{\pi(x)}\right). \]

一个数值例子

若候选状态的目标密度只有当前状态的一半,而且提议对称:

\[\alpha=\min(1,0.5)=0.5. \]

这意味着候选虽然更差,但仍有一半概率被接受。允许“偶尔下坡”正是 MCMC 能穿过低概率区域的原因之一。

5.4 为什么这个接受率能工作:详细平衡

一个常见的充分条件是详细平衡:

\[\pi(x)K(x,x') = \pi(x')K(x',x), \]

其中 \(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 采样依次更新每个分量:

\[x_j^{(t+1)} \sim p\left( x_j \mid x_1^{(t+1)},\ldots,x_{j-1}^{(t+1)}, x_{j+1}^{(t)},\ldots,x_d^{(t)} \right). \]

怎么读

更新第 \(j\) 个变量时:

  • 已经在本轮更新过的变量,使用新值;
  • 尚未更新的变量,使用上一轮的旧值;
  • 只随机生成当前这个变量。

直觉解释

直接在高维联合分布里抽样可能很难,但固定其他变量后,每一个条件分布可能很简单。Gibbs 采样通过“逐坐标刷新”完成整体采样。

5.7 HMC:把随机游走变成有方向的运动

HMC 给参数 \(x\) 引入辅助动量 \(r\),定义哈密顿量:

\[H(x,r)=U(x)+K(r). \]

其中:

\[U(x)=-\log\pi(x), \]

\[K(r)=\frac{1}{2}r^{\mathsf T}M^{-1}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

概念上的链是:

\[x^{(0)}\rightarrow x^{(1)}\rightarrow\cdots\rightarrow x^{(T)}. \]

但程序中通常使用连续数组或张量,而不是链表。这里的“链”指概率依赖关系,不是编程语言中的链式存储。

5.9 自相关与 MCMC 有效样本量

MCMC 相邻样本相关,因此 \(T\) 次迭代不等于 \(T\) 个独立样本。一个常见关系是:

\[\operatorname{ESS}_{\mathrm{MCMC}} \approx \frac{T} {1+2\sum_{k=1}^{\infty}\rho_k}, \]

其中 \(\rho_k\) 是间隔 \(k\) 的自相关。

怎么读

分母是“相关性惩罚”。若自相关都接近零,ESS 接近 \(T\);若正相关很强,ESS 会远小于 \(T\)

实际软件会使用截断、谱估计或多链方法稳定估计 ESS,因此不应直接把无限和式机械实现。

5.10 MCMC 的关键风险

  • 预热不足:链仍受初始状态强烈影响;
  • 混合缓慢:链在参数空间移动太慢;
  • 多峰困境:链长期停留在一个模态;
  • 步长不合适:过小则移动慢,过大则接受率低;
  • 只看轨迹图不够:还应结合多链诊断、ESS 和模型检查。

6. SMC 与粒子滤波:用一群带权粒子追踪分布

6.1 状态空间模型

动态系统通常包含看不见的隐藏状态 \(x_t\) 和可观察数据 \(y_t\)

\[x_t\sim p(x_t\mid x_{t-1}), \]

\[y_t\sim p(y_t\mid x_t). \]

怎么读

  • 第一式:下一时刻状态由上一时刻状态随机演化;
  • 第二式:观测由当前隐藏状态随机生成。

目标是在看到 \(y_1,\ldots,y_t\) 后估计:

\[p(x_t\mid y_{1:t}). \]

这里 \(y_{1:t}\) 表示从第 \(1\) 时刻到第 \(t\) 时刻的全部观测。

6.2 用粒子表示一个分布

粒子滤波使用带权点集合近似后验分布:

\[p(x_t\mid y_{1:t}) \approx \sum_{i=1}^{N} w_t^{(i)}\, \delta\left(x_t-x_t^{(i)}\right). \]

这个公式怎么理解

它看起来抽象,其实表达的是:

“我们不再保存一条连续的概率密度曲线,而是用 \(N\) 个带权代表点来近似它。”

符号 含义
\(x_t^{(i)}\) \(i\) 个粒子在时刻 \(t\) 的状态
\(w_t^{(i)}\) 这个粒子的可信权重,通常满足权重和为 \(1\)
\(\delta(\cdot)\) 狄拉克点质量,表示概率集中在该粒子位置

不要把 \(\delta\) 理解成普通函数。这里它只是一种紧凑写法,用于表达“离散加权点集所代表的经验分布”。

6.3 Bootstrap 粒子滤波的权重更新

先传播粒子:

\[x_t^{(i)} \sim p(x_t\mid x_{t-1}^{(i)}). \]

然后用当前观测计算未归一化权重:

\[\tilde w_t^{(i)} = w_{t-1}^{(i)} p(y_t\mid x_t^{(i)}). \]

归一化:

\[w_t^{(i)} = \frac{\tilde w_t^{(i)}} {\sum_{j=1}^{N}\tilde w_t^{(j)}}. \]

直觉解释

  • 状态转移模型负责“预测下一步可能在哪里”;
  • 观测似然负责“根据新证据给预测打分”;
  • 归一化后,所有粒子的权重重新组成一个概率分布。

6.4 为什么需要重采样

多次更新后,常出现少数粒子权重很大,大多数权重接近零。此时计算资源被浪费在几乎没有贡献的粒子上。

使用:

\[\operatorname{ESS} = \frac{1}{\sum_i(w_t^{(i)})^2} \]

判断退化程度。当 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 通过重复模拟,把计算预算集中到更有希望的分支。

一次迭代分为:

  1. 选择:沿当前树选择值得继续探索的节点;
  2. 扩展:创建一个尚未尝试的子节点;
  3. 模拟:从新节点出发,用随机或启发式策略获得结果;
  4. 反向传播:把结果沿路径更新回去。

7.2 节点统计量

对节点或动作边,常保存:

\[Q_j=\frac{W_j}{N_j}. \]

符号 含义
\(N_j\) 被访问的次数
\(W_j\) 所有模拟收益之和
\(Q_j\) 平均模拟收益

直觉解释

\(Q_j\) 是“目前看起来有多好”,\(N_j\) 是“这个判断建立在多少次试验上”。一个高收益但只试过一次的节点,可信度通常不如一个试过很多次且收益稳定的节点。

7.3 UCT:利用与探索的平衡

UCT 评分为:

\[\operatorname{UCT}_j = \underbrace{\frac{W_j}{N_j}}_{\text{利用:目前平均收益}} + \underbrace{C\sqrt{\frac{\ln N_{\mathrm{parent}}}{N_j}}}_{\text{探索:对少访问节点的奖励}}. \]

逐项解释

利用项

\[\frac{W_j}{N_j} \]

优先选择目前平均收益高的节点。

探索项

\[C\sqrt{\frac{\ln N_{\mathrm{parent}}}{N_j}} \]

  • 父节点被访问得越多,仍未充分尝试的子节点越值得补试;
  • 子节点访问次数 \(N_j\) 越小,探索奖励越大;
  • \(C\) 控制探索强度。

为什么要取对数和平方根

这来自多臂老虎 JI的上置信界思想:探索奖励应随总试验次数缓慢增长,并随某个选项自身试验次数增加而衰减。

未访问节点怎么办

\(N_j=0\) 时公式会除以零。实现中通常把未访问节点的分数视为正无穷,或在进入 UCT 比较前优先访问每个动作至少一次。

7.4 PUCT:加入策略先验

PUCT 常写为:

\[\operatorname{PUCT}(s,a) = Q(s,a) + c_{\mathrm{puct}}P(s,a) \frac{\sqrt{N(s)}}{1+N(s,a)}. \]

符号含义

符号 含义
\(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\) 上:

\[I=\int_{[0,1]^d}f(x)\,dx, \]

QMC 使用确定性点列 \(x_1,\ldots,x_N\)

\[\hat I_N^{\mathrm{QMC}} = \frac{1}{N}\sum_{i=1}^{N}f(x_i). \]

公式和普通蒙特卡洛的样本平均几乎相同,区别在于点的产生方式:

  • MC:点是独立伪随机样本;
  • QMC:点是按低差异规则构造的有序序列。

8.3 “低差异”是什么意思

可以把单位区域切成许多小盒子。理想情况下,每个盒子里的点数应该接近“盒子体积 × 总点数”。差异度衡量实际点数与理想点数之间的最大偏差。

一个经典理论界是 Koksma–Hlawka 不等式:

\[\left| \frac{1}{N}\sum_{i=1}^{N}f(x_i) - \int_{[0,1]^d}f(x)\,dx \right| \le V_{\mathrm{HK}}(f)D_N^*. \]

不必被符号吓到

它只表达一件事:

积分误差不超过“函数有多难积分”乘以“点分布有多不均匀”。

  • \(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_1-P_0) +(P_2-P_1) +\cdots +(P_L-P_{L-1}). \]

中间项会相互抵消。例如前三层:

\[P_0+(P_1-P_0)+(P_2-P_1)=P_2. \]

对两边取期望:

\[\mathbb{E}[P_L] = \mathbb{E}[P_0] + \sum_{\ell=1}^{L} \mathbb{E}[P_\ell-P_{\ell-1}]. \]

直觉解释

不直接反复估计昂贵的 \(P_L\),而是:

  1. 大量估计便宜的粗层 \(P_0\)
  2. 分别估计相邻精度之间的小修正;
  3. 把所有修正加起来。

若相邻层使用耦合随机数,使 \(P_\ell\)\(P_{\ell-1}\) 高度相关,那么差值 \(P_\ell-P_{\ell-1}\) 的方差通常很小,因此高精度层只需要少量样本。

9.3 MLMC 估计量

令:

\[Y_0=P_0, \qquad Y_\ell=P_\ell-P_{\ell-1}\quad(\ell\ge1). \]

每层用 \(N_\ell\) 个样本估计:

\[\hat Y_\ell = \frac{1}{N_\ell} \sum_{i=1}^{N_\ell}Y_\ell^{(i)}. \]

总估计为:

\[\hat P_{\mathrm{MLMC}} = \sum_{\ell=0}^{L}\hat Y_\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

加权均值:

\[\bar x_w=\sum_{i=1}^{N}\bar w_i x_i, \qquad \sum_i\bar w_i=1. \]

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. 加权抽样的辅助数据结构

重要性采样和粒子重采样经常需要按离散权重:

\[w_1,w_2,\ldots,w_N, \qquad w_i\ge0, \]

随机选择索引 \(i\)

11.1 前缀和数组

构造:

\[C_i=\sum_{j=1}^{i}w_j. \]

生成 \(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

把每个离散概率拆成“主选项 + 备用选项”,使采样只需:

  1. 均匀选择一个桶;
  2. 再做一次伯努利判断。
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 对数域权重

直接计算:

\[w_i=\frac{p(x_i)}{q(x_i)} \]

可能因概率极小而下溢。改为:

\[\log w_i = \log p(x_i)-\log q(x_i). \]

归一化时使用 Log-Sum-Exp。设 \(a_i=\log w_i\),取:

\[m=\max_i a_i. \]

则:

\[\log\sum_i e^{a_i} = m+ \log\sum_i e^{a_i-m}. \]

为什么有效

\(a_i-m\le0\),所以指数值不会爆炸;最大的那一项变为 \(e^0=1\),也不容易全部下溢为零。

归一化对数权重可写为:

\[\log\bar w_i = a_i- \operatorname{LSE}(a_1,\ldots,a_N). \]

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 算法层面

直接与蒙特卡洛相关的主要算法是:

\[\text{普通 MC} + \text{重要性与拒绝采样} + \text{MCMC} + \text{SMC} + \text{MCTS} + \text{QMC} + \text{MLMC}. \]

它们的本质区别可以压缩成一句话:

  • 普通 MC:生成独立样本后求平均;
  • 重要性采样:从别处分布抽样,再用权重修正;
  • 拒绝采样:从包络分布提候选,只保留合格样本;
  • MCMC:通过有依赖的随机链间接获得目标分布样本;
  • SMC:让一群带权粒子随时间传播、加权和重采样;
  • MCTS:把随机模拟结果累计在搜索树节点或边上;
  • QMC:用覆盖更均匀的序列计算样本平均;
  • MLMC:把昂贵高精度答案拆成便宜粗层和小修正。

16.2 数据结构层面

直接绑定这些算法的核心结构是:

\[\text{样本集合} + \text{带权样本集合} + \text{MCMC 轨迹} + \text{粒子集合} + \text{搜索树或图} + \text{低差异序列状态} + \text{分层统计}. \]

最重要的对应关系:

分支 主要对象 必须保存的核心信息
普通 MC 样本或在线统计 数量、均值、方差
重要性采样 带权样本 对数权重、归一化权重、ESS
MCMC 有序轨迹 状态、接受标记、概率与诊断
SMC 粒子群 状态、权重、祖先、时间步
MCTS 搜索节点和边 访问次数、累计价值、先验
QMC 序列生成器 维度、索引、方向数、扰乱状态
MLMC 多层统计 各层样本数、差值方差、成本

最终可以把“蒙特卡洛数据结构”理解为:

为了生成、保存、更新、重采样和汇总蒙特卡洛样本而设计的算法状态结构。

它不是一种单独的数据结构类别,而是一组与不同蒙特卡洛算法紧密绑定的表示方式。

posted @ 2026-08-03 21:05  aixueforever  阅读(7)  评论(0)    收藏  举报