采样——逆变换、重参数化与蒙特卡洛

上一课我们让三条线汇合：交叉熵 = 最大似然 = 最小化 \(D_{\mathrm{KL}}(p\|q)\)。那一课我们一直在"算期望、算 KL"，但有一个问题被悄悄跳过了：这些期望到底怎么算出来？ 当分布复杂到积不出闭式时，我们就只能采样（sampling）——造一堆样本，用样本平均替代积分。这一课就讲清楚：怎么从一个分布里造样本，怎么用样本估计期望，以及一个让 VAE / 扩散模型能反向传播的关键花招——重参数化（reparameterization）。

一、为什么要采样

采样这件事，在机器学习里干三种活：

• 生成样本：你训练了一个图像/文本/分子的分布模型，想"画"出新东西——就得从分布里采样。LLM 每吐一个 token 都是一次采样。
• 估计期望：损失、梯度、价值函数、KL……机器学习里几乎所有目标都是某个期望 \(\mathbb{E}[f(X)]\)。分布一复杂，积分就没有闭式，只能用样本平均近似。
• 随机算法：随机初始化、dropout、SGD 的 mini-batch、强化学习里按策略 \(\pi(a|s)\) 探索动作——本质都是采样。

二、蒙特卡洛估计：用平均代替积分

回忆本模块第1课：期望就是"按概率加权的平均"。连续情形 \(\mathbb{E}[f(X)]=\int f(x)\,p(x)\,dx\)。如果我们能从 \(p\) 里独立采出 \(x_1,\dots,x_N\)，那么自然的近似是直接求样本平均：

这就是蒙特卡洛（Monte Carlo, MC）估计。它为什么对？因为大数定律保证 \(\hat\mu_N\) 收敛到真值。但工程上更关心的是：收敛有多快？

三、逆变换采样：用 CDF 的逆把均匀数"掰"成任意分布

蒙特卡洛要求"能从 \(p\) 采样"。可计算机的随机数发生器只会吐均匀随机数 \(U\sim\mathrm{Uniform}(0,1)\)。怎么把它变成指数分布、高斯分布的样本？

先想机制，别背公式。 回忆本模块第1课的累积分布函数（CDF）\(F(x)=P(X\le x)\)：它是一条从 0 单调爬到 1 的曲线。关键观察：CDF 的取值本身就均匀分布在 \([0,1]\)。所以反过来——如果我先在 \(y\) 轴上均匀地取一个高度 \(U\)，再沿水平线打到曲线上，竖直落下读出 \(x\)，那么曲线陡的地方（概率密度大）落点密，曲线平的地方落点稀——这正好就是 \(p\) 想要的样本分布！

四、★重参数化技巧：让梯度穿过采样

这是本课的"主菜"，也是后面 VAE / 扩散模型反向传播能工作的根。先讲清它解决的痛点。

假设我们有一个高斯 \(\mathcal N(\mu,\sigma^2)\)，参数 \(\mu,\sigma\) 是网络要学的。训练目标常常长这样：

我们想用反向传播（模块1第10课）求 \(\partial L/\partial\mu\)。但问题来了：计算图里有一个"采样"节点——\(x\) 是从 \(\mathcal N(\mu,\sigma^2)\) 随机抽的。采样是个掷骰子的离散动作，不可导，梯度一到这里就断了，传不回 \(\mu,\sigma\)。

具体导数（呼应链式法则，模块1第8课）：把 \(\varepsilon\) 当常数，对 \(x=\mu+\sigma\varepsilon\) 求偏导，

这两个数清清楚楚、处处可导，链式法则可以接力：\(\dfrac{\partial f}{\partial\mu}=f'(x)\cdot 1\)，\(\dfrac{\partial f}{\partial\sigma}=f'(x)\cdot\varepsilon\)。

五、采样工具箱的其余三件：拒绝、重要性、MCMC

拒绝采样（rejection sampling）——"撒点，扔掉超出去的"

想从一个只知道形状 \(p(x)\)（甚至没归一化）的分布采样，但 \(F^{-1}\) 求不出来。办法：找一个容易采样的"罩子" \(M\,q(x)\ge p(x)\)，从 \(q\) 采一个候选 \(x\)，再以概率 \(p(x)/(M q(x))\) 接受它，否则扔掉重来。直觉上，当罩子 \(q\) 本身是均匀分布时，这就像往罩子覆盖的矩形里均匀撒点，落在 \(p\) 曲线下面的留、上面的扔，留下的点就服从 \(p\)（下面练习 3 正是这个最简特例）。\(q\) 非均匀时，"均匀撒点"的图景要换成"先按 \(q\) 撒、再按比值接受"，但留下的样本同样服从 \(p\)。缺点：罩子套得松（\(M\) 大）时拒绝率高、很浪费，高维尤其糟。

重要性采样（importance sampling）——"借别人的样本，按权重折算"

有时我们根本不需要 \(p\) 的样本，只想算 \(\mathbb E_p[f(X)]\)；偏偏从 \(p\) 采样很难，但有个好采样的 \(q\)。做个恒等变形：

于是从 \(q\) 采样，每个样本乘上重要性权重 \(w(x)=p(x)/q(x)\) 再平均即可：

直觉：\(q\) 在某处采多了（相对 \(p\)），权重 \(p/q<1\) 就把它调低；采少了就调高。校正"用错了分布采样"这件事。

MCMC（马尔可夫链蒙特卡洛）——高维相关采样

当分布是高维、强相关、且只知道未归一化形状时，前面几招都失灵。MCMC 换个思路：不直接采独立样本，而是造一条随机游走的链，让它的长期停留分布恰好是目标 \(p\)。Metropolis–Hastings：从当前点提议一个邻近点，按 \(\min(1,\,p_{\text{新}}/p_{\text{旧}})\) 的概率接受（注意只用到比值，未知的归一化常数在分子分母上下相消，自动消掉！）。Gibbs 采样：每次固定其他维、只按条件分布更新一维，轮流刷。代价是相邻样本相关、需要"烧入期（burn-in）"、且要判断收敛——这一块我们只点到为止。

六、随机种子与可复现，以及 LLM 解码也是采样

所有这些方法都吃随机数。为了让实验可复现（reproducible），要固定随机种子（seed）：同一个 seed 喂给伪随机数发生器，吐出的随机序列完全一样。本课代码统一用 np.random.default_rng(seed) 建一个独立的发生器——比全局 np.random.seed 更干净、不会被别处偷偷改状态。

调一调，观察现象

下面每个小实验都只用 numpy + print，几秒跑完。先猜，再跑，对不上就回头想机制。

import numpy as np
rng = np.random.default_rng(0)
lam = 1.5
for N in [100, 400, 1600, 6400]:
    ests = []
    for _ in range(2000):                 # 重复多次让 std 估计更稳
        U = rng.uniform(0, 1, N)
        X = -np.log(1 - U) / lam          # 逆变换采样
        ests.append(X.mean())
    print(f"N={N:>5}  std(estimate)={np.std(ests):.4f}  "
          f"1/sqrtN*sigma={(1/lam)/np.sqrt(N):.4f}")

import numpy as np
rng = np.random.default_rng(0)
lam = 1.5
U = rng.uniform(0, 1, 200000)
X = -np.log(1 - U) / lam
edges = np.linspace(0, 5, 11)
counts, _ = np.histogram(X, bins=edges, density=True)
centers = 0.5 * (edges[:-1] + edges[1:])
theory = lam * np.exp(-lam * centers)
print("x      :", np.round(centers, 2))
print("empir. :", np.round(counts, 3))
print("theory :", np.round(theory, 3))

import numpy as np
rng = np.random.default_rng(1)
eps = rng.standard_normal(5)
print("fixed eps:", np.round(eps, 3))
sigma = 1.0
for mu in [0.0, 1.0, 2.0]:
    print(f"mu={mu}: x =", np.round(mu + sigma * eps, 3))

# 同一批大 eps 上做数值微分，验证 dE[x]/dmu=1, dE[x]/dsigma=mean(eps)
eps = rng.standard_normal(100000)
mu, sigma, h = 2.0, 0.5, 1e-3
base = (mu + sigma * eps).mean()
print("dE/dmu    ~", round(((mu + h + sigma * eps).mean() - base) / h, 4))
print("dE/dsigma ~", round(((mu + (sigma + h) * eps).mean() - base) / h, 4),
      " (mean(eps)=", round(eps.mean(), 4), ")")   # 两数应相等

动手练习

1. 逆变换造别的分布。 标准 Logistic 分布的 CDF 是 \(F(x)=\frac{1}{1+e^{-x}}\)（就是 sigmoid！）。推导 \(F^{-1}(U)\)，并用下面骨架采样、对比经验均值≈0。

   import numpy as np
   rng = np.random.default_rng(0)
   U = rng.uniform(0, 1, 100000)
   X = np.log(U / (1 - U))   # F^{-1}(U) = logit(U)，自己推一遍
   print("mean:", round(X.mean(), 3), " median:", round(np.median(X), 3))
   

2. 重要性采样动手。 估计 \(\mathbb E_p[X^2]\)，其中 \(p=\mathcal N(3,1)\)（真值 \(=\mu^2+\sigma^2=10\)），但只从 \(q=\mathcal N(0,2^2)\) 采样。补全权重 \(p/q\)。

   import numpy as np
   rng = np.random.default_rng(1)
   p = lambda x: np.exp(-0.5*(x-3)**2)/np.sqrt(2*np.pi)
   q = lambda x: np.exp(-0.5*(x/2)**2)/(2*np.sqrt(2*np.pi))
   xs = rng.normal(0, 2, 200000)
   w  = p(xs) / q(xs)                 # 重要性权重
   print("IS est:", round(np.mean(w * xs**2), 3), " (theory 10)")
   

   现在把 \(q\) 改成 \(\mathcal N(0,0.5^2)\)（尾巴太薄），观察估计变得多么不稳——体会"权重爆炸"。
3. 拒绝采样最小实现。 从三角形密度 \(p(x)=2x,\ x\in[0,1]\) 采样：用 \(q=\mathrm{Uniform}(0,1)\)、罩子 \(M=2\)。提议 \(x\sim U(0,1)\)，再抽 \(u\sim U(0,1)\)，当 \(u
4. 验证 \(1/\sqrt N\) 的斜率。 把调一调 1 的 \((\,\ln N,\ \ln\mathrm{std}\,)\) 打印出来，估计斜率，应接近 \(-0.5\)。
5. 种子可复现。 用相同 seed 跑两次同一段采样，验证两次结果逐位相同；换 seed 则不同。

掌握自检

• 为什么蒙特卡洛误差是 \(1/\sqrt N\) 而不是 \(1/N\)？这和"随机游走 \(N\) 步约走 \(\sqrt N\)"是同一回事吗？
• 逆变换采样的一行证明里，哪一步用到了"\(U\) 均匀 \(\Rightarrow P(U\le t)=t\)"？
• 为什么"直接从 \(\mathcal N(\mu,\sigma^2)\) 采样"会切断梯度，而 \(x=\mu+\sigma\varepsilon\) 不会？随机性被挪到哪了？参数从"权重"挪到了"被积函数"，这为何让求导能与求期望交换？
• 不查公式，写出 \(\partial x/\partial\mu\) 和 \(\partial x/\partial\sigma\)。
• 重要性采样里权重 \(p/q\) 在"\(q\) 采多了的地方"是大于还是小于 1？为什么 \(q\) 尾巴要比 \(p\) 厚？
• Metropolis–Hastings 的接受概率 \(\min(1,p_{\text{新}}/p_{\text{旧}})\) 为什么只用到 \(p\) 的比值？这让我们能对哪类分布（少知道什么）也能采样？
• LLM 的 temperature 调大，输出更随机还是更确定？对应分布变尖还是变平？

下一课我们正式进入最优化：凸性与梯度下降的几何——为什么沿负梯度走能下山，Hessian \(H\) 的特征值/条件数（呼应模块1）又如何决定下山的快慢与形状。本课的采样会和它在后面的随机优化、变分推断里重逢。