M9-T1 — Monte Carlo Belief and Importance Sampling

Learning objectives

After this lesson, you should be able to:

• Approximate an expectation by a Monte Carlo average and state how its error scales with the number of samples.
• Derive importance sampling and self-normalised weights for a target known up to a constant.
• Compute the effective sample size and explain why a poor proposal wastes samples.
• Choose a proposal and recognise the choice made by the bootstrap particle filter.

The intuition

The Bayes filter of module 6 multiplies a prior by a likelihood and normalises. On a grid that is easy in one dimension and impossible in six. A Gaussian filter does it in closed form but only for Gaussians. A particle filter takes a third route: represent the belief by samples and do the multiplication with weights.

Suppose a robot's prior belief is N(0, 2²) m and one reading, z = 3 m with 0.5 m of noise, arrives. The posterior is the prior times the likelihood, narrow and centred near 2.8 m. We cannot sample from it directly (in general we can only evaluate prior times likelihood at a point). But we can sample from the prior, and for each sample ask how well it explains the reading. A sample at 0.1 m explains it very badly; one at 2.9 m explains it very well. Give each sample a weight equal to its likelihood, and the weighted cloud is the posterior.

The left panel shows the result with 100 samples. Most of the prior samples are wasted: they lie where the posterior has no mass, and their weights are almost zero. A dozen or so samples near 2 to 4 m carry the answer. The effective sample size counts these: 12.4 of 100 for the prior proposal. The right panel shows how the choice of the proposal changes it: a proposal that is a good guess of the posterior (a Gaussian at 2.8 m, 0.6 m wide) has 93.7, a proposal that is too narrow and in the wrong place (0.5 m around the origin) has 2.2: the whole answer rests on two samples.

Despite the waste, the estimate is good: the posterior mean from 100 prior samples is within 0.11 m of the true 2.82 m, and with 10 000 the error is 1 cm. The error falls as 1/√N, for any dimension: that is why Monte Carlo is the method of choice where grids fail.

Importance sampling turns the problem of sampling from a distribution you can only evaluate into the problem of sampling from one you can draw from, at the price of weights. The closer the proposal is to the target, the less is wasted.

Definitions

Monte Carlo estimate

$\mathbb{E}[f(X)] \approx \frac1N\sum_i f(x_i)$ with $x_i$ drawn from the distribution.

Target $p$, proposal $q$

The distribution we want to represent, and the one we sample from. $q$ must be positive wherever $p$ is.

Importance weight $w_i$

$p(x_i)/q(x_i)$: how much more (or less) the target likes the sample than the proposal did.

Self-normalised weights $\tilde w_i$

$w_i/\sum_j w_j$: work when $p$ is known only up to a constant, as a posterior is.

Effective sample size

$N_\mathrm{eff} = 1/\sum_i\tilde w_i^2$: $N$ for equal weights and 1 when one sample has all the weight.

Bootstrap choice

The proposal is the prior (the motion model) and the weight the likelihood: $w_i = p(z\mid x_i)$. Simple, and what the standard particle filter does.

Derivation: why the weighted samples represent p

Monte Carlo error

1. $$\hat\mu = \frac1N\sum_i f(x_i),\quad x_i\ \text{independent}$$

   The average of N independent copies of $f(X)$.

2. $$\mathbb{E}[\hat\mu]=\mu,\qquad \mathrm{Var}(\hat\mu)=\frac{\mathrm{Var}(f)}{N}$$

   The estimate is unbiased, and the variance of an average of independent terms is the variance of one divided by N (the variance lesson of module 4).

3. $$\mathrm{std}(\hat\mu)=\frac{\sigma_f}{\sqrt N}$$

   Hence the 1/√N: to halve the error, take four times as many samples. The dimension of x does not appear.

Importance sampling

1. $$\mathbb{E}_p[f]=\int f(x)\,p(x)\,dx=\int f(x)\,\frac{p(x)}{q(x)}\,q(x)\,dx=\mathbb{E}_q\!\left[f(X)\,w(X)\right]$$

   Multiply and divide by $q$ (needs $q>0$ wherever $p>0$). The expectation under p becomes one under q of the weighted function.

2. $$\mathbb{E}_p[f]\approx\frac1N\sum_i w_i f(x_i)$$

   The Monte Carlo average of the weighted function over samples from q. Unbiased when p is normalised.

3. $$\mathbb{E}_p[f]=\frac{\mathbb{E}_q[w\,f]}{\mathbb{E}_q[w]}\approx\frac{\sum_i w_if(x_i)}{\sum_i w_i}$$

   If p is known only up to a constant, $p = \tilde p/Z$, use $\tilde p/q$ as the weight: $Z$ appears in numerator and denominator and cancels. Slightly biased for finite N, consistent as N grows.

The effective sample size

The variance of the weighted estimate is governed by the variance of the weights. For equal weights $\tilde w_i = 1/N$ the sum $\sum\tilde w_i^2 = 1/N$ and $N_\mathrm{eff}=N$. If one weight is 1 and the others 0, $N_\mathrm{eff}=1$. In between, $N_\mathrm{eff}$ is the number of equally weighted samples that would give the same variance (to a good approximation).

The bootstrap choice

For the posterior $p(x\mid z)\propto p(z\mid x)\,p(x)$ with the prior as proposal, $w\propto p(z\mid x)p(x)/p(x) = p(z\mid x)$: the weight is simply the likelihood. No density of the prior is ever needed, only the ability to sample it, which is why the motion model of a robot (easy to sample, hard to write as a density) fits.

The importance weight of a prior sample is its likelihood. The prior samples are where the robot might be; the weights say how well each place explains the reading.

A worked example

The posterior of the figure: prior N(0, 2²), one reading z = 3 with 0.5 m of noise. Exactly, the posterior has mean 2.8235 and standard deviation 0.485 (fusion of Gaussians). We estimate its mean by importance sampling with four proposals, 500 repetitions of N = 100 samples each, and then study the error against N.

1. Proposals. The prior N(0, 2) (the bootstrap choice) gives an rms error of 0.114 m and an effective sample size of 12.4. The wider N(0, 5) is no better: 0.109 and 11.8 (it also wastes samples far out). The narrow N(0, 0.5) is badly wrong: 1.63 m and 2.2, since almost no sample reaches the posterior. The good proposal N(2.8, 0.6) gives 0.043 m and 93.7.
2. Number of samples. With the prior as proposal, N = 10, 100, 1 000 and 10 000 give rms errors of 0.67, 0.109, 0.031 and 0.011 m. Multiplying by √N gives 2.1, 1.09, 0.97 and 1.07: the error falls as 1/√N, as the derivation says (the first value is high because with ten samples and an effective size of 1.8 the self-normalised estimate is poor).
3. Effective size. The effective sample size is about 12% of N at every N (12.3 of 100, 121 of 1 000, 1 211 of 10 000): the waste is a fixed fraction, set by how much narrower the posterior is than the proposal.

import math
from robotprob import Rng, gaussian_pdf, normalize_weights, effective_sample_size, mean

# target: the posterior of a robot position with prior N(0, 2^2) and one reading z = 3 from a sensor with noise 0.5
prior_sd, z, noise = 2.0, 3.0, 0.5
unnorm = lambda x: gaussian_pdf(x, 0.0, prior_sd) * gaussian_pdf(z, x, noise)      # prior times likelihood: known up to a constant
prec = 1 / prior_sd**2 + 1 / noise**2
exact_mean, exact_sd = (z / noise**2) / prec, math.sqrt(1 / prec)
print('exact posterior: mean', round(exact_mean, 4), ' std', round(exact_sd, 4))

def estimate(rng, n, q_mean, q_sd):
    xs = [rng.gauss(q_mean, q_sd) for _ in range(n)]
    w = [unnorm(x) / gaussian_pdf(x, q_mean, q_sd) for x in xs]       # importance weights p / q
    w = normalize_weights(w)
    return sum(wi * xi for wi, xi in zip(w, xs)), effective_sample_size(w)

print('proposal q            N     rms error of the mean   mean ESS')
for name, qm, qs in [('prior N(0, 2)', 0.0, 2.0), ('wide N(0, 5)', 0.0, 5.0), ('narrow N(0, 0.5)', 0.0, 0.5), ('good N(2.8, 0.6)', 2.8, 0.6)]:
    rng = Rng(3)
    for n in (100,):
        res = [estimate(rng, n, qm, qs) for _ in range(500)]
        rms = math.sqrt(mean([(m - exact_mean) ** 2 for m, e in res]))
        print(name.ljust(20), str(n).rjust(5), str(round(rms, 4)).rjust(16), str(round(mean([e for m, e in res]), 1)).rjust(14))

print('error against the number of samples (proposal = the prior):')
rng = Rng(4)
for n in (10, 100, 1000, 10000):
    reps = 200 if n < 10000 else 50
    res = [estimate(rng, n, 0.0, prior_sd) for _ in range(reps)]
    rms = math.sqrt(mean([(m - exact_mean) ** 2 for m, e in res]))
    print(str(n).rjust(7), ' rms error', str(round(rms, 4)).rjust(7), ' rms * sqrt(N)', str(round(rms * math.sqrt(n), 2)).rjust(6), ' mean ESS', round(mean([e for m, e in res]), 1))

exact posterior: mean 2.8235  std 0.4851
proposal q            N     rms error of the mean   mean ESS
prior N(0, 2)          100           0.1142           12.4
wide N(0, 5)           100           0.1092           11.8
narrow N(0, 0.5)       100           1.6279            2.2
good N(2.8, 0.6)       100           0.0429           93.7
error against the number of samples (proposal = the prior):
     10  rms error   0.674  rms * sqrt(N)   2.13  mean ESS 1.8
    100  rms error  0.1094  rms * sqrt(N)   1.09  mean ESS 12.3
   1000  rms error  0.0306  rms * sqrt(N)   0.97  mean ESS 121.4
  10000  rms error  0.0107  rms * sqrt(N)   1.07  mean ESS 1210.7

In robotics

Particle filters

The localisation filter of the next lessons is this method applied at every step: the belief at time t−1 is a set of weighted samples; each is moved by the motion model (the proposal, so it is the prior), weighted by the likelihood of the new reading, and the set is resampled so that the weights return to equal. The Navigation course shows the same filter at work in its localisation lessons.

Monte Carlo in robotics beyond filtering

• Estimating the probability of a collision along a path with uncertain obstacles.
• Expected cost in model-predictive control by sampling noise trajectories.
• Evaluating the risk of a plan under a belief (the sample average of the cost).
• Bayesian inference of parameters when the posterior has no closed form.

Choosing and checking the proposal

The prior is the easiest proposal and wastes samples when the likelihood is sharper than the prior (an accurate sensor): the effective size falls. Better proposals mix in the measurement (the “optimal” proposal and its approximations, Module 9 lesson 2) at the price of a more complicated derivation. Always monitor the effective sample size: a value near 1 means the filter has effectively collapsed.

Heavy tails

The proposal should have heavier tails than the target. A proposal that is too narrow (the 0.5 m one) misses the regions where the target has mass, with no warning: the weights of the samples you do have are fine, but the estimate is badly biased. The failure is silent, which is one reason to use broad priors, to inject random samples, and to compare against the sensor.

Common misconceptions

Key takeaways

1. A Monte Carlo average of N independent samples is unbiased with an error that falls as 1/√N, whatever the dimension.
2. Importance sampling draws from a proposal q and weights by p/q. Self-normalised weights work when p is known only up to a constant.
3. The effective sample size 1/Σw̃² counts the samples doing the work. It falls when the proposal is far from the target.
4. In the bootstrap choice the proposal is the prior (the motion model) and the weight is the likelihood.
5. Heavy-tailed, broad proposals are safe; narrow ones in the wrong place fail silently.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• S. Thrun, W. Burgard and D. Fox, Probabilistic Robotics, MIT Press, 2005, section 4.3 (particle filters) and 4.3.1 (importance sampling). Paid book.
• S. Särkkä and L. Svensson, Bayesian Filtering and Smoothing, 2nd edition, 2023 , chapter 11 (particle filtering), including importance sampling. Free to read as the authors’ online PDF; the book is copyrighted.
• T. D. Barfoot, State Estimation for Robotics, 2nd edition, 2024 , the section on the particle filter in chapter 4. Free to read as the author’s PDF; the book is copyrighted.
• C. Stachniss, “Particle Filter and Monte Carlo Localization” (video lecture) . Free to watch.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Now the algorithm: apply importance sampling at every time step, resample to keep the weights usable, and see the particle filter localise a robot.

Related content

Nothing below is required to finish this lesson.

Other courses on this site treat neighbouring topics from their own angle:

• SLAM M3-T8: Particle Filters and FastSLAM
• SLAM M3-L4: Explore Particle-Filter Localization