M12-T2 — Monte Carlo Dispersion Analysis

Learning objectives

After this lesson, you should be able to:

• Compute the number of failure-free runs that demonstrates a probability at a stated confidence.
• Compute an exact lower confidence bound on a success probability from a count of failures.
• Say why the observed fraction of successes is not the demonstrated probability, and how many runs separate them.
• Explain the effect of stopping a campaign as soon as it passes, and read the tail of a result rather than a Gaussian fit.

The intuition

Probability M1-T3: Simulating Randomness: Sampling and Monte Carlo shows how to draw random numbers and estimate an average. A GNC campaign uses the same idea for a different purpose. The thing estimated is not an average miss but the probability of meeting a requirement, which is a count: of 1000 flights, 990 landed on the pad. The requirement for Hopper reads "within 3 m in at least 99% of flights", so the campaign must produce a statement of the form "with 95% confidence, the probability exceeds 99%". That sentence has three numbers, and each of them costs runs.

What to disperse

Every quantity that differs between two real flights is drawn at random in each simulated flight: the wind speed and direction, the initial position error, the sensor biases and noise, the mass and its distribution, the engine thrust scale, the actuator lag. The distribution of each comes from a measurement, a datasheet or a stated assumption, and the assumption is written down, because a campaign gives an answer about those distributions and no others. Dispersed together, they produce the combinations that nobody thought to test by hand, such as a high headwind, a low engine thrust and a late navigation fix on the same flight.

Faults are not dispersed in the same draw. A beacon failure with probability 0.4% per flight is a separate event with its own probability, and the campaign needs a stated rule for how it counts toward the requirement. In the example it counts as an ordinary flight with a rare cause.

Why the count matters more than the fraction

Suppose 1000 flights all land. The observed fraction is 100%, but the data are equally consistent with a true probability of 99.7%, because at that probability the chance of seeing no failure in 1000 flights is still 5%. So the honest statement is the bound, 99.70%, and it is the bound that is compared with the requirement. If five flights fail, the observed fraction is 99.5% (above the requirement), but the bound is 98.95% (below it): the campaign has not demonstrated the requirement, and it is not evidence that the requirement fails either. It says that more runs are needed.

The figure carries the practical consequence. A vehicle that never fails demonstrates 99% after 299 runs. A vehicle whose failure fraction is 0.3% needs about 650 runs to be shown to comply, at 0.6% about 1700, and at 0.8% about 6500. Each step toward the line multiplies the size of the campaign, and at 0.99% (on the line) no number of runs is enough. The smaller the margin, the larger the campaign.

A campaign demonstrates a bound, not a fraction. The number of runs is chosen before it starts, from the probability, the confidence and the margin that the design is expected to have.

Definitions

Requirement probability $p_0$

The success probability to be demonstrated: 0.99 for Hopper.

Confidence $1-\alpha$

The probability that the method of calculation gives a bound that is correct. 0.95 here, so $\alpha = 0.05$.

Failures $k$ of $n$ runs

A flight that lands farther than 3 m from the pad, or fails to land, is a failure.

Lower confidence bound $p_L$

A number such that, whatever the true probability, the bound exceeds it with probability at most $\alpha$. The requirement is demonstrated when $p_L \ge p_0$.

Dispersion

Drawing each uncertain input from its distribution for each run.

The count of failures in $n$ independent runs, each failing with probability $q = 1 - p$, follows the binomial law. The exact lower bound $p_L = 1 - q_U$ uses the largest failure probability that still makes the observed count likely:

With no failures the equation is a power and has a closed form. The number of failure-free runs that demonstrates $p_0$ at confidence $1-\alpha$ is:

Derivation: the failure-free count and the stopping problem

The failure-free count first.

1. $$P(\text{no failure in } n \text{ runs}) = p^{\,n}$$

   The runs are independent and each succeeds with probability $p$.

2. $$p^{\,n} \le \alpha \quad\text{for every } p \le p_0 \iff p_0^{\,n} \le \alpha$$

   To demonstrate $p_0$, a vehicle with any probability at or below $p_0$ must have a chance of at most $\alpha$ of passing. The worst case is $p = p_0$, since $p^n$ grows with $p$.

3. $$n \ge \frac{\ln\alpha}{\ln p_0}$$

   Take logarithms; $\ln p_0$ is negative, which turns the inequality around. For $p_0 = 0.99$ and $\alpha = 0.05$: 298.07, so 299 runs.

Now the effect of repeated looks. A test with a single look at $n$ runs accepts a non-compliant vehicle with probability at most $\alpha$. A test that looks at $n_1 < n_2 < \dots < n_L$ and accepts as soon as any look passes accepts if any of the looks accepts.

1. $$A_j = \{\text{look } j \text{ passes}\},\qquad P(A_j) \le \alpha$$

   Each look on its own has a false-acceptance probability of at most $\alpha$ for a vehicle at the requirement.

2. $$P(\text{accepted}) = P\Bigl(\bigcup_{j=1}^{L} A_j\Bigr) \;\ge\; \max_j P(A_j)$$

   The union of events is at least as probable as the most probable of them, and at most the sum of probabilities. The looks are strongly dependent (they share runs), so the union lies between $\alpha$ and $L\alpha$.

3. $$\text{so the true confidence is below } 1-\alpha$$

   A campaign stopped at the first pass has a false-acceptance rate above the stated 5%. In the example the rate rises from 1.2% (one look) to 6.3% (15 looks).

The remedy is to fix the number of runs before the campaign, or to use a sequential test designed for repeated looks with a stricter bound at each.

For zero failures, $n \ge \ln\alpha/\ln p_0$ runs. For any count, the bound comes from the binomial law, not from the observed fraction. And the number of looks is part of the test.

A worked example

A dispersed Hopper campaign: the linear-loop landing miss has a spread of 0.95 m per axis (inside 3 m with probability 99.32% when Gaussian), and a beacon fault in 0.4% of flights adds a 2.5 m offset in a random direction.

1. Failure-free runs. To demonstrate 95%, 99%, 99.9% at a confidence of 95%: 59, 299 and 2995 runs; the rule of three gives 60, 300 and 3000.
2. The count against the fraction. After 1000 runs with 0, 5 and 10 failures the bound on the success probability is 99.70%, 98.95% and 98.31%.
3. The campaign. The simulated campaign reaches 98.00% after 100 runs, 98.33% after 300, 99.00% after 1000, 99.30% after 3000 and 99.06% after 10 000, with lower bounds of 93.8%, 96.5%, 98.3%, 99.0% and 98.9%: at none of these sizes does it demonstrate 99%.
4. Why. Integrating the model gives a true success probability of 99.17%, 0.17 points above the requirement; the 10 000-run estimate (94 failures against 82.6 expected) was a little unlucky. At that margin a campaign of 10 000 runs demonstrates the requirement only 55% of the time, and one of 30 000 runs 93%.
5. The tail. The radius containing 99% of the 10 000 flights is 2.97 m (2.88 m without faults); for 99.9% it is 3.90 m against 3.53 m. A Gaussian fit to the bulk underestimates the far tail by 10%.
6. Failures allowed. At 95% confidence the requirement is demonstrated by at most 0 failures in 300 runs, 4 in 1000, 20 in 3000 and 83 in 10 000.
7. Repeated looks. A vehicle at 98.9% in 2000 simulated campaigns of up to 3000 runs: accepted by a single look at 3000 runs in 1.2% of them; by stopping at the first of 15 looks that passes, in 6.3%.

import math
import random

R = 3.0                          # requirement: land within 3 m ...
P_REQ = 0.99                     # ... in at least 99% of flights
CONF = 0.95                      # the confidence the campaign must reach

def binom_cdf(k, n, p):
    """P(X <= k) for X ~ Binomial(n, p), summed in logs so that large n is safe."""
    if p <= 0: return 1.0
    if p >= 1: return 1.0 if k >= n else 0.0
    s = 0.0
    for i in range(k + 1):
        s += math.exp(math.lgamma(n + 1) - math.lgamma(i + 1) - math.lgamma(n - i + 1) + i * math.log(p) + (n - i) * math.log1p(-p))
    return min(1.0, s)

def bisect(f, lo, hi):
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        if f(mid) > 0: lo = mid
        else: hi = mid
    return 0.5 * (lo + hi)

def success_lower(k_fail, n, conf=CONF):
    """One-sided Clopper-Pearson lower bound on the success probability, given k_fail failures in n runs."""
    # the largest failure probability q for which seeing k_fail or fewer failures still has probability >= 1 - conf
    q = bisect(lambda q: binom_cdf(k_fail, n, q) - (1 - conf), 0.0, 1.0)
    return 1 - q

# how many failure-free runs demonstrate a probability at a confidence
print("failure-free runs needed to demonstrate a probability with 95% confidence:")
for p in (0.95, 0.99, 0.999):
    n = math.ceil(math.log(1 - CONF) / math.log(p))
    print(f"  {p:.3f}: {n} runs  (rule of three: 3/(1-p) = {3 / (1 - p):.0f})")

# the lower bound after zero, one and two failures in 1000 runs
print("lower bound on the success probability after 1000 runs:")
for k in (0, 5, 10):
    print(f"  {k:2d} failures: estimate {1 - k / 1000:.3f}, lower bound {success_lower(k, 1000):.4f}")

# a dispersed Hopper campaign: linear-loop miss per axis (sigma 0.95 m), plus a beacon fault in 0.4% of flights that adds a 2.5 m offset
SIG, P_FAULT, OFFSET = 0.95, 0.004, 2.5
rng = random.Random(2024)
def flight():
    x, y = rng.gauss(0, SIG), rng.gauss(0, SIG)
    if rng.random() < P_FAULT:
        a = rng.uniform(0, 2 * math.pi)
        x += OFFSET * math.cos(a); y += OFFSET * math.sin(a)
    return math.hypot(x, y)

p_gauss = 1 - math.exp(-R * R / (2 * SIG * SIG))
print(f"the Gaussian linear-loop prediction without faults: {p_gauss:.4f} inside {R:.0f} m")
misses = []
print("   runs  failures  estimate  95% lower bound  decision")
marks = (100, 300, 1000, 3000, 10000)
for n in range(1, 10001):
    misses.append(flight())
    if n in marks:
        k = sum(1 for d in misses if d > R)
        lb = success_lower(k, n)
        print(f"{n:7d}  {k:8d}  {1 - k / n:8.4f}  {lb:15.4f}  {'requirement demonstrated' if lb >= P_REQ else 'not demonstrated'}")

def inside(cx, cy, N=300):
    """Probability that a 2-D Gaussian of spread SIG centred at (cx, cy) falls within R of the origin, by midpoint integration."""
    h = 2 * R / N; tot = 0.0
    for i in range(N):
        x = -R + (i + 0.5) * h
        for j in range(N):
            y = -R + (j + 0.5) * h
            if x * x + y * y <= R * R:
                tot += math.exp(-((x - cx) ** 2 + (y - cy) ** 2) / (2 * SIG * SIG)) * h * h / (2 * math.pi * SIG * SIG)
    return tot
p_true = (1 - P_FAULT) * inside(0, 0) + P_FAULT * inside(OFFSET, 0)
print(f"the true success probability with the faults, by integration: {p_true:.4f}; expected failures in 10000 runs {10000 * (1 - p_true):.1f}, observed {sum(1 for d in misses if d > R)}")
srt = sorted(misses)
r99 = srt[int(0.99 * len(srt))]
r999 = srt[int(0.999 * len(srt))]
print(f"radius containing 99% of 10000 flights: {r99:.2f} m (Gaussian without faults {SIG * math.sqrt(-2 * math.log(0.01)):.2f} m)")
print(f"radius containing 99.9%: {r999:.2f} m (Gaussian without faults {SIG * math.sqrt(-2 * math.log(0.001)):.2f} m)")

# the largest number of failures that still demonstrates the requirement, for a campaign of n runs
def k_allowed(n):
    k = -1
    while binom_cdf(k + 1, n, 1 - P_REQ) <= 1 - CONF: k += 1
    return k
print("failures allowed for the requirement to be demonstrated at 95% confidence:")
for n in (300, 1000, 3000, 10000):
    print(f"  {n:6d} runs: at most {k_allowed(n)} failures ({k_allowed(n) / n:.4f} of the runs)")

# looking again and again and stopping at the first success is not the same test: a vehicle that truly succeeds 98.9% of the time
looks = list(range(200, 3001, 200))
allowed = {n: k_allowed(n) for n in looks}
rng = random.Random(99)
trials, accepted_once, accepted_final = 2000, 0, 0
for _ in range(trials):
    fails = 0; n = 0; hit = False
    for nl in looks:
        while n < nl:
            n += 1
            if rng.random() < 0.011: fails += 1
        if fails <= allowed[nl]: hit = True
        last = fails <= allowed[nl]
    accepted_once += hit; accepted_final += last
print(f"a vehicle that truly meets only 98.9%, {trials} campaigns of up to 3000 runs: accepted by a single look at 3000 runs {accepted_final / trials:.1%}, by stopping at the first of {len(looks)} looks that passes {accepted_once / trials:.1%} (the confidence promises 5% at most)")

for n in (10000, 30000):
    kk = k_allowed(n)
    print(f"a campaign of {n} runs of the true {p_true:.4f} vehicle demonstrates the requirement with probability {binom_cdf(kk, n, 1 - p_true):.2f}")

failure-free runs needed to demonstrate a probability with 95% confidence:
  0.950: 59 runs  (rule of three: 3/(1-p) = 60)
  0.990: 299 runs  (rule of three: 3/(1-p) = 300)
  0.999: 2995 runs  (rule of three: 3/(1-p) = 3000)
lower bound on the success probability after 1000 runs:
   0 failures: estimate 1.000, lower bound 0.9970
   5 failures: estimate 0.995, lower bound 0.9895
  10 failures: estimate 0.990, lower bound 0.9831
the Gaussian linear-loop prediction without faults: 0.9932 inside 3 m
   runs  failures  estimate  95% lower bound  decision
    100         2    0.9800           0.9384  not demonstrated
    300         5    0.9833           0.9653  not demonstrated
   1000        10    0.9900           0.9831  not demonstrated
   3000        21    0.9930           0.9899  not demonstrated
  10000        94    0.9906           0.9889  not demonstrated
the true success probability with the faults, by integration: 0.9917; expected failures in 10000 runs 82.6, observed 94
radius containing 99% of 10000 flights: 2.97 m (Gaussian without faults 2.88 m)
radius containing 99.9%: 3.90 m (Gaussian without faults 3.53 m)
failures allowed for the requirement to be demonstrated at 95% confidence:
     300 runs: at most 0 failures (0.0000 of the runs)
    1000 runs: at most 4 failures (0.0040 of the runs)
    3000 runs: at most 20 failures (0.0067 of the runs)
   10000 runs: at most 83 failures (0.0083 of the runs)
a vehicle that truly meets only 98.9%, 2000 campaigns of up to 3000 runs: accepted by a single look at 3000 runs 1.2%, by stopping at the first of 15 looks that passes 6.3% (the confidence promises 5% at most)
a campaign of 10000 runs of the true 0.9917 vehicle demonstrates the requirement with probability 0.55
a campaign of 30000 runs of the true 0.9917 vehicle demonstrates the requirement with probability 0.93

In robotics

These rules decide what a campaign report may claim.

• Acceptance campaigns. An acceptance criterion is written as a bound ("99% at 95% confidence"), the number of runs is fixed in advance, and the table of allowed failures is part of the plan. The same table serves as a go or no-go for a new software version: 4 failures in 1000 runs passes, 5 does not.
• Margin determines cost. A design with 99.9% true success is demonstrated with a few thousand runs; one at 99.17% needs tens of thousands. Improving the design is cheaper than extending the campaign.
• Reading the tail. The fault flights make the far tail heavier than the Gaussian prediction: a Gaussian fit to 10 000 flights predicts 3.53 m for 99.9%, the run shows 3.90 m. A requirement stated on a tail, such as 99.9% containment, must be read from the runs and not from a fitted curve, and needs many more of them.
• Seeds and reproducibility. Each run is made reproducible by its seed, so a failing flight can be replayed with logging and examined: the first use of a campaign is the list of failures, and the second the number. The seeds of the failures are archived with the report.
• Combining with the analysis. The covariance analysis of M12-T1: Error Budgets and Covariance Analysis predicts which inputs to disperse widely (those with the largest variance share) and gives an expected result for the linear case to check the simulation against.

The statistics of seeds and runs for learned policies are in Robot Learning M11-T2: Seeds, Variance and Statistics, and the planning of the flight tests that follow in UAV M8-T5: Flight-test planning.

Common misconceptions

Key takeaways

1. A campaign demonstrates a lower confidence bound on the success probability, not the observed fraction.
2. Failure-free runs needed: $\ln\alpha/\ln p_0$, about 3/(1−p₀) at 95%: 299 for 99%, 2995 for 99.9%.
3. The bound for any failure count comes from the exact binomial sum; 5 failures in 1000 runs gives 98.95%, below 99%.
4. The smaller the margin between the true probability and the requirement, the more runs; at a margin of 0.17 points, 10 000 runs decide only about half the time.
5. Fix the number of runs and looks in advance: stopping at the first pass took the false-acceptance rate from 1.2% to 6.3%.
6. Read the tail from the runs, not from a Gaussian fit, and archive the seeds of the failures.

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.

• C. J. Clopper and E. S. Pearson, "The use of confidence or fiducial limits illustrated in the case of the binomial", Biometrika 26(4), 1934, the source of the exact bound. Paid (Oxford University Press).
• L. D. Brown, T. T. Cai and A. DasGupta, "Interval estimation for a binomial proportion", Statistical Science 16(2), 2001, a comparison of the bounds. Free to read from Project Euclid.
• W. Feller, An Introduction to Probability Theory and Its Applications, vol. 1, Wiley, 1968, for the binomial law. Paid.
• MIT 18.05 Introduction to Probability and Statistics (OpenCourseWare) , for confidence intervals and the binomial. Free; licence CC BY-NC-SA 4.0.

Knowledge check

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

What’s next

Related content

Nothing below is required to finish this lesson.

• Before
  M12-T1: Error Budgets and Covariance Analysis

  Allocate a requirement over error sources, propagate their covariance through a linearized loop, and find the source that dominates.

• Lab
  M12-L1: How Many Runs Are Enough?

  A requirement says 99% of landings within 3 m. Build the error budget, run a dispersion campaign, and state with what confidence the requirement is met.

• Lab
  M12-L2: Capstone: Land Hopper

  Write the navigation filter, the guidance law and the controller of Hopper, and pass an acceptance campaign of dispersed flights with wind, sensor errors, mass uncertainty and one beacon fault.

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

• Probability M1-T3: Simulating Randomness: Sampling and Monte Carlo
• Probability M4-T4: Propagating Uncertainty Through Transformations
• Linear Algebra M9-T4: Covariance Propagation
• UAV M8-T9: Fault detection, isolation and recovery
• UAV M8-T7: FMEA
• UAV M8-T3: Verification and validation
• UAV M8-T4: SIL, HIL and ground testing
• UAV M6-T14: Software-in-the-Loop simulation
• UAV M7-T7: Jitter and deadline monitoring
• UAV M8-T1: System requirements engineering
• UAV M8-T5: Flight-test planning
• Navigation M10-T2: Engineering Trade-offs and Safety Cases
• Robot Learning M11-T2: Seeds, Variance and Statistics

This course leaves the following to the pages that already teach them:

• Requirements engineering, FMEA and fault-tree analysis as methods
• SIL and HIL test benches, real-time scheduling and time synchronization
• Flight-test planning and incident investigation