Part I · 第 4 章 概率 · 4.9Part I · Ch. 4 Probability · 4.9

采样方法Sampling Methods

计算机只会掷一种骰子:[0,1) 上的均匀随机数。可机器人要的是"带噪声的里程计会把我带到哪"、"这个读数有多可能"、"一个没法积分的后验的期望"。采样就是把均匀随机数变成任意分布的样本,再用样本的平均代替积分。本节讲四种思路:逆 CDF、拒绝采样、重要性采样、MCMC,最后把它们用到机器人运动模型上,为 4.10 节的粒子滤波打下基础。 A computer rolls only one kind of die: a uniform random number on [0,1). Robots, however, need "where will noisy odometry take me?", "how likely is this reading?" and "the expectation under a posterior that cannot be integrated". Sampling turns uniform random numbers into samples from any distribution and replaces integrals by sample averages. This lesson covers four ideas, the inverse CDF, rejection sampling, importance sampling and MCMC, and then applies them to a robot's motion model, preparing the ground for the particle filter of Lesson 4.10.

约 25 分钟~25 min 前置:4.1 随机变量与分布 · 4.4 高斯分布 · 4.7 贝叶斯滤波Prereqs: 4.1 Random Variables · 4.4 Gaussians · 4.7 Bayes Filter 5 个交互5 interactives

    1为什么机器人需要它Why robots need this

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

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

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

    Robots keep running into three jobs that no formula finishes in one step:

    • Simulating uncertainty: odometry is noisy at every step, so where might the robot be after 5 m? After a few turns the distribution is no longer an ellipse but a "banana". The most honest approach is to let hundreds of "possible robots" each walk with their own random noise and see where they land.
    • Approximating expectations and integrals: planning needs an "average cost", learning an "expected return", both of the form \(\mathbb E[f(x)]=\int f(x)p(x)\,dx\). When the integral cannot be done, draw \(N\) samples and average; the error shrinks like \(1/\sqrt N\), independently of the dimension.
    • Drawing from complicated distributions: a posterior is often known only "up to a constant" (the denominator of Bayes' rule is intractable). Algorithms that sample from such distributions are the engine of modern probabilistic inference, and underlie the particle filter (4.10), sampling-based optimization (6.11) and the exploration and importance weights of reinforcement learning.

    The library already gives us uniform random numbers (the randn used in Lesson 4.8 is built from them). This lesson answers how to get from there to the distribution we want.

    2直觉:把均匀随机数"挤"成想要的形状Intuition: squeeze uniform numbers into the shape you want

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

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

    Draw the distribution's cumulative distribution function (CDF) \(F(x)=P(X\le x)\), rising monotonically from 0 to 1. Pick a height \(u\) uniformly on the vertical axis \([0,1]\), walk horizontally to the curve, then drop vertically to the horizontal axis. The landing point \(x=F^{-1}(u)\) is a sample.

    Why does it work? Where the CDF is steep the density is high; an equal-length piece of the vertical axis covers only a narrow piece of the horizontal axis there. The same number of "raindrops" on the vertical axis get squeezed into high-density regions, so the samples pile up in the shape of the distribution.

    FIG 1逆 CDF 采样Inverse-CDF sampling
    上方:蓝色是 CDF \(F(x)\)。拖左侧的圆点选一个 \(u\),紫色线水平走到曲线、再垂直落下得到 \(x\)。"十分位"按钮画出 \(u=0.05,0.15,\dots,0.95\) 的十条线:纵轴上等距,横轴上却在高密度处挤在一起。下方:黄线是概率密度 \(p(x)\),柱子是已抽样本的直方图(密度归一)。 Top: the blue curve is the CDF \(F(x)\). Drag the dot on the left to choose \(u\); the purple line walks horizontally to the curve and then drops to give \(x\). The "Deciles" button draws ten lines at \(u=0.05,0.15,\dots,0.95\): evenly spaced on the vertical axis, yet bunched together where the density is high. Bottom: the amber line is the density \(p(x)\); the bars are the histogram of the samples drawn so far (density-normalized).
    试一试Try this
    • 拖动 \(u\),观察 \(x\) 怎么变:在 CDF 陡的地方,\(u\) 动一大段,\(x\) 动多少?Drag \(u\) and watch \(x\): where the CDF is steep, \(u\) moves a lot while \(x\) moves how much?
    • 选"双峰"分布,抽取至少 500 个样本,让直方图贴近黄色曲线。Choose the "bimodal" distribution and draw at least 500 samples so the histogram hugs the amber curve.
    核心洞见Key insight

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

    Inverse-CDF sampling is "exact": the distribution of \(F^{-1}(U)\) is exactly \(p\). The catch is that you must be able to compute (or tabulate) \(F^{-1}\). The Gaussian has no closed-form inverse, so libraries use transforms such as Box–Muller; for high-dimensional or complicated distributions there is no hope. The next three methods relax this requirement.

    3动手玩:拒绝采样,在"盒子"里撒点Play: rejection sampling, throwing darts into a box

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

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

    The target \(\tilde p(x)\) (only its shape is needed, no normalization) cannot be sampled directly. Find an easy proposal distribution \(q(x)\) (uniform or Gaussian will do) and a constant \(M\) such that the curve \(Mq(x)\) covers \(\tilde p(x)\) everywhere.

    Then throw darts: draw a horizontal position \(x\) from \(q\), then a height \(y\) uniformly in \([0,\,Mq(x)]\). If the dart lands below the \(\tilde p\) curve accept it, otherwise reject. The \(x\) of the accepted darts follow \(\tilde p\). The closer the envelope hugs the target, the less is wasted.

    FIG 2拒绝采样:包络、接受率与"盖不住"的后果Rejection sampling: envelope, acceptance rate and what happens when it fails to cover
    飞镖:绿 = 接受,红 = 拒绝DARTS: GREEN = ACCEPTED, RED = REJECTED
    被接受样本的直方图 vs 真实目标HISTOGRAM OF ACCEPTED SAMPLES VS THE TRUE TARGET
    黄色是目标 \(\tilde p\)(两个凸起),紫色虚线是包络 \(Mq(x)\)。把 \(M\) 调小到包络压到目标下面,会发现直方图的峰被"削平"了,因为那些本该被接受的高飞镖根本没有机会出现。包络不盖住目标时,样本是错的,而且没有任何警告。 The amber curve is the target \(\tilde p\) (two bumps), the purple dashed curve the envelope \(Mq(x)\). Lower \(M\) until the envelope dips below the target and the histogram's peak is "shaved off", because the high darts that should have been accepted never get the chance to exist. When the envelope fails to cover the target, the samples are wrong, with no warning at all.
    试一试Try this
    • 省飞镖:让 \(M\) 恰好盖住目标(不低于最小值),同时接受率超过 40%。需要换哪个提议分布?Save darts: make \(M\) just cover the target (not below the minimum) while the acceptance rate exceeds 40%. Which proposal do you need?
    • 制造偏差:把 \(M\) 调到最小有效值的 70% 以下,抽至少 300 个飞镖,观察被削平的峰。Create bias: set \(M\) below 70% of the minimum valid value and throw at least 300 darts to see the shaved peak.

    4第二个视角:不丢样本,给它们加权A second view: keep every sample, but weight it

    拒绝采样会扔掉大部分样本。重要性采样 (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\) 个加权样本,相当于几个理想样本"。

    Rejection sampling throws most samples away. Importance sampling throws none away: it draws straight from the proposal \(q\), then gives each sample a weight \(w=p(x)/q(x)\) that compensates for "sampling where we do not really care". To estimate \(\mathbb E_p[f]\), just take a weighted average.

    The key is choosing a good \(q\). If \(q\) barely visits the important region of \(p\), a few samples get huge weights and the rest almost none, so the estimate hangs on a handful of samples and is very noisy. The measure of this is the effective sample size \(\mathrm{ESS}=(\sum w)^2/\sum w^2\), which says "these \(N\) weighted samples are worth this many ideal samples".

    FIG 3重要性采样:提议分布好不好,看 ESSImportance sampling: is the proposal any good? Look at ESS
    目标 p、提议 q 与带权样本TARGET p, PROPOSAL q AND WEIGHTED SAMPLES
    同样设置重复 40 次的估计值ESTIMATES FROM 40 REPEATS OF THE SAME SETUP
    黄 = 目标 \(p\)(双峰),青 = 提议 \(q\)(一个高斯:拖山顶改均值,拖腰部改宽度)。底部紫色圆点是 100 个样本,圆点越大权重越大。下方图把"同一设置独立重复 40 次"得到的估计值画出来,绿线是真值:点越散,估计越不可靠。 Amber = target \(p\) (bimodal), cyan = proposal \(q\) (one Gaussian: drag the peak for its mean, the flank for its width). The purple dots at the bottom are 100 samples; a bigger dot means a bigger weight. The lower plot shows the estimates from 40 independent repeats of the same setup; the green line is the truth. The more the dots scatter, the less reliable the estimate.
    试一试Try this
    • 好的提议:估计均值时,让 ESS/N 超过 0.6。A good proposal: when estimating the mean, get ESS/N above 0.6.
    • 坏的提议:把 \(q\) 拖到离 \(p\) 很远或收得很窄,让 ESS/N 低于 0.1。A bad proposal: drag \(q\) far from \(p\) or squeeze it narrow so ESS/N drops below 0.1.
    • 稀有事件:切到「尾部概率 P(x>8)」,找一个 \(q\),让 40 次估计的标准差小于 0.024,且平均值与真值相差不到 0.01。(提示:只把 \(q\) 停在右侧尾部,会得到一个很稳定、但完全错误的答案,因为左边那个山头没有被覆盖。)A rare event: switch to the tail probability P(x>8) and find a \(q\) such that the 40 estimates have a standard deviation below 0.024 and an average within 0.01 of the truth. (Hint: parking \(q\) on the right tail alone gives a very stable but wrong answer, because the left hill is left uncovered.)

    5第三个视角:让样本"走"起来,MCMCA third view: let the sample walk, MCMC

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

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

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

    In high dimensions both methods above break down: no \(q\) is both easy to sample and tight around the target. Markov chain Monte Carlo (MCMC) takes another route: instead of drawing independent samples one at a time, let a point wander randomly through space under rules that guarantee that, in the long run, the time it spends in each region is proportional to the target's probability.

    The simplest is the Metropolis algorithm: standing at \(x\), propose a neighbour \(x'=x+\text{noise}\). If \(\tilde p(x')\ge\tilde p(x)\) (uphill), always accept; otherwise accept with probability \(\tilde p(x')/\tilde p(x)\) (an occasional step downhill); if rejected, stay and count \(x\) again. Note that only a ratio is used, so \(\tilde p\) need not be normalized.

    The one thing you must tune is the step size. Too small: almost every step is accepted but the walk crawls, samples are highly correlated, and the second hill is hardly ever reached. Too large: each step lands in a low-density region and nearly all are rejected, so the chain stands still. The target below has two hills.

    FIG 4Metropolis 随机游走:步长决定一切Random-walk Metropolis: the step size decides everything
    背景的青色深浅是目标密度(两个山头)。黄色圆圈是起点,可拖动(改变后需"重置")。紫色折线是最近 60 步,红色小叉是最近被拒绝的提议,淡点是最近 2000 个样本。"ESS"是考虑自相关之后相当于多少个独立样本。 The cyan shading in the background is the target density (two hills). The amber ring is the starting point (draggable; press "Reset" afterwards). The purple path shows the last 60 steps, red crosses are recently rejected proposals, faint dots are the last 2000 samples. "ESS" is how many independent samples the chain is worth after accounting for autocorrelation.
    试一试Try this
    • 蜗牛:步长 ≤ 0.2,跑到 2000 步以上:几乎 100% 接受,却只待在一个山头里。The snail: step size ≤ 0.2, run 2000+ steps: nearly 100% accepted, yet stuck on a single hill.
    • 原地踏步:步长 ≥ 3,跑 500 步以上:接受率低于 10%。Standing still: step size ≥ 3, run 500+ steps: acceptance below 10%.
    • 刚刚好:找一个步长,接受率在 20%–60%,两个山头各占至少 25% 的样本(至少 3000 步)。Just right: find a step size with acceptance 20%–60% and at least 25% of the samples on each hill (3000+ steps).
    接受率不是越高越好A higher acceptance rate is not better

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

    An acceptance rate near 100% usually means the steps are too small, not that the sampler is excellent. As a rule of thumb the optimal random-walk acceptance rate for high-dimensional Gaussians is around 23%, somewhat higher in one or two dimensions. The real quality measure is ESS: how many nearly independent samples you get for the same CPU time.

    6再看数学The math

    提示:悬停下面公式中带底色的符号,对应的图中曲线会高亮(逆 CDF 对应 FIG 1,包络对应 FIG 2,权重对应 FIG 3)。Tip: hover the shaded symbols below to light up the matching curve (inverse CDF → FIG 1, envelope → FIG 2, weights → FIG 3).

    蒙特卡洛:用平均代替积分Monte Carlo: averages instead of integrals

    $$ \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 倍样本。好处是这个速度不依赖维度,而网格积分的代价随维度指数增长。

    The error falls like \(1/\sqrt N\): to cut the error by 10 you need 100 times more samples. The good news is that this rate is independent of the dimension, whereas grid integration grows exponentially with dimension.

    逆变换采样Inverse transform sampling

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

    If \(U\sim\mathrm{Uniform}(0,1)\), then \(X=\hl{icdf}{F^{-1}(U)}\) has distribution function exactly \(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 节线性变换传播协方差的反向使用。

    Example: the exponential has \(F(x)=1-e^{-\lambda x}\), inverted to \(x=-\ln(1-u)/\lambda\). The Gaussian has no closed-form inverse, so libraries use Box–Muller: take two independent uniforms \(u_1,u_2\), and \(z=\sqrt{-2\ln u_1}\,\cos(2\pi u_2)\) is a standard normal sample (this is what R101.M.randn on this site does). Then \(x=\mu+\sigma z\) gives \(\mathcal N(\mu,\sigma^2)\). For a multivariate \(\mathcal N(\boldsymbol\mu,\Sigma)\): factor \(\Sigma=LL^\top\) (Cholesky) and set \(\mathbf x=\boldsymbol\mu+L\mathbf z\), which is Lesson 4.5's rule for transforming covariances, used in reverse.

    拒绝采样Rejection sampling

    $$ 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\) 维,让包络贴紧目标越来越难,接受率常常随维度指数衰减。

    The accepted \(x\) have density proportional to \(q(x)\cdot\tilde p(x)/(Mq(x))\propto\tilde p(x)\). The acceptance rate is \(Z/M\), with \(Z=\int\tilde p\,dx\) the target's "area": the fraction of the envelope's area that lies under the target, exactly the "box" picture of FIG 2. In \(d\) dimensions it becomes harder and harder to hug the target and the acceptance rate typically decays exponentially with \(d\).

    重要性采样Importance sampling

    $$ \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\) 更重,否则偶尔出现的巨大权重会让方差无穷大。

    After dividing by \(\sum w_i\) ("self-normalization") \(p\) need only be known up to a constant. \(\mathrm{ESS}\in[1,N]\): it equals \(N\) when all weights are equal and 1 when a single sample carries all the weight. \(q\) must be non-zero wherever \(p\) has mass, and ideally have heavier tails than \(p\); otherwise the occasional huge weight makes the variance infinite.

    Metropolis–HastingsMetropolis–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)\)。

    A chain built this way satisfies detailed balance \(\tilde p(x)\,T(x\to x')=\tilde p(x')\,T(x'\to x)\), so \(p\) is its stationary distribution; after convergence the samples approximately follow \(p\). The unconverged beginning is called "burn-in" and is usually discarded. Successive samples are autocorrelated with coefficients \(\rho_k\), and the effective sample size is about \(N/(1+2\sum_{k\ge1}\rho_k)\).

    方法Method需要什么Needs优点Strength弱点Weakness
    逆 CDFInverse CDF\(F^{-1}\)精确、独立、不浪费exact, independent, no waste只适用于一维或结构特别的分布1-D or specially structured distributions only
    拒绝采样Rejection\(\tilde p,\,q,\,M\)精确、独立、\(\tilde p\) 不必归一化exact, independent, \(\tilde p\) unnormalized高维时接受率崩溃acceptance collapses in high dimensions
    重要性采样Importance\(\tilde p,\,q\)不丢样本;可估稀有事件;粒子滤波的基础keeps all samples; handles rare events; basis of particle filters\(q\) 差则权重退化 (ESS 低)weights degenerate (low ESS) when \(q\) is poor
    MCMC\(\tilde p,\,g\)高维也能用;只要比值works in high dimensions; only a ratio needed样本相关、要预热、难判断是否收敛correlated samples, burn-in, hard to judge convergence

    7机器人中的它:从运动模型采样出"香蕉"In robotics: sampling the motion model gives a "banana"

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

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

    A robot drives along a circular arc; every small step has a random heading error and a slight proportional error in length. Apply the one-step noise model 20 times and you get the pose distribution after the whole trip, with no closed form to write down. The simplest approach is to sample: let 400 "possible robots" each walk with their own independent noise.

    This is exactly sample_motion_model in the predict step of a particle filter (Thrun et al., Probabilistic Robotics, Ch. 5). We also draw the approach of Lesson 4.5: linearize with Jacobians along the nominal path and propagate the covariance \(F\Sigma F^\top+Q\) step by step, giving a 2σ ellipse. With small heading noise the ellipse is accurate; with large heading noise and a long drive the samples bend into a banana that no ellipse can enclose.

    FIG 5采样运动模型 vs 线性化高斯Sampling the motion model vs the linearized Gaussian
    黄点是 400 个采样位姿,青色虚线是名义(无噪声)路径,紫色椭圆是沿名义路径线性化传播得到的 2σ 椭圆(4.5 节)。若分布真是高斯,约 86.5% 的点应落在 2σ 椭圆内;读数里的"椭圆内比例"偏离这个值越多,高斯假设越差。 The amber dots are 400 sampled poses, the cyan dashed line is the nominal (noise-free) path, and the purple ellipse is the 2σ ellipse from linearized propagation along the nominal path (Lesson 4.5). If the distribution were truly Gaussian, about 86.5% of the points would fall inside the 2σ ellipse; the further the "inside ellipse" readout is from that, the worse the Gaussian assumption.
    试一试Try this
    • 把航向噪声调到 3°/√m 以下:样本云几乎是个椭圆,高斯近似很好(椭圆内比例约 80%–90%)。Set the heading noise below 3°/√m: the cloud is nearly an ellipse and the Gaussian approximation is fine (inside-ellipse fraction roughly 80%–90%).
    • 把航向噪声调大(≥ 12°/√m)、距离拉到 7 m 以上:样本弯成"香蕉",椭圆内比例掉到 78% 以下,而且样本均值偏离了名义终点(椭圆始终以名义终点为中心)。Raise the heading noise (≥ 12°/√m) and the distance to 7 m or more: the samples bend into a "banana", the inside-ellipse fraction falls below 78%, and the sample mean drifts away from the nominal endpoint (where the ellipse stays centered).
    真实系统里In real systems

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

    "Sample from a model for which you can write a forward simulation but not a density" is the most common use of sampling in robotics: the predict step of a particle filter (4.10); domain randomization, sampling masses, friction and latency in simulation to train policies that are more robust in reality; sampling-based model-predictive control such as MPPI and CEM, which in every control cycle sample hundreds of noisy trajectories: MPPI weights them by \(\exp(-\text{cost}/\lambda)\), which are importance weights, while CEM keeps the lowest-cost few and refits its sampling distribution to them; RRT, which samples the configuration space; and Monte-Carlo tree search, which estimates values by random rollouts.

    8常见误区Common pitfalls

    误区 1:包络没盖住目标Pitfall 1: an envelope that does not cover the target

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

    FIG 2 showed it: where \(Mq(x)<\tilde p(x)\) the sample distribution is silently truncated while the program still outputs normal-looking numbers. Rejection sampling needs the upper bound of \(\tilde p/q\) verified over the whole domain, not just at a few points you happened to try.

    误区 2:重要性权重退化Pitfall 2: degenerate importance weights

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

    When \(q\) and \(p\) differ greatly ESS collapses to about 1: of 1000 samples only one matters, and the estimate looks "converged" while it is really that single sample. Always report ESS alongside the estimate. In particle filters this is weight degeneracy, the problem that resampling exists to fix (Lesson 4.10).

    误区 3:把 MCMC 的相邻样本当作独立样本Pitfall 3: treating consecutive MCMC samples as independent

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

    Neighbouring samples are strongly correlated; 10,000 steps may have an ESS of only a few dozen, and error bars computed as if they were independent badly understate the uncertainty. Beware also of "looks converged": a snail-step chain sits comfortably on one hill, every single-chain diagnostic looks fine, and the other hill is never visited. The remedy: run several chains from different starting points and compare them.

    误区 4:随机数用得不规范Pitfall 4: careless random numbers

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

    In robot systems fix and log the random seed, or you cannot reproduce a failed simulation; do not let parallel threads share one global generator, and do not give every thread the same seed, or all "independent" particles walk the same path. The pages of this course use seeded generators, so a figure shows the same data every time you open it.

    9小练习Check yourself

    1. 指数分布 \(F(x)=1-e^{-\lambda x}\),\(\lambda=2\)。用逆 CDF 法,均匀数 \(u=0.5\) 对应的样本是?1. Exponential distribution \(F(x)=1-e^{-\lambda x}\) with \(\lambda=2\). With inverse-CDF sampling, which sample does the uniform number \(u=0.5\) give?
    \(x=-\ln(1-u)/\lambda=-\ln(0.5)/2=0.693/2\approx0.347\)。它正好是这个分布的中位数,因为 \(u=0.5\) 就是 CDF 的 50% 处。\(x=-\ln(1-u)/\lambda=-\ln(0.5)/2=0.693/2\approx0.347\). It is the median of the distribution, since \(u=0.5\) is exactly the 50% level of the CDF.
    2. 目标曲线 \(\tilde p\) 下的面积是 2;提议是 \([0,10]\) 上的均匀分布,包络是高度 0.5 的水平线(\(\tilde p\le 0.5\) 处处成立)。接受率是多少?2. The area under the target \(\tilde p\) is 2; the proposal is uniform on \([0,10]\) and the envelope is a horizontal line at height 0.5 (\(\tilde p\le0.5\) everywhere). What is the acceptance rate?
    包络下的面积是 \(0.5\times10=5\),目标下的面积是 2,飞镖落在目标下的比例 \(=2/5=40\%\)。The area under the envelope is \(0.5\times10=5\), the area under the target is 2, so the fraction of darts landing under the target is \(2/5=40\%\).
    3. 四个重要性权重(已归一化)为 \((0.7,\,0.1,\,0.1,\,0.1)\)。有效样本数 ESS 约为?3. Four (normalized) importance weights are \((0.7,\,0.1,\,0.1,\,0.1)\). The effective sample size ESS is about?
    权重和为 1,\(\sum w^2=0.49+0.01\times3=0.52\),\(\mathrm{ESS}=1/0.52\approx1.92\)。4 个样本只相当于不到 2 个理想样本,因为 0.7 的那个独占了大部分权重。The weights sum to 1, \(\sum w^2=0.49+0.01\times3=0.52\), so \(\mathrm{ESS}=1/0.52\approx1.92\). Four samples are worth fewer than two ideal ones, because the 0.7 sample hogs most of the weight.
    4. 一条 Metropolis 链的接受率是 99%。最可能的情况是?4. A Metropolis chain has a 99% acceptance rate. What is most likely going on?
    步子太小时,\(x'\) 与 \(x\) 的密度几乎相同,比值接近 1,当然几乎全被接受,可是走得慢。(目标恰好是均匀分布也会得到 100%,但那是个特例。)要的是"又远又被接受"的折中,一般让接受率落在 20%–50%。With tiny steps the densities at \(x'\) and \(x\) are nearly equal, the ratio is close to 1, and almost everything is accepted, but the walk is slow. (A uniform target would also give 100%, but that is a special case.) You want the compromise "far and still accepted", typically an acceptance rate of 20%–50%.
    5. 用 \(N=100\) 个独立样本做蒙特卡洛估计,标准误差约 0.1。改用 \(N=10{,}000\),标准误差约为?5. A Monte-Carlo estimate with \(N=100\) independent samples has standard error about 0.1. With \(N=10{,}000\), the standard error is about?
    误差 \(\propto1/\sqrt N\)。样本多 100 倍,误差只缩小 \(\sqrt{100}=10\) 倍,得 0.01。想再缩小 10 倍,又要多 100 倍的样本,所以蒙特卡洛适合"够用就行"的精度,不适合高精度。The error scales as \(1/\sqrt N\). With 100 times more samples the error shrinks only \(\sqrt{100}=10\) times, to 0.01. Another factor of 10 costs another factor of 100 in samples, which is why Monte Carlo is good for "good enough" accuracy, not for high precision.

    10延伸阅读Further reading

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