Part I · 第 4 章概率 · 4.9

采样方法

计算机只会掷一种骰子：[0,1) 上的均匀随机数。可机器人要的是"带噪声的里程计会把我带到哪"、"这个读数有多可能"、"一个没法积分的后验的期望"。采样就是把均匀随机数变成任意分布的样本，再用样本的平均代替积分。本节讲四种思路：逆 CDF、拒绝采样、重要性采样、MCMC，最后把它们用到机器人运动模型上，为 4.10 节的粒子滤波打下基础。

约 25 分钟前置：4.1 随机变量与分布 · 4.4 高斯分布 · 4.7 贝叶斯滤波 5 个交互

1为什么机器人需要它

有三类事机器人经常要做，却没法用公式一步算完：

模拟不确定性：里程计每一步都带噪声，走 5 米以后机器人可能在哪？转弯多了，分布不再是椭圆，而是一根"香蕉"。最诚实的办法是让成百上千个"可能的机器人"各自带着随机噪声走一遍，看它们落在哪。
近似期望与积分：规划里要算"平均代价"，学习里要算"期望回报"，都是 $\mathbb E[f(x)]=\int f(x)p(x)\,dx$。积分算不出来时，抽 $N$ 个样本取平均，误差大约按 $1/\sqrt N$ 减小，而且与维度无关。
从复杂分布里抽样：后验往往只知道一个"差个常数"的形式（贝叶斯公式的分母没法算）。能从这样的分布里抽样的算法，是现代概率推断的发动机，也是粒子滤波 (4.10)、随机优化 (6.11) 和强化学习里的探索与重要性权重的底层工具。

程序库已经给了我们均匀随机数（4.8 节里的 randn 也由它变来）。本节要回答的，就是怎样从它出发，得到想要的分布。

2直觉：把均匀随机数"挤"成想要的形状

画出分布的累积分布函数 (CDF) $F(x)=P(X\le x)$，它从 0 单调爬到 1。从纵轴 $[0,1]$ 上均匀地选一个高度 $u$，水平走到曲线，再垂直落到横轴，落点 $x=F^{-1}(u)$ 就是一个样本。

为什么对？CDF 陡的地方是概率密度高的地方；纵轴上等长的一段，在陡的地方只对应横轴上很窄的一段。也就是说，同样多的"雨点"落在纵轴上，被挤进了密度高的区域：样本就按照分布的形状堆积起来了。

FIG 1逆 CDF 采样

上方：蓝色是 CDF $F(x)$。拖左侧的圆点选一个 $u$，紫色线水平走到曲线、再垂直落下得到 $x$。"十分位"按钮画出 $u=0.05,0.15,\dots,0.95$ 的十条线：纵轴上等距，横轴上却在高密度处挤在一起。下方：黄线是概率密度 $p(x)$，柱子是已抽样本的直方图（密度归一）。

试一试

拖动 $u$，观察 $x$ 怎么变：在 CDF 陡的地方，$u$ 动一大段，$x$ 动多少？
选"双峰"分布，抽取至少 500 个样本，让直方图贴近黄色曲线。

核心洞见

逆 CDF 采样是"精确"的：$F^{-1}(U)$ 的分布恰好是 $p$。可惜它要求你能算（或查表）$F^{-1}$。高斯没有解析的反函数，所以标准库用 Box–Muller 之类的变换；高维、复杂的分布就更没有办法了。下面三种方法放宽了这个要求。

3动手玩：拒绝采样，在"盒子"里撒点

目标分布 $\tilde p(x)$（只要知道形状，不必归一化）没法直接抽。找一个容易抽的提议分布 $q(x)$（均匀、高斯都行），再选常数 $M$，使得曲线 $Mq(x)$ 在每一点都盖住 $\tilde p(x)$。

做法是撒飞镖：先从 $q$ 抽一个横坐标 $x$，再在 $[0,\,Mq(x)]$ 上均匀抽一个高度 $y$。飞镖落在 $\tilde p$ 曲线以下就接受，否则拒绝。被接受的点的 $x$ 就服从 $\tilde p$。包络越贴近目标，浪费越少。

FIG 2拒绝采样：包络、接受率与"盖不住"的后果

飞镖：绿 = 接受，红 = 拒绝

被接受样本的直方图 vs 真实目标

黄色是目标 $\tilde p$（两个凸起），紫色虚线是包络 $Mq(x)$。把 $M$ 调小到包络压到目标下面，会发现直方图的峰被"削平"了，因为那些本该被接受的高飞镖根本没有机会出现。包络不盖住目标时，样本是错的，而且没有任何警告。

试一试

省飞镖：让 $M$ 恰好盖住目标（不低于最小值），同时接受率超过 40%。需要换哪个提议分布？
制造偏差：把 $M$ 调到最小有效值的 70% 以下，抽至少 300 个飞镖，观察被削平的峰。

4第二个视角：不丢样本，给它们加权

拒绝采样会扔掉大部分样本。重要性采样 (importance sampling) 一个都不扔：直接从提议分布 $q$ 抽样，再给每个样本一个权重 $w=p(x)/q(x)$，补偿"抽样的地方和真正关心的地方不一致"。想求 $\mathbb E_p[f]$，只需做加权平均。

关键是 $q$ 要选得好。如果 $q$ 几乎碰不到 $p$ 的重要区域，少数几个样本的权重会巨大，其余几乎为 0，估计就由寥寥几个样本决定，噪声很大。衡量这件事的指标是有效样本数 $\mathrm{ESS}=(\sum w)^2/\sum w^2$，它告诉你"这 $N$ 个加权样本，相当于几个理想样本"。

FIG 3重要性采样：提议分布好不好，看 ESS

目标 p、提议 q 与带权样本

同样设置重复 40 次的估计值

黄 = 目标 $p$（双峰），青 = 提议 $q$（一个高斯：拖山顶改均值，拖腰部改宽度）。底部紫色圆点是 100 个样本，圆点越大权重越大。下方图把"同一设置独立重复 40 次"得到的估计值画出来，绿线是真值：点越散，估计越不可靠。

试一试

好的提议：估计均值时，让 ESS/N 超过 0.6。
坏的提议：把 $q$ 拖到离 $p$ 很远或收得很窄，让 ESS/N 低于 0.1。
稀有事件：切到「尾部概率 P(x>8)」，找一个 $q$，让 40 次估计的标准差小于 0.024，且平均值与真值相差不到 0.01。（提示：只把 $q$ 停在右侧尾部，会得到一个很稳定、但完全错误的答案，因为左边那个山头没有被覆盖。）

5第三个视角：让样本"走"起来，MCMC

高维时，上面两种方法都会失效：找不到一个既容易抽、又贴着目标的 $q$。马尔可夫链蒙特卡洛 (MCMC) 换了思路：不再一次独立地抽一个，而是让一个点在空间里随机游走，游走的规则保证长期来看它停留在每个区域的时间，正比于目标的概率。

最简单的 Metropolis 算法：当前在 $x$，随机提议一个邻居 $x'=x+\text{噪声}$。若 $\tilde p(x')\ge\tilde p(x)$（爬坡）就一定接受；否则以概率 $\tilde p(x')/\tilde p(x)$ 接受（偶尔下坡），被拒绝就原地不动、再记一次。注意：只用到比值，所以 $\tilde p$ 不需要归一化。

唯一要你调的是步长。太小：几乎每步都被接受，但走得像蜗牛，样本高度相关，很难走到第二个山头。太大：步子一迈就掉到低密度区，几乎全被拒绝，原地踏步。下面的目标有两个山头。

FIG 4Metropolis 随机游走：步长决定一切

背景的青色深浅是目标密度（两个山头）。黄色圆圈是起点，可拖动（改变后需"重置"）。紫色折线是最近 60 步，红色小叉是最近被拒绝的提议，淡点是最近 2000 个样本。"ESS"是考虑自相关之后相当于多少个独立样本。

试一试

蜗牛：步长 ≤ 0.2，跑到 2000 步以上：几乎 100% 接受，却只待在一个山头里。
原地踏步：步长 ≥ 3，跑 500 步以上：接受率低于 10%。
刚刚好：找一个步长，接受率在 20%–60%，两个山头各占至少 25% 的样本（至少 3000 步）。

接受率不是越高越好

接受率接近 100% 通常说明步子太小，而不是采样器很棒。经验上，高维高斯的最优随机游走接受率约 23%；一维、二维可略高。真正的质量指标是 ESS：同样的 CPU 时间，能得到多少个"几乎独立"的样本。

6再看数学

提示：悬停下面公式中带底色的符号，对应的图中曲线会高亮（逆 CDF 对应 FIG 1，包络对应 FIG 2，权重对应 FIG 3）。

蒙特卡洛：用平均代替积分

$$ \mathbb E_p[f(x)] = \int f(x)\,p(x)\,dx \;\approx\; \frac1N\sum_{i=1}^N f(x_i),\quad x_i\overset{\text{iid}}{\sim} p,\qquad \text{std. error} = \frac{\sigma_f}{\sqrt N}. $$

误差按 $1/\sqrt N$ 下降：想让误差缩小到 1/10，要多抽 100 倍样本。好处是这个速度不依赖维度，而网格积分的代价随维度指数增长。

逆变换采样

若 $U\sim\mathrm{Uniform}(0,1)$，则 $X=\hl{icdf}{F^{-1}(U)}$ 的分布函数恰为 $F$：

$$ P(X\le x) = P\big(F^{-1}(U)\le x\big) = P\big(U\le F(x)\big) = F(x). $$

例：指数分布 $F(x)=1-e^{-\lambda x}$，求逆得 $x=-\ln(1-u)/\lambda$。高斯没有解析反函数，标准库用 Box–Muller：取两个独立均匀数 $u_1,u_2$，则 $z=\sqrt{-2\ln u_1}\,\cos(2\pi u_2)$ 就是标准正态样本（本站 R101.M.randn 用的就是它）。再令 $x=\mu+\sigma z$ 即得 $\mathcal N(\mu,\sigma^2)$。多维高斯 $\mathcal N(\boldsymbol\mu,\Sigma)$：取 $\Sigma=LL^\top$（Cholesky 分解），$\mathbf x=\boldsymbol\mu+L\mathbf z$，这正是 4.5 节线性变换传播协方差的反向使用。

拒绝采样

$$ x\sim q,\quad u\sim\mathrm{U}(0,1),\quad \text{accept if } u\le \frac{\tilde p(x)}{\hl{env}{Mq(x)}},\qquad \tilde p(x)\le Mq(x)\;\;\forall x. $$

被接受的 $x$ 的密度正比于 $q(x)\cdot\tilde p(x)/(Mq(x))\propto\tilde p(x)$。接受率 $=Z/M$，其中 $Z=\int\tilde p\,dx$ 是目标的"面积"，也就是目标下的面积占包络下面积的比例，和 FIG 2 的"盒子"直观一致。在 $d$ 维，让包络贴紧目标越来越难，接受率常常随维度指数衰减。

重要性采样

$$ \mathbb E_p[f] = \int f(x)\frac{p(x)}{q(x)}\,q(x)\,dx \approx \frac{\sum_i \hl{wts}{w_i}\,f(x_i)}{\sum_i w_i},\qquad x_i\sim q,\;\; w_i=\frac{\tilde p(x_i)}{q(x_i)}, $$ $$ \mathrm{ESS}=\frac{(\sum_i w_i)^2}{\sum_i w_i^2}. $$

用分母 $\sum w_i$ 归一化（"自归一化"）后，$p$ 也只需要已知到相差一个常数。$\mathrm{ESS}\in[1,N]$：权重全相等时等于 $N$，一个样本独占全部权重时等于 1。要求 $q$ 在 $p$ 有质量的地方都不为 0，而且最好尾巴比 $p$ 更重，否则偶尔出现的巨大权重会让方差无穷大。

Metropolis–Hastings

$$ x'\sim g(x'\mid x),\qquad \alpha = \min\Big(1,\;\frac{\tilde p(x')\,g(x\mid x')}{\tilde p(x)\,g(x'\mid x)}\Big)\;\overset{g\text{ symmetric}}{=}\;\min\Big(1,\frac{\tilde p(x')}{\tilde p(x)}\Big). $$

这样构造的链满足细致平衡 $\tilde p(x)\,T(x\to x')=\tilde p(x')\,T(x'\to x)$，所以 $p$ 是它的平稳分布；链收敛以后，样本就近似服从 $p$。开头那一段还没收敛的样本叫"预热 (burn-in)"，通常丢掉。样本之间有自相关 $\rho_k$，有效样本数约为 $N/(1+2\sum_{k\ge1}\rho_k)$。

方法	需要什么	优点	弱点
逆 CDF	$F^{-1}$	精确、独立、不浪费	只适用于一维或结构特别的分布
拒绝采样	$\tilde p,\,q,\,M$	精确、独立、$\tilde p$ 不必归一化	高维时接受率崩溃
重要性采样	$\tilde p,\,q$	不丢样本；可估稀有事件；粒子滤波的基础	$q$ 差则权重退化 (ESS 低)
MCMC	$\tilde p,\,g$	高维也能用；只要比值	样本相关、要预热、难判断是否收敛

7机器人中的它：从运动模型采样出"香蕉"

机器人沿一条圆弧前进，每一小步的航向有随机偏差，步长也有一点比例误差。把"一步"的噪声模型反复套用 20 次，就得到走完全程的位姿分布，却写不出它的解析形式。最简单的办法：采样。让 400 个"可能的机器人"各自带着独立噪声走一遍。

这正是粒子滤波预测步骤里的 sample_motion_model (Thrun 等《概率机器人》第 5 章)。同时我们也画出 4.5 节的做法：沿名义路径用雅可比线性化、逐步传播协方差 $F\Sigma F^\top+Q$，得到的 2σ 椭圆。航向噪声小时椭圆很准；航向噪声大、走得远时，样本弯成一根香蕉，任何椭圆都罩不住。

FIG 5采样运动模型 vs 线性化高斯

黄点是 400 个采样位姿，青色虚线是名义（无噪声）路径，紫色椭圆是沿名义路径线性化传播得到的 2σ 椭圆（4.5 节）。若分布真是高斯，约 86.5% 的点应落在 2σ 椭圆内；读数里的"椭圆内比例"偏离这个值越多，高斯假设越差。

试一试

把航向噪声调到 3°/√m 以下：样本云几乎是个椭圆，高斯近似很好（椭圆内比例约 80%–90%）。
把航向噪声调大（≥ 12°/√m）、距离拉到 7 m 以上：样本弯成"香蕉"，椭圆内比例掉到 78% 以下，而且样本均值偏离了名义终点（椭圆始终以名义终点为中心）。

真实系统里

"从一个你写得出前向仿真、却写不出密度的模型里采样"是机器人里最常用的采样用法：粒子滤波的预测步 (4.10)；域随机化里在仿真中采样质量、摩擦、延迟，训练出对现实更鲁棒的策略；MPPI 与 CEM 这类采样式模型预测控制，每个控制周期采样成百上千条带噪声的轨迹：MPPI 按 $\exp(-\text{代价}/\lambda)$ 给它们加权，这正是重要性采样里的权重，CEM 则只保留代价最低的一批来重新拟合采样分布；RRT 在构型空间里随机采样；蒙特卡洛树搜索靠随机展开估计价值。

8常见误区

误区 1：包络没盖住目标

FIG 2 已经演示：$Mq(x)<\tilde p(x)$ 的地方，样本分布被悄悄截断，而程序照样输出看起来正常的数字。拒绝采样要在全域上验证 $\tilde p/q$ 的上界，不能只在你看到的几个点上试。

误区 2：重要性权重退化

$q$ 和 $p$ 差得远时，ESS 会塌到 1 左右：1000 个样本只剩一个在起作用，估计值看似"收敛"，其实只是那一个样本。务必同时报告 ESS。粒子滤波里这叫权重退化，重采样就是为它准备的（4.10 节）。

误区 3：把 MCMC 的相邻样本当作独立样本

相邻样本高度相关，10000 步的 ESS 可能只有几十。用它们算误差线会严重低估不确定性。还要警惕"看起来收敛了"：蜗牛步长的链安安稳稳待在一个山头里，任何单链的诊断图都显示一切正常，另一个山头却从未被访问。办法：多条不同起点的链，比较它们的结果。

误区 4：随机数用得不规范

机器人系统里要固定并记录随机种子，否则你无法复现一次失败的仿真；并行线程不要共享同一个全局生成器，也不要给每个线程用同一个种子，否则所有"独立"的粒子其实在走同一条路。本课网页里用了带种子的生成器，所以每次打开同一个图，数据都相同。

9小练习

1. 指数分布 $F(x)=1-e^{-\lambda x}$，$\lambda=2$。用逆 CDF 法，均匀数 $u=0.5$ 对应的样本是？

$x=-\ln(1-u)/\lambda=-\ln(0.5)/2=0.693/2\approx0.347$。它正好是这个分布的中位数，因为 $u=0.5$ 就是 CDF 的 50% 处。

2. 目标曲线 $\tilde p$ 下的面积是 2；提议是 $[0,10]$ 上的均匀分布，包络是高度 0.5 的水平线（$\tilde p\le 0.5$ 处处成立）。接受率是多少？

包络下的面积是 $0.5\times10=5$，目标下的面积是 2，飞镖落在目标下的比例 $=2/5=40\%$。

3. 四个重要性权重（已归一化）为 $(0.7,\,0.1,\,0.1,\,0.1)$。有效样本数 ESS 约为？

权重和为 1，$\sum w^2=0.49+0.01\times3=0.52$，$\mathrm{ESS}=1/0.52\approx1.92$。4 个样本只相当于不到 2 个理想样本，因为 0.7 的那个独占了大部分权重。

4. 一条 Metropolis 链的接受率是 99%。最可能的情况是？

步子太小时，$x'$ 与 $x$ 的密度几乎相同，比值接近 1，当然几乎全被接受，可是走得慢。（目标恰好是均匀分布也会得到 100%，但那是个特例。）要的是"又远又被接受"的折中，一般让接受率落在 20%–50%。

5. 用 $N=100$ 个独立样本做蒙特卡洛估计，标准误差约 0.1。改用 $N=10{,}000$，标准误差约为？

误差 $\propto1/\sqrt N$。样本多 100 倍，误差只缩小 $\sqrt{100}=10$ 倍，得 0.01。想再缩小 10 倍，又要多 100 倍的样本，所以蒙特卡洛适合"够用就行"的精度，不适合高精度。

10延伸阅读

C. M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006, Ch. 11 — 采样方法：逆变换、拒绝、重要性、Metropolis–Hastings 与 Gibbs
D. J. C. MacKay, Information Theory, Inference, and Learning Algorithms, Cambridge University Press, 2003, Ch. 29 — 蒙特卡洛方法，讲得非常直观，附免费电子版
S. Thrun, W. Burgard & D. Fox, Probabilistic Robotics, MIT Press, 2005, Ch. 5 & 4.3 — 从运动模型中采样 (sample_motion_model) 与粒子滤波
C. P. Robert & G. Casella, Monte Carlo Statistical Methods, Springer, 2004 — 系统而严谨的蒙特卡洛参考书
N. Metropolis et al., "Equation of State Calculations by Fast Computing Machines", Journal of Chemical Physics, 1953; W. K. Hastings, Biometrika, 1970 — Metropolis 算法的原始论文与推广

方法	需要什么	优点	弱点
逆 CDF	\(F^{-1}\)	精确、独立、不浪费	只适用于一维或结构特别的分布
拒绝采样	\(\tilde p,\,q,\,M\)	精确、独立、\(\tilde p\) 不必归一化	高维时接受率崩溃
重要性采样	\(\tilde p,\,q\)	不丢样本；可估稀有事件；粒子滤波的基础	\(q\) 差则权重退化 (ESS 低)
MCMC	\(\tilde p,\,g\)	高维也能用；只要比值	样本相关、要预热、难判断是否收敛