采样方法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.
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.
- 拖动 \(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.
逆 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.
- 省飞镖:让 \(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".
- 好的提议:估计均值时,让 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.
- 蜗牛:步长 ≤ 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).
接受率接近 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
误差按 \(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\):
例:指数分布 \(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\) 的密度正比于 \(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
用分母 \(\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
这样构造的链满足细致平衡 \(\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.
- 把航向噪声调到 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).
"从一个你写得出前向仿真、却写不出密度的模型里采样"是机器人里最常用的采样用法:粒子滤波的预测步 (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
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.
\(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).
相邻样本高度相关,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.
机器人系统里要固定并记录随机种子,否则你无法复现一次失败的仿真;并行线程不要共享同一个全局生成器,也不要给每个线程用同一个种子,否则所有"独立"的粒子其实在走同一条路。本课网页里用了带种子的生成器,所以每次打开同一个图,数据都相同。
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
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