Failure Distribution

Chapter 1: The Unnormalized Trick

Here is the key mathematical observation that makes everything in this chapter possible. We want to sample from p(τ | failure), but we don't know the normalizing constant p(failure). No problem. Define the unnormalized failure density:

This is just the trajectory probability p(τ) multiplied by the failure indicator. For non-failure trajectories, it's zero. For failure trajectories, it equals p(τ). The key identity:

The normalizing constant Z is literally p(failure) — the total probability mass assigned to failure trajectories. We don't know Z, but here is the beautiful part: MCMC doesn't need Z. Metropolis-Hastings only uses ratios of the target density, and Z cancels in the ratio:

As long as both τ and τ' are failure trajectories, this ratio is simply p(τ') / p(τ) — the likelihood ratio, which we can always compute.

Let's ground this with a simple example. Suppose the system has a 1D disturbance x ~ N(0, 1), and failure occurs when |x| > 2. The failure distribution is:

This is a standard normal truncated to the tails beyond ±2. The unnormalized version is just the normal density times the indicator — two bumps, one on each side of zero, decaying outward. We don't need to know how much probability is in the tails to sample from this shape.

Chapter 2: Rejection Sampling

Before reaching for MCMC, let's try the simplest possible approach: rejection sampling. The idea is charmingly straightforward. Sample trajectories from the nominal distribution p(τ), and keep only those that fail. Throw away the rest. The survivors are exact samples from p(τ | failure).

Think of throwing darts at a dartboard. The board is the full trajectory space. The bullseye is the failure region. If you throw darts uniformly (well, proportional to p(τ)), the ones that hit the bullseye are distributed exactly as p(τ | failure). Simple, exact, elegant.

More generally, rejection sampling works with any proposal distribution q(τ) that covers the support of the target. We need a constant c such that c · q(τ) ≥ p̄(τ | failure) everywhere. The algorithm:

The acceptance rate is 1/c. If q = p (the nominal distribution) and the target is p(τ | failure), then c = 1/p(failure). The acceptance rate equals the failure probability itself.

When the failure threshold is at 2.0, about 4.6% of standard normal samples exceed it — so roughly 1 in 22 proposals is accepted. Manageable. But push the threshold to 3.0 (0.27% acceptance) or 4.0 (0.006% acceptance) and you'll be waiting a long time. For rare failures, rejection sampling is hopelessly inefficient.

Chapter 3: When Rejection Fails

Rejection sampling breaks down for the same reason that naive Monte Carlo estimation of rare events breaks down: the acceptance rate is proportional to the failure probability. If failures happen once in a million trajectories, we need roughly a million proposals to get one accepted sample.

But the real killer is dimensionality. Consider an inverted pendulum with 100 timesteps of 3D disturbances. The trajectory space is 300-dimensional. Even if the marginal failure probability isn't astronomically small, the constant c that bounds the density ratio grows exponentially with dimension.

Why is c so large? Because the proposal q(τ) spreads its probability mass over the entire trajectory space — including the vast majority that doesn't lead to failure. The failure-conditional distribution concentrates its mass in a tiny region. The ratio p̄(τ | failure) / q(τ) spikes enormously at failure trajectories, forcing c to be enormous to bound it.

There is one more subtle problem. Even if we could magically improve the acceptance rate, rejection sampling produces independent samples. That sounds good, but it means we cannot use information from previous samples to guide future ones. MCMC trades independence for adaptivity: each sample informs the next, creating a chain that efficiently explores the target distribution.

Chapter 4: MCMC Basics

Rejection sampling tries to generate independent samples and mostly fails. Markov chain Monte Carlo (MCMC) takes a different philosophy: generate correlated samples that form a chain, where each sample is a small perturbation of the previous one. The chain is designed so that its long-run behavior matches the target distribution.

Think of a drunk person wandering around a city. If the bars are in the failure region and the wanderer has a slight tendency to walk toward bars, then over time, the wanderer's position traces out the distribution of bar locations. Individual positions are correlated (each step is near the previous), but the collection of all positions visited gives us a map of where bars are.

A Markov chain is defined by a transition kernel g(τ' | τ) that says: given the current sample τ, where does the chain go next? The chain has the correct stationary distribution if running it long enough produces samples from p(τ | failure), regardless of where it started.

The sufficient condition for correct stationary distribution is detailed balance:

where T is the overall transition kernel (proposal + accept/reject). This says: the flow of probability from τ to τ' equals the flow from τ' to τ. If this holds for all pairs, the distribution is in equilibrium.

Key concepts to watch for in the simulation:

Chapter 5: Metropolis-Hastings

How do we build an MCMC chain with the right stationary distribution? The Metropolis-Hastings algorithm gives a general recipe. Start with any proposal kernel g(τ' | τ), and add an accept/reject step that corrects for the proposal's bias.

Given the current sample τ, propose τ' ~ g(τ' | τ). Accept with probability:

Let's unpack this. The ratio p̄(τ') / p̄(τ) measures whether the proposal is in a higher-density region of the target. The ratio g(τ | τ') / g(τ' | τ) corrects for asymmetry in the proposal. If the proposal is more likely to jump from τ to τ' than vice versa, we need to accept less often to compensate.

For a symmetric proposal (like a Gaussian random walk where g(τ' | τ) = g(τ | τ')), the correction term cancels, and the acceptance probability simplifies to:

With a symmetric random walk proposal, this becomes: accept if the proposed trajectory is more probable than the current one. If the proposal is less probable, accept with a probability equal to the probability ratio. This means we mostly move "uphill" toward high-density regions but occasionally accept "downhill" moves, which prevents the chain from getting permanently stuck at a single mode.

What if the proposed trajectory is not a failure? Then 1{fail(τ')} = 0, so p̄(τ' | fail) = 0, and the acceptance probability is zero. The chain never leaves the failure region once it enters. This is correct — we're sampling the conditional distribution given failure.

pseudocode
# Metropolis-Hastings for failure distribution
Input: initial failure trajectory τ0, proposal kernel g, num_samples
chain = [τ0]

for i = 1 to num_samples:
  τ' ~ g(τ' | chain[i-1])
  if is_failure(τ'):
    α = min(1, p(τ') / p(chain[i-1]))  # symmetric g
  else:
    α = 0  # never accept non-failure
  if random() < α:
    chain.append(τ')  # accept
  else:
    chain.append(chain[i-1])  # reject: repeat current

return chain[burn_in:]  # discard burn-in

Chapter 6: The Smoothing Trick

There is a nasty problem with the MH algorithm as described. The failure indicator 1{fail(τ)} is a hard boundary. A trajectory either fails or it doesn't — there is no middle ground. If the current trajectory is barely a failure (robustness ρ = −0.001) and the proposal perturbs it slightly into non-failure territory (ρ = +0.001), the proposal is instantly rejected. The chain is trapped on one side of the failure boundary.

This is catastrophic when there are multiple failure modes separated by non-failure regions. The chain enters one failure mode and can never cross the safe region to reach the other. It's like being on an island — the MCMC walker can explore every beach on this island, but it can never swim to the next island because the ocean (non-failure region) has zero density.

The fix is elegant: replace the hard indicator with a smooth approximation. Instead of 1{fail(τ)}, use a Gaussian kernel centered at the failure boundary:

Here Δ(τ) is the robustness — positive means safe, negative means failure, zero is the boundary. The Gaussian kernel gives high weight to trajectories near the failure boundary (|Δ| small) and decaying weight to trajectories farther away on either side. Crucially, it assigns nonzero weight to safe trajectories near the boundary. This lets the chain cross through safe regions to reach other failure modes.

The parameter ε controls the smoothing width. Small ε approximates the hard indicator closely (good fidelity, poor mixing). Large ε smooths aggressively (poor fidelity, good mixing). The practical recipe: start with large ε, let the chain explore all modes, then gradually decrease ε toward zero to sharpen the distribution.

The effect is dramatic. Set ε = 0 (no smoothing) and run 200 steps — the chain stays entirely in one mode. Now set ε = 0.5 and reset — the chain crosses between modes freely. The trace plot tells the story: without smoothing, a flat line in one mode; with smoothing, the trace jumps between modes.

Chapter 7: Multimodal Challenges

The smoothing trick enables mode crossing, but multimodal failure distributions remain one of the hardest problems in computational validation. Let's understand why, and what strategies help beyond smoothing.

Consider an autonomous vehicle. One failure mode might be "rear-ended by a fast-following vehicle." Another might be "sudden pedestrian crossing." A third might be "sensor occlusion in heavy rain." These failure modes occupy completely different regions of the disturbance space. The failure distribution has three (or more) well-separated peaks.

Several strategies address this:

Parallel tempering is the most principled approach. Run K chains in parallel, each targeting the smoothed distribution with a different ε. The chain with the largest ε moves freely across the whole space. Periodically, propose swaps between adjacent chains. If the high-ε chain has discovered a new mode, the swap propagates this discovery down to the low-ε chains. The Metropolis criterion for swaps ensures detailed balance is maintained.

In practice, the most common approach for validation is the simplest: run the falsification methods from Chapter 5 (MCTS, RL adversary) to find diverse failure modes, then initialize an MCMC chain in each mode and explore locally. This hybrid approach uses each tool for what it does best — search algorithms for discovery, MCMC for characterization.

Chapter 8: Gradient-Based MCMC

The random walk MH algorithm proposes steps in random directions. In high dimensions, most random directions are "sideways" — they neither increase nor decrease the target density. The chain shuffles along the level set instead of moving efficiently toward high-density regions.

If we can compute the gradient of the log-density, we can bias proposals toward high-density directions. This is the idea behind gradient-based MCMC.

The Metropolis-Adjusted Langevin Algorithm (MALA) uses a proposal that follows the gradient:

The first term is a gradient step: the proposal drifts toward higher density. The second term is random noise for exploration. The step size ε balances the two. Because the proposal is no longer symmetric (it's biased toward the gradient direction), MH must correct with the full acceptance ratio including the proposal asymmetry.

Hamiltonian Monte Carlo (HMC) goes further. Instead of one gradient step, it simulates a physical system: the current position τ gets an initial momentum, then follows Hamiltonian dynamics (gradient + momentum) for L steps. The result is a proposal that can be far from the current point while maintaining high acceptance probability. HMC explores high-dimensional distributions dramatically faster than random walk MH.

The No-U-Turn Sampler (NUTS) automates HMC's tuning parameters. It simulates the trajectory forward and backward in time, stopping when the trajectory starts to double back on itself (the "U-turn"). NUTS is the default sampler in modern probabilistic programming systems (Stan, Turing.jl, NumPyro) and is the current gold standard for gradient-based MCMC.

Chapter 9: Probabilistic Programming

Everything in this chapter — MCMC, smoothing, gradients — can be implemented by hand. But doing so for complex systems is tedious and error-prone. Probabilistic programming languages (PPLs) automate the heavy lifting. You write the model, and the PPL handles the inference.

The key idea: encode the validation problem as a probabilistic program. The disturbances are random variables. The simulator is a deterministic function of those variables. The failure condition is an observation. The PPL's inference engine (typically NUTS) samples from the posterior — which is exactly the failure distribution.

julia
using Turing

@model function failure_distribution(simulator, ε)
    # Prior: disturbances from nominal distribution
    x ~ MvNormal(zeros(300), I)

    # Run simulator (deterministic given x)
    τ = simulator(x)

    # Compute robustness
    Δ = robustness(τ)

    # Smoothed failure indicator via @addlogprob!
    Turing.@addlogprob!(-Δ^2 / (2 * ε^2))
end

# Sample from failure distribution
model = failure_distribution(my_simulator, 0.1)
chain = sample(model, NUTS(), 1000)

The elegance is striking. The entire failure distribution sampling problem — disturbance model, simulator, robustness metric, smoothing, MCMC inference — fits in about 10 lines of code. The PPL handles all the machinery: gradient computation (via automatic differentiation), NUTS tuning, burn-in, diagnostics, and sample storage.

The requirement is that the simulator must be differentiable (for NUTS) or at least evaluable (for random walk MH). Most modern PPLs support automatic differentiation through arbitrary Julia/Python code, including array operations, conditionals, and loops. If your simulator is written in a supported language, it "just works."

What comes next: Sampling the failure distribution tells us what failures look like. Chapter 7 asks the remaining question: how likely is failure? This requires estimating p(failure) — the normalizing constant Z that we carefully avoided throughout this chapter. Importance sampling, adaptive methods, and sequential Monte Carlo provide the answer.