第十七章：蒙特卡洛方法

1. 采样与概率基础

1.1 蒙特卡洛方法的哲学

蒙特卡洛方法（Monte Carlo Methods）是一类通过随机采样来近似计算确定性问题的方法。其核心思想源自大数定律：当样本数量趋于无穷时，样本均值将收敛到期望值。这一特性使得我们能够通过有限样本近似难以精确求解的概率和积分问题。

在深度学习中，蒙特卡洛方法扮演着至关重要的角色。许多深度生成模型（如变分自编码器、玻尔兹曼机）的训练过程本质上涉及对难以解析计算的概率分布进行采样与估计。

1.2 期望的蒙特卡洛估计

给定一个随机变量 \(\mathbf{x} \sim p(\mathbf{x})\)，我们希望计算函数 \(f(\mathbf{x})\) 的期望：

\[\mathbb{E}_{p(\mathbf{x})}[f(\mathbf{x})] = \int f(\mathbf{x}) p(\mathbf{x}) \, d\mathbf{x}\]

蒙特卡洛估计器通过从 \(p(\mathbf{x})\) 中独立抽取 \(N\) 个样本 \(\{\mathbf{x}^{(i)}\}_{i=1}^{N}\) 来近似这个期望：

\[\widehat{\mathbb{E}}_{N} = \frac{1}{N} \sum_{i=1}^{N} f(\mathbf{x}^{(i)})\]

根据大数定律，当 \(N \to \infty\) 时，\(\widehat{\mathbb{E}}_{N} \to \mathbb{E}[f(\mathbf{x})]\) 几乎必然收敛。估计的方差为：

\[\text{Var}(\widehat{\mathbb{E}}_{N}) = \frac{\text{Var}(f(\mathbf{x}))}{N}\]

这表明估计误差的衰减速度是 \(O(1/\sqrt{N})\)——样本数量增加到原来的4倍，误差才减半。

1.3 采样问题的本质

采样的核心任务是：从已知但难以直接处理的概率分布 \(p(\mathbf{x})\) 中生成独立同分布的样本。直接采样的障碍通常来自：

• 归一化常数未知：\(p(\mathbf{x}) = \tilde{p}(\mathbf{x}) / Z\)，其中 \(Z\) 难以计算
• 多峰分布：简单的拒绝采样效率极低
• 高维空间：维数灾难使传统方法失效

---

2. 随机数生成方法

2.1 逆变换法（Inverse Transform Method）

逆变换法是最基本的采样技术，适用于已知累积分布函数（CDF）的情况。

设随机变量 \(X\) 的CDF为 \(F_X(x) = P(X \leq x)\)，其逆函数为 \(F_X^{-1}(u)\)。若 \(U \sim \text{Uniform}(0, 1)\)，则：

\[X = F_X^{-1}(U)\]

推导：由于 \(P(F_X^{-1}(U) \leq x) = P(U \leq F_X(x)) = F_X(x)\)，故 \(F_X^{-1}(U)\) 的分布正是 \(X\) 的分布。

示例——指数分布： 指数分布 \(p(x) = \lambda e^{-\lambda x}\)（\(x \geq 0\)）的CDF为 \(F(x) = 1 - e^{-\lambda x}\)，其逆函数为：

\[F^{-1}(u) = -\frac{\ln(1-u)}{\lambda}\]

当 \(u \sim \text{Uniform}(0,1)\) 时，\(x = -\ln(1-u)/\lambda \sim \text{Exp}(\lambda)\)。

适用范围：需要知道逆CDF的解析形式，对于复杂分布往往不可行。

2.2 拒绝采样（Rejection Sampling）

拒绝采样允许我们从难以直接采样的分布 \(p(x)\) 中获取样本，条件是存在一个易采样的建议分布（proposal distribution） \(q(x)\)，且满足：

\[p(x) \leq M \cdot q(x), \quad \forall x\]

其中常数 \(M \geq 1\)（通常希望 \(M\) 尽可能小以提高效率）。

算法步骤：

1. 从 \(q(x)\) 中抽取样本 \(x\)
2. 从 \(\text{Uniform}(0, 1)\) 中抽取 \(u\)
3. 若 \(u < p(x) / (M q(x))\)，接受 \(x\)；否则拒绝并回到步骤1

接受概率：每次迭代的接受概率为 \(1/M\)（理论上），实际为：

\[P(\text{accept}) = \int \frac{p(x)}{M q(x)} q(x) dx = \frac{1}{M} \int p(x) dx = \frac{1}{M}\]

效率分析：当 \(p(x)\) 和 \(q(x)\) 差异较大时，\(M\) 必须很大，导致接受率极低。这在高维问题中尤为严重——\(M\) 往往随维数指数增长，使得拒绝采样实际不可用。

2.3 重要性采样（Importance Sampling）

重要性采样通过加权平均来估计期望，而非直接拒绝样本。

设真实分布为 \(p(\mathbf{x})\)，建议分布为 \(q(\mathbf{x})\)，则：

\[\mathbb{E}_p[f(\mathbf{x})] = \int f(\mathbf{x}) p(\mathbf{x}) d\mathbf{x} = \int f(\mathbf{x}) \frac{p(\mathbf{x})}{q(\mathbf{x})} q(\mathbf{x}) d\mathbf{x}\]

定义重要性权重：\(w(\mathbf{x}) = p(\mathbf{x}) / q(\mathbf{x})\)，则：

\[\mathbb{E}_p[f(\mathbf{x})] \approx \sum_{i=1}^{N} \frac{w(\mathbf{x}^{(i)})}{\sum_{j=1}^{N} w(\mathbf{x}^{(j)})} f(\mathbf{x}^{(i)})\]

自归一化重要性采样（SNIS）：当 \(p(\mathbf{x})\) 只知其未归一化形式 \(\tilde{p}(\mathbf{x})\) 时：

\[\widehat{\mathbb{E}}_{p}[f] = \sum_{i=1}^{N} \tilde{w}(\mathbf{x}^{(i)}) f(\mathbf{x}^{(i)}) / \sum_{j=1}^{N} \tilde{w}(\mathbf{x}^{(j)})\]

其中 \(\tilde{w}(\mathbf{x}^{(i)}) = \tilde{p}(\mathbf{x}^{(i)}) / q(\mathbf{x}^{(i)})\)。

最优建议分布：理论上，最优建议分布为 \(q^*(\mathbf{x}) \propto |f(\mathbf{x})| p(\mathbf{x})\)，此时方差为零——但这需要从 \(p\) 本身采样，本质上是循环论证。

应用场景：重要性采样广泛用于深度学习中的权重更新和策略梯度估计，特别是在强化学习中。

---

3. 马尔可夫链蒙特卡洛方法

3.1 马尔可夫链基础

马尔可夫链蒙特卡洛（MCMC）是一类通过构建马尔可夫链来获取高维分布样本的方法。其核心思想是：设计一个马尔可夫链，使其平稳分布恰好是我们希望采样的目标分布 \(p(\mathbf{x})\)。

马尔可夫链定义：考虑离散状态空间 \(\mathcal{X}\)，转移概率 \(T(\mathbf{x}'|\mathbf{x})\) 表示从状态 \(\mathbf{x}\) 转移到 \(\mathbf{x}'\) 的概率。

平稳分布：若分布 \(\pi(\mathbf{x})\) 满足：

\[\pi(\mathbf{x}') = \sum_{\mathbf{x}} T(\mathbf{x}'|\mathbf{x}) \pi(\mathbf{x})\]

则称 \(\pi\) 为该马尔可夫链的平稳分布。

细致平衡条件（Detailed Balance）：若对所有 \(\mathbf{x}, \mathbf{x}'\) 有：

\[T(\mathbf{x}'|\mathbf{x}) \pi(\mathbf{x}) = T(\mathbf{x}|\mathbf{x}') \pi(\mathbf{x}')\]

则平稳条件必然满足。这是构造MCMC方法的关键。

3.2 Metropolis-Hastings算法

Metropolis-Hastings算法是最广泛使用的MCMC方法之一，通过接受-拒绝机制确保细致平衡。

算法步骤：

1. 给定当前状态 \(\mathbf{x}^{(t)}\)
2. 从建议分布 \(q(\mathbf{x}'|\mathbf{x}^{(t)})\) 中抽取候选样本 \(\mathbf{x}'\)
3. 计算接受概率：

\[A = \min\left(1, \frac{p(\mathbf{x}') q(\mathbf{x}^{(t)}|\mathbf{x}')}{p(\mathbf{x}^{(t)}) q(\mathbf{x}'|\mathbf{x}^{(t)})}\right)\]

1. 从 \(\text{Uniform}(0,1)\) 抽取 \(u\)，若 \(u < A\) 则接受 \(\mathbf{x}'\)，否则保持 \(\mathbf{x}^{(t)}\)

对称建议分布：当 \(q(\mathbf{x}'|\mathbf{x}) = q(\mathbf{x}|\mathbf{x}')\)（对称）时，接受概率简化为：

\[A = \min\left(1, \frac{p(\mathbf{x}')}{p(\mathbf{x}^{(t)})}\right)\]

这即是原始Metropolis算法的形式。

遍历性条件：为确保马尔可夫链收敛到唯一平稳分布，需要满足： - 不可约性：从任意状态可到达任意其他状态 - 非周期性：避免周期性振荡

3.3 吉布斯采样（Gibbs Sampling）

吉布斯采样是一种特殊的Metropolis-Hastings算法，其接受率始终为1，适合于多维分布的条件采样。

核心思想：在第 \(t\) 次迭代时，依次对每个维度进行条件采样：

\[x_1^{(t+1)} \sim p(x_1 | x_2^{(t)}, x_3^{(t)}, \ldots, x_d^{(t)})$$ $$x_2^{(t+1)} \sim p(x_2 | x_1^{(t+1)}, x_3^{(t)}, \ldots, x_d^{(t)})$$ $$\vdots$$ $$x_d^{(t+1)} \sim p(x_d | x_1^{(t+1)}, x_2^{(t+1)}, \ldots, x_{d-1}^{(t+1)})\]

接受率证明：令 \(\mathbf{x} = (x_i, \mathbf{x}_{-i})\)，建议分布为条件分布 \(q(\mathbf{x}'|\mathbf{x}) = p(x_i'|\mathbf{x}_{-i})\)，则：

\[\frac{p(\mathbf{x}') q(\mathbf{x}|\mathbf{x}')}{p(\mathbf{x}) q(\mathbf{x}'|\mathbf{x})} = \frac{p(x_i'|\mathbf{x}_{-i}) p(\mathbf{x}_{-i})}{p(x_i|\mathbf{x}_{-i}) p(\mathbf{x}_{-i})} \cdot \frac{p(x_i|\mathbf{x}_{-i})}{p(x_i'|\mathbf{x}_{-i})} = 1\]

因此 \(A = \min(1,1) = 1\)，接受率为100%。

应用场景：吉布斯采样在贝叶斯统计中广泛应用，尤其当条件分布具有简洁形式时。

3.4 哈密顿蒙特卡洛（Hamiltonian Monte Carlo）

哈密顿蒙特卡洛（HMC）利用梯度信息来提高采样效率，减少随机游走行为。

物理类比：将采样过程类比为哈密顿动力学系统。引入辅助动量变量 \(\mathbf{r}\)，定义哈密顿量：

\[H(\mathbf{x}, \mathbf{r}) = U(\mathbf{x}) + K(\mathbf{r}) = -\log p(\mathbf{x}) + \frac{1}{2} \mathbf{r}^T \mathbf{M}^{-1} \mathbf{r}\]

其中 \(U(\mathbf{x})\) 为势能，\(K(\mathbf{r})\) 为动能，\(\mathbf{M}\) 为质量矩阵（通常取单位矩阵）。

更新方程（Leapfrog积分）：

\[\mathbf{r}(t + \epsilon/2) = \mathbf{r}(t) - \frac{\epsilon}{2} \nabla U(\mathbf{x}(t))$$ $$\mathbf{x}(t + \epsilon) = \mathbf{x}(t) + \epsilon \mathbf{M}^{-1} \mathbf{r}(t + \epsilon/2)$$ $$\mathbf{r}(t + \epsilon) = \mathbf{r}(t + \epsilon/2) - \frac{\epsilon}{2} \nabla U(\mathbf{x}(t + \epsilon))\]

算法步骤：

1. 从标准正态分布抽取动量 \(\mathbf{r}\)
2. 执行 \(L\) 步Leapfrog积分得到候选状态 \((\mathbf{x}^*, \mathbf{r}^*)\)
3. 以概率 \(\min(1, \exp(H(\mathbf{x}, \mathbf{r}) - H(\mathbf{x}^*, \mathbf{r}^*)))\) 接受

优势：HMC利用梯度信息在高维空间中有效探索，避免随机游走的低效性。在深度学习的贝叶斯后验采样中表现优异。

---

4. 马尔可夫链诊断

4.1 混合时间（Mixing Time）

混合时间衡量马尔可夫链从初始分布收敛到平稳分布所需的时间。

总变差距离：

\[d(t) = \max_{\mathbf{x}} \| P^t(\mathbf{x}, \cdot) - \pi(\cdot) \|_{\text{TV}} = \frac{1}{2} \sum_{\mathbf{x}'} |P^t(\mathbf{x}, \mathbf{x}') - \pi(\mathbf{x}')|\]

其中 \(P^t(\mathbf{x}, \mathbf{x}')\) 为 \(t\) 步转移概率。

\(\epsilon\)-混合时间定义为：

\[\tau_\epsilon = \min\{t : d(t) \leq \epsilon\}\]

影响因素：

• 建议分布选择：好的建议分布加速混合
• 维度：高维空间通常导致更长的混合时间
• 分布形态：多峰分布、窄通道等问题

实践中的处理：通常需要"预烧期（burn-in）"，丢弃前 \(\tau_\epsilon\) 个样本以确保链已收敛。

4.2 自相关函数（Autocorrelation Function）

由于MCMC生成的样本是相关的，不能视为独立同分布。自相关函数衡量滞后 \(k\) 的样本间的相关性。

自相关函数定义：

\[\rho_k = \frac{\text{Cov}(X_t, X_{t+k})}{\text{Var}(X_t)} = \frac{\sum_{t=1}^{N-k}(X_t - \bar{X})(X_{t+k} - \bar{X})}{\sum_{t=1}^{N}(X_t - \bar{X})^2}\]

有效样本量的关系：自相关时间 \(\tau_{\text{AC}}\) 定义为使 \(\rho_k\) 衰减到接近零的 \(k\) 值：

\[\tau_{\text{AC}} = 1 + 2\sum_{k=1}^{\infty} \rho_k\]

样本间的有效独立样本数量约为 \(N / \tau_{\text{AC}}\)。

4.3 有效样本大小（Effective Sample Size, ESS）

有效样本大小是衡量MCMC采样效率的关键指标，它将相关样本数折算为等效的独立样本数。

定义：

\[\text{ESS} = \frac{N}{\tau_{\text{AC}}} = N \cdot \frac{1}{1 + 2\sum_{k=1}^{\infty} \rho_k}\]

实践估计：由于无法计算无穷级数，通常使用截断估计：

\[\widehat{\text{ESS}} = \frac{N}{1 + 2\sum_{k=1}^{K} \rho_k}\]

其中 \(K\) 通常取使 \(\rho_k < 0.05\) 的最小值。

判断标准：

ESS值	评估
\(> 1000\)	充足
\(100 \sim 1000\)	可接受
\(< 100\)	不足，需改进采样

提高ESS的方法：减小样本间相关性的关键是改进马尔可夫链的混合速度——可考虑调整建议分布、使用HMC或并行链技术。

---

5. 方差缩减技术

5.1 对偶变量法（Antithetic Variates）

对偶变量法通过引入负相关的样本对来减少估计方差。

原理：若两个随机变量 \(X\) 和 \(X'\) 满足 \(\text{Cov}(X, X') < 0\)，则：

\[\text{Var}\left(\frac{X + X'}{2}\right) = \frac{\text{Var}(X) + \text{Var}(X') + 2\text{Cov}(X, X')}{4} < \frac{\text{Var}(X)}{2}\]

实现方式：对于逆变换法，若 \(U \sim \text{Uniform}(0,1)\)，则 \(U' = 1 - U\) 也服从均匀分布，且两者负相关。样本对为：

\[X = F^{-1}(U), \quad X' = F^{-1}(1-U)\]

估计器：

\[\widehat{\mu} = \frac{1}{2N} \sum_{i=1}^{N} \left[ F^{-1}(U_i) + F^{-1}(1-U_i) \right]\]

方差缩减效果：理论上可将方差减少数倍，但效果取决于分布的具体形式。

5.2 控制变量法（Control Variates）

控制变量法利用与目标变量相关的已知量来减少估计方差。

设 \(f\) 为目标函数，\(g\) 为控制变量，且 \(\mathbb{E}[g(\mathbf{x})]\) 已知。构造估计器：

\[\widehat{\mathbb{E}}[f] = \frac{1}{N} \sum_{i=1}^{N} \left[ f(\mathbf{x}^{(i)}) - c \cdot (g(\mathbf{x}^{(i)}) - \mathbb{E}[g])\right]\]

其中 \(c\) 为最优系数：

\[c^* = \frac{\text{Cov}(f, g)}{\text{Var}(g)}\]

最优方差：

\[\text{Var}(\widehat{\mathbb{E}}^*[f]) = \text{Var}(f) - \frac{\text{Cov}(f, g)^2}{\text{Var}(g)} = \text{Var}(f)(1 - \rho_{fg}^2)\]

其中 \(\rho_{fg}\) 为 \(f\) 与 \(g\) 的相关系数。当 \(|\rho_{fg}|\) 接近1时，方差可大幅缩减。

应用要点：控制变量的期望必须已知，且与目标函数高度相关。在金融衍生品定价和深度学习训练中应用广泛。

5.3 重要性加权（Importance Weighting）

如第二章所述，重要性采样本质上也是一种方差缩减技术。关键在于选择合适的建议分布 \(q(\mathbf{x})\)。

最优建议分布：理论上 \(q^*(\mathbf{x}) \propto |f(\mathbf{x})| p(\mathbf{x})\) 给出零方差，但不可行。实践中需在方差缩减与采样难度间权衡。

自适应重要性采样：通过迭代改进建议分布，使其逐步接近最优分布。常见方法包括：

• 基于样本密度估计：用核方法或高斯混合模型拟合加权样本
• 正则化流动：使用可逆神经网络变换简单分布

---

6. 公式汇总表

方法类别	公式名称	核心公式	应用场景
采样基础	蒙特卡洛估计	\(\widehat{\mathbb{E}}[f] = \frac{1}{N}\sum_{i=1}^{N} f(\mathbf{x}^{(i)})\)	期望估计
采样基础	估计方差	\(\text{Var}(\widehat{\mathbb{E}}_N) = \text{Var}(f)/N\)	误差分析
逆变换	逆CDF变换	\(X = F_X^{-1}(U),\ U\sim\text{U}(0,1)\)	简单分布采样
拒绝采样	接受概率	\(A = \min(1, p(\mathbf{x}')/(Mq(\mathbf{x}')))\)	一般分布采样
重要性采样	权重公式	\(w(\mathbf{x}) = p(\mathbf{x})/q(\mathbf{x})\)	期望估计
M-H算法	接受率	$A = \min(1, \frac{p(\mathbf{x}')q(\mathbf{x}	\mathbf{x}')}{p(\mathbf{x})q(\mathbf{x}'
吉布斯采样	条件采样	$x_i \sim p(x_i	\mathbf{x}_{-i})$
HMC	哈密顿量	\(H = -\log p(\mathbf{x}) + \frac{1}{2}\mathbf{r}^T\mathbf{M}^{-1}\mathbf{r}\)	高效采样
诊断	自相关时间	\(\tau_{\text{AC}} = 1+2\sum_{k=1}^{\infty}\rho_k\)	样本相关性
诊断	有效样本量	\(\text{ESS} = N/\tau_{\text{AC}}\)	采样效率
对偶变量	估计器	\(\widehat{\mu} = \frac{1}{2N}\sum[F^{-1}(U_i)+F^{-1}(1-U_i)]\)	方差缩减
控制变量	最优系数	\(c^* = \text{Cov}(f,g)/\text{Var}(g)\)	方差缩减

---

7. 深度学习应用

7.1 深度玻尔兹曼机（DBM）的训练

深度玻尔兹曼机（DBM）是一种基于能量的概率图模型，其配分函数 \(Z\) 通常难以计算。训练涉及对数似然的梯度：

\[\nabla_\theta \log p(\mathbf{v}; \theta) = \mathbb{E}_{p(\mathbf{h}|\mathbf{v})}[\nabla_\theta \log \tilde{p}(\mathbf{v}, \mathbf{h}; \theta)] - \mathbb{E}_{p(\mathbf{h}, \mathbf{v}; \theta)}[\nabla_\theta \log \tilde{p}(\mathbf{v}, \mathbf{h}; \theta)]\]

第一项为正相期望，可通过给定 \(\mathbf{v}\) 条件下的隐藏层采样高效计算；第二项为负相期望，涉及对 \(p(\mathbf{h}, \mathbf{v})\) 的采样，通常使用MCMC（特别是CD-k或PCD）近似。

对比散度（Contrastive Divergence, CD）：用少量MCMC步骤（通常 \(k=1\)）近似负相期望：

\[\mathbf{w} \leftarrow \mathbf{w} + \epsilon \left( \mathbb{E}_{p_{\text{data}}}[\mathbf{v}\mathbf{h}^T] - \mathbb{E}_{p_{\text{model}}^{(k)}}[\mathbf{v}\mathbf{h}^T] \right)\]

持续对比散度（PCD）：维护一个马尔可夫链的持久样本，避免每轮重启动，提高负相估计质量。

7.2 变分自编码器（VAE）

变分自编码器（VAE）通过变分推断近似后验分布 \(p(\mathbf{z}|\mathbf{x})\)，涉及对潜在变量 \(\mathbf{z}\) 的采样。

ELBO（证据下界）：

\[\mathcal{L}(\mathbf{x}; \theta, \phi) = \mathbb{E}_{q_\phi(\mathbf{z}|\mathbf{x})}[\log p_\theta(\mathbf{x}|\mathbf{z})] - D_{\text{KL}}(q_\phi(\mathbf{z}|\mathbf{x}) \| p(\mathbf{z}))\]

重参数化技巧：为使采样过程可导，引入辅助噪声 \(\boldsymbol{\epsilon} \sim \mathcal{N}(0, \mathbf{I})\)：

\[\mathbf{z} = \boldsymbol{\mu}_\phi(\mathbf{x}) + \boldsymbol{\sigma}_\phi(\mathbf{x}) \odot \boldsymbol{\epsilon}\]

这将随机采样嵌入到确定性的神经网络前向传播中，实现端到端梯度训练。

重要性采样增强：在计算重构似然时，可使用重要性采样提高估计精度：

\[\log p_\theta(\mathbf{x}) \approx \log \frac{1}{K} \sum_{k=1}^{K} \frac{p_\theta(\mathbf{x}, \mathbf{z}^{(k)})}{q_\phi(\mathbf{z}^{(k)}|\mathbf{x})}\]

7.3 Helmholtz机与唤醒-睡眠算法

Helmholtz机是最早的深度变分模型之一，引入辅助隐变量进行协同工作。

双循环结构： - 唤醒相（Wake phase）：固定 bottom-up 识别网络 \(q(\mathbf{h}|\mathbf{v})\)，采样 \(\mathbf{h}\)，更新生成网络 \(p(\mathbf{v}|\mathbf{h})\) - 睡眠相（Sleep phase）：固定 top-down 生成网络 \(p(\mathbf{v}, \mathbf{h})\)，采样 \((\mathbf{v}, \mathbf{h})\)，更新识别网络 \(q(\mathbf{h}|\mathbf{v})\)

训练目标：最大化似然的变分下界，等价于最小化识别分布与生成分布间的KL散度。

梯度估计：Helmholtz机使用采样估计梯度，涉及从 \(p(\mathbf{h}|\mathbf{v})\) 和 \(p(\mathbf{v}, \mathbf{h})\) 的双重采样。蒙特卡洛方法在此发挥关键作用。

与DBM的关系：Helmholtz机可视为DBM的前身，但其训练算法（唤醒-睡眠）不如DBM的CD算法高效。Wake-Sleep算法存在"睡觉问题"——识别网络可能学习到错误的生成分布。

---

结语

蒙特卡洛方法是深度学习不可或缺的理论工具，从最初的随机采样发展到如今精密的MCMC框架，涵盖了丰富的方差缩减与诊断技术。在深度生成模型日益重要的今天，理解这些方法的原理与局限性，对于设计和实现高效、可信的机器学习系统至关重要。随着计算资源的增长和算法的发展，蒙特卡洛方法将继续在贝叶斯深度学习、概率编程和生成模型等领域发挥核心作用。