M1-T3 — Simulating Randomness: Sampling and Monte Carlo

Learning objectives

After this lesson, you should be able to:

• Explain what a pseudo-random number generator and a seed are, and why a simulation must be reproducible.
• Estimate a probability by Monte Carlo simulation and state the estimate with its standard error.
• Derive the standard error sqrt(p(1 − p)/n) of an estimated probability.
• Compute how many runs are needed for a required accuracy, and recognise what more runs cannot fix.

The intuition

The previous lesson read a probability off a table of 200 real runs. Real runs are slow. If we have a model of the robot (the command, plus noise of a known size), a computer can produce a run in microseconds, and the recipe for any probability becomes:

1. Simulate the experiment once, drawing every noisy quantity at random.
2. Record whether the event of interest happened.
3. Repeat $n$ times and report the fraction of runs in which it happened.

This is a Monte Carlo method: computing by random sampling. It works for events of any shape and for models of any complexity, which is why it is everywhere in robotics. Its one weakness is visible in the figure: the answer is itself random.

With 10 runs the three estimates disagree wildly. With 10 000 they agree to within about half a percentage point. Two facts explain the picture, and the derivation below proves both:

• The estimate converges to the truth. This is the law of large numbers, and it is what gives the frequency reading of probability its meaning.
• It converges slowly. The typical error shrinks like $1/\sqrt{n}$: four times the runs for half the error.

Random numbers from a deterministic machine

A computer has no dice. A pseudo-random number generator is an ordinary function that turns an internal number, its state, into a new state and an output between 0 and 1, in a way so scrambled that the outputs pass statistical tests for randomness. The starting state is the seed. The same seed always yields the same sequence; a different seed yields an unrelated one.

That determinism is a feature. A simulated failure that you cannot reproduce cannot be debugged. With a seed, “run 4 812 hit the wall” is a test case you can replay under a debugger, hand to a colleague, and keep as a regression test after the fix.

Fix the seed to make an experiment repeatable. Change the seed to see how much of your result was luck. A conclusion that holds for one seed only is not a conclusion.

Definitions

Sample

One value drawn from a distribution: one simulated heading error, one simulated end point.

Seed

The number that initializes a generator. Rng(42) always produces the same stream.

Indicator $I_i$

A variable that is 1 if the event happened in run $i$ and 0 if not. Counting successes is adding indicators.

Estimator $\hat{p}$

The rule that turns data into an estimate. Here $\hat{p} = k/n$, the fraction of the $n$ runs in which the event happened. Because the runs are random, $\hat{p}$ is random too.

Standard error

The standard deviation of an estimator: the typical distance between the estimate and the true value over repetitions of the whole experiment.

Bias

A systematic difference between what the estimator converges to and the quantity you wanted. More runs do not reduce it.

Derivation: why the error shrinks like 1/√n

The derivation uses two summaries of a random quantity that M4-T1: Expectation and Variance treats in full. For now, two facts are enough. The expectation $\mathbb{E}[X]$ is the long-run average of $X$. The variance $\mathrm{Var}(X)$ is the long-run average of the squared distance of $X$ from its expectation; its square root is the standard deviation, the typical size of a deviation.

1. $$\mathbb{E}[I_i] = 1 \cdot p + 0 \cdot (1-p) = p$$

   An indicator is 1 with probability $p$ and 0 otherwise, so its long-run average is $p$.

2. $$\mathrm{Var}(I_i) = (1-p)^2\,p + (0-p)^2\,(1-p) = p\,(1-p)$$

   Squared distance from the mean, weighted by probability, then simplified. It is largest at $p = 0.5$ (the outcome is least predictable) and zero at $p = 0$ or 1.

3. $$\mathbb{E}[\hat{p}] = \frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[I_i] = p$$

   The expectation of a sum is the sum of the expectations (linearity). So the estimator is unbiased: on average it is exactly right.

4. $$\mathrm{Var}\!\left(\sum_{i=1}^{n} I_i\right) = \sum_{i=1}^{n}\mathrm{Var}(I_i) = n\,p\,(1-p)$$

   For independent runs, variances add. This is where the requirement enters that each run uses fresh random numbers and does not influence another.

5. $$\mathrm{Var}(\hat{p}) = \frac{1}{n^{2}}\, n\,p\,(1-p) = \frac{p\,(1-p)}{n}$$

   Multiplying a quantity by a constant $c$ multiplies its variance by $c^{2}$; here $c = 1/n$.

6. $$\mathrm{se}(\hat{p}) = \sqrt{\mathrm{Var}(\hat{p})} = \sqrt{\frac{p\,(1-p)}{n}}$$

   The standard error is the standard deviation of the estimator.

Steps 3 and 6 together are the law of large numbers for this case: the estimator is centred on $p$ and its spread goes to zero as $n$ grows, so it converges to $p$.

How many runs?

To make the 95% interval no wider than $\pm h$, set $1.96\,\mathrm{se} = h$ and solve for $n$:

Monte Carlo error falls as 1/√n: to gain one more decimal digit of accuracy you need one hundred times as many runs.

A worked example

Rover turns towards a dock 2.00 m away and drives. Its heading error is Gaussian with standard deviation 1.5°. To keep the example exact, ignore the distance noise: the sideways error at the dock is then $y = 2\sin\theta$, and Rover docks when $|y| \le 0.10$ m. What is the probability of docking?

1. Exact answer, for comparison. $|2\sin\theta| \le 0.10$ means $|\theta| \le \arcsin(0.05) = 2.866^{\circ}$, which is $2.866/1.5 = 1.911$ standard deviations. The area of a Gaussian within ±1.911 standard deviations is 0.9440. (Most real problems have no such formula; this one was chosen so that the simulation can be checked.)
2. Simulate. Draw $\theta$ from the generator, compute $y$, count the run if $|y| \le 0.10$.
3. Report with a standard error. With 100 runs: 95 docked, $\hat{p} = 0.95$, $\mathrm{se} = \sqrt{0.95 \times 0.05 / 100} = 0.022$. The honest statement is “0.95 ± 0.04”: compatible with a true value anywhere from 0.91 to 0.99.
4. Increase n. With 10 000 runs: $\hat{p} = 0.9447$, $\mathrm{se} = 0.0023$, and the actual error is 0.0007.
5. Check the standard error itself. Repeat the whole experiment with 200 different seeds and measure the spread of the 200 estimates. For 100, 400 and 1600 runs it is 0.0234, 0.0117 and 0.0057: each fourfold increase in runs halves it, and each value matches $\sqrt{p(1-p)/n}$.
6. Plan. To know the docking probability within ±0.01: $n = 1.96^{2} \times 0.944 \times 0.056 / 0.01^{2} \approx 2030$ runs (the program, which keeps every digit of the exact probability, prints 2033).

import math
from robotprob import Rng, gaussian_cdf

D, sigma_theta, edge = 2.0, math.radians(1.5), 0.10
exact = 2 * gaussian_cdf(math.asin(edge / D), 0.0, sigma_theta) - 1

rng = Rng(42)                                   # the seed makes this experiment repeatable
for n in (100, 1000, 10000):
    k = 0
    for _ in range(n):
        theta = rng.gauss(0.0, sigma_theta)     # one simulated run
        if abs(D * math.sin(theta)) <= edge:
            k += 1
    p = k / n
    se = math.sqrt(p * (1 - p) / n)
    print(n, 'runs: p =', round(p, 4), ' standard error =', round(se, 4), ' actual error =', round(p - exact, 4))
print('exact probability =', round(exact, 4))
print('runs for +-0.01    =', math.ceil(1.96 ** 2 * exact * (1 - exact) / 0.01 ** 2))

100 runs: p = 0.95  standard error = 0.0218  actual error = 0.006
1000 runs: p = 0.926  standard error = 0.0083  actual error = -0.018
10000 runs: p = 0.9447  standard error = 0.0023  actual error = 0.0007
exact probability = 0.944
runs for +-0.01    = 2033

The second program repeats the experiment under 200 seeds to measure the spread of the estimate directly.

import math
from robotprob import Rng, std

D, sigma_theta, edge, exact = 2.0, math.radians(1.5), 0.10, 0.9440
for n in (100, 400, 1600):
    estimates = []
    for seed in range(1, 201):
        rng = Rng(seed)
        k = sum(1 for _ in range(n) if abs(D * math.sin(rng.gauss(0.0, sigma_theta))) <= edge)
        estimates.append(k / n)
    print(n, 'runs: spread of 200 estimates =', round(std(estimates, 0), 4), ' formula =', round(math.sqrt(exact * (1 - exact) / n), 4))

100 runs: spread of 200 estimates = 0.0234  formula = 0.023
400 runs: spread of 200 estimates = 0.0117  formula = 0.0115
1600 runs: spread of 200 estimates = 0.0057  formula = 0.0057

In robotics

Particle filters are Monte Carlo inside the robot

A particle filter represents the belief about the robot's pose by a set of samples, typically a few hundred to a few thousand, instead of a formula. Everything derived here applies to it: its estimates carry a $1/\sqrt{n}$ error, more particles cost more computation for slowly improving accuracy, and its random choices should be seeded so that a failure can be replayed. Module 9 builds it; M9-T1: Monte Carlo Belief and Importance Sampling starts from this lesson.

Evaluating a robot in simulation

“The planner succeeded in 47 of 50 simulated trials” is a Monte Carlo estimate of a success probability: $\hat{p} = 0.94$ with a standard error of 0.034. It cannot distinguish a 90% planner from a 97% one. The Robot Learning course treats this in Evaluating a Policy: held-out scenarios, intervals and paired comparisons.

Rare failures need many runs

For a small probability $p$, the relative error of the estimate is $\mathrm{se}/p \approx 1/\sqrt{n\,p}$. To know a failure probability of 1 in 10 000 to within 10% takes about $100/p = 10^{6}$ runs, and zero failures in 1000 runs only shows that the probability is probably below about 0.3%. This is why safety cases cannot rest on plain simulation counts alone.

Randomized algorithms

Sampling-based planners (PRM and RRT) and robust estimators (RANSAC) draw random samples by design. Their results vary from run to run, so they are compared over many seeds, and their random generators are seeded in tests.

Ground truth for approximations

When a filter approximates how uncertainty passes through a nonlinear function, the reference answer is a large Monte Carlo run. The tests of this course's library do exactly that: the unscented transform of module 8 is accepted because it matches 400 000 samples where plain linearization does not.

Common misconceptions

Key takeaways

1. Monte Carlo estimates a probability as the fraction of simulated runs in which the event happens. It needs only a model you can sample from.
2. A seeded generator makes a random experiment exactly repeatable. Keep the seed of every failure.
3. The estimate is unbiased and its standard error is sqrt(p(1 − p)/n). Always report it, or a 95% interval of about ±2 standard errors.
4. Error falls as 1/√n: half the error costs four times the runs. The runs needed for a half-width h are 1.96² p(1 − p)/h².
5. Sampling error and model error are different. Runs reduce the first; only a better model or real data reduce the second.

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 (the particle filter as a sample-based representation). Paid book.
• J. Tsitsiklis and P. Jaillet, MIT RES.6-012 Introduction to Probability , lecture 18 (inequalities, convergence and the weak law of large numbers). Free videos and notes, licence CC BY-NC-SA 4.0.
• J. Blitzstein and J. Hwang, Introduction to Probability, 2nd edition, CRC Press, 2019 , chapter 10 (inequalities and limit theorems). Free to read online; the book is copyrighted.
• N. Metropolis and S. Ulam, “The Monte Carlo Method”, Journal of the American Statistical Association 44 (247), 1949 . The paper that named the method. Paid journal article.
• E. B. Wilson, “Probable Inference, the Law of Succession, and Statistical Inference”, Journal of the American Statistical Association 22 (158), 1927 . The Wilson interval. Paid journal article.

Knowledge check

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

What’s next

Module 1 has given you the vocabulary, the rules and an experimental method. Module 2 turns to the objects those rules act on: random variables and their distributions.

Related content

Nothing below is required to finish this lesson.

• Before
  M1-T2: Events, Sample Spaces and the Rules of Probability

  Describe an uncertain situation as a sample space and events, and use the three axioms, the complement rule and the addition rule.

• Lab
  M1-L1: Noisy Robot Motion Simulator

  The same command gives a different end position every run. Simulate noisy motion, measure the spread, and predict how often the robot stops inside the docking zone.

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

• Physical AI: State, observation, belief
• SLAM M2-T4: Probability and Uncertainty