Importance Sampling for Rare Events

PyApprox Tutorial Library

When Monte Carlo needs millions of samples to see a single failure, importance sampling deliberately targets the tail and corrects with likelihood ratios.

TipDownload Notebook

Download as Jupyter Notebook

Learning Objectives

After completing this tutorial, you will be able to:

• Explain why Monte Carlo is inefficient for small failure probabilities and quantify the cost
• Formulate a failure probability as an integral and identify the variance bottleneck
• Derive the importance sampling estimator and explain the likelihood ratio correction
• Demonstrate empirically that importance sampling reduces variance by orders of magnitude for rare events
• Recognize the risk of weight degeneracy from a poor biasing distribution

Prerequisites

Complete Monte Carlo Sampling and Estimator Accuracy and MSE before this tutorial.

Why Rare Events Break Monte Carlo

Estimating the probability of a rare but critical event — a bridge deflecting past its safety limit, a circuit failing under thermal stress — is one of the hardest problems in uncertainty quantification. The difficulty is simple to state: if failures happen one time in ten thousand, then a Monte Carlo experiment with a few hundred samples will almost certainly see zero failures. The estimate \(\hat{p} = 0\) is technically unbiased (averaged over infinitely many experiments), but any single experiment is useless.

The QMC tutorial showed that structured point sets accelerate convergence for smooth integrands. But failure probabilities involve the indicator function \(\mathbb{1}[f > c]\), which is discontinuous — exactly the case where QMC’s advantage disappears.

We need a fundamentally different idea. Before deriving anything, let’s see the problem.

Seeing the Problem: Where Do the Samples Land?

Consider estimating \(p = P(x > 2)\) where \(x \sim \mathcal{N}(0,1)\). The true probability is about 2.3%. Figure 1 shows 120 samples drawn from the standard normal, with the failure region \(x > 2\) shaded.

Figure 1: 120 samples from the standard normal \(\pi\). The failure region \(x > 2\) (shaded red) sits in the far tail where \(\pi\) has almost no mass. Only 1 sample (red tick mark) lands past the threshold — the other 119 contribute nothing to the failure probability estimate.

The picture tells the whole story: the target density \(\pi\) concentrates almost all its mass between \(-2\) and \(2\). The failure region \(x > 2\) sits in the far tail where \(\pi\) is negligible — samples simply don’t land there. Of 120 draws, only 1 exceeds the threshold. The other 119 contribute nothing: each returns \(\mathbb{1}[x_i > 2] = 0\).

MC Estimation: The Consequence

We just saw that a single MC experiment with 120 samples produced only 1 failure. What happens if we repeat this experiment hundreds of times? Figure 2 runs 500 independent experiments at two sample sizes.

Figure 2: MC estimates of \(P(x > 2)\) from 500 independent experiments. Left: \(M = 500\) — many experiments see zero failures, producing \(\hat{p} = 0\) (the tall spike at the origin). Right: \(M = 5{,}000\) — fewer zeros, but estimates still scatter widely around the true value (red dashed line).

The spike at zero in the left panel is exactly what Figure 1 predicted: experiments that drew hundreds of samples without a single one reaching the tail. The estimator is unbiased — averaged over all 500 experiments the mean is close to the true \(p\) — but any individual estimate is nearly useless.

We can now quantify the pain. The MC estimator \(\hat{p}_M = \frac{1}{M}\sum \mathbb{1}[x_i > c]\) is a sample mean of Bernoulli(\(p\)) random variables. Its variance is \(p(1-p)/M \approx p/M\) when \(p\) is small, so the relative standard error is:

\[ \frac{\text{SE}[\hat{p}_M]}{p} = \frac{1}{\sqrt{Mp}} \]

To achieve 10% relative error we need \(M \approx 100/p\). For \(p = 10^{-3}\) that is \(100{,}000\) samples; for \(p = 10^{-6}\) it is \(100\) million. This is the formula behind the spike at zero: with \(M = 500\) and \(p \approx 0.023\), we need roughly \(M \approx 4{,}000\) samples for 10% accuracy, so most 500-sample experiments miss the tail entirely.

What if we could sample from a different distribution — one that puts its mass right where failures occur?

The Key Idea: Sample Where Failures Happen

The fix is to sample from a different distribution \(q\) that places mass where failures occur, then correct for the fact that we sampled from the wrong distribution. Figure 3 shows both steps.

Figure 3: The importance sampling fix, continuing the 1D example from Figure 1. (Left) A shifted proposal \(q\) (blue, centered at 2.5) places most of its mass past the threshold. Of 120 samples, the majority now land in the failure region — compare to the single hit in Figure 1. (Right) But we can’t just count them: they were drawn from \(q\), not \(\pi\). Each sample is placed on the target density \(\pi\) and sized by its importance weight \(w_i = \pi(x_i) / q(x_i)\). Samples near the threshold carry the largest weights (dark red, large markers) because they are relatively most probable under \(\pi\). Samples deep in the tail, where \(\pi\) is negligible, are downweighted (yellow, small markers).

The left panel is the “aha”: by shifting the sampling distribution to \(q\), the majority of samples now exceed the threshold. But these samples were drawn from \(q\), not \(\pi\) — we can’t simply count them. The right panel shows the correction: each sample is placed on the target density \(\pi\) and scaled by its importance weight \(w_i = \pi(x_i)/q(x_i)\). Samples near the threshold are the most valuable — they are relatively probable under \(\pi\) and carry the largest weights. Samples deep in the tail, where \(q \gg \pi\), contribute almost nothing after reweighting.

This is the entire importance sampling mechanism: sample where it matters, then correct for the bias.

Formalizing the Idea

The importance sampling estimator

The figure above showed the mechanism; now we write it down. The failure probability is an integral:

\[ p = P\bigl(f(\boldsymbol{\theta}) > c\bigr) = \int \mathbb{1}\bigl[f(\boldsymbol{\theta}) > c\bigr]\, \pi(\boldsymbol{\theta})\, d\boldsymbol{\theta} \]

Multiply and divide by the proposal density \(q\):

\[ p = \int \mathbb{1}[f(\boldsymbol{\theta}) > c]\, \frac{\pi(\boldsymbol{\theta})}{q(\boldsymbol{\theta})}\, q(\boldsymbol{\theta})\, d\boldsymbol{\theta} \]

Now sample from \(q\) instead of \(\pi\):

\[ \hat{p}_M^{\text{IS}} = \frac{1}{M} \sum_{i=1}^{M} \mathbb{1}[f(\boldsymbol{\theta}^{(i)}) > c]\, \frac{\pi(\boldsymbol{\theta}^{(i)})}{q(\boldsymbol{\theta}^{(i)})}, \qquad \boldsymbol{\theta}^{(i)} \sim q \]

The ratio \(w(\boldsymbol{\theta}) = \pi(\boldsymbol{\theta}) / q(\boldsymbol{\theta})\) is the likelihood ratio or importance weight. Each term maps directly to something in Figure 3:

• \(\boldsymbol{\theta}^{(i)} \sim q\) — left panel: samples drawn from the shifted proposal
• \(\mathbb{1}[f(\boldsymbol{\theta}^{(i)}) > c]\) — the red tick marks: which samples exceed the threshold
• \(w_i = \pi / q\) — right panel: the marker sizes that correct for the distributional shift

The estimator is unbiased for any valid \(q\) (any distribution positive wherever \(\pi \cdot \mathbb{1}[f > c] > 0\)). The choice of \(q\) affects only the variance — and a good choice can reduce it by orders of magnitude.

Setup: The Beam Model with a Failure Threshold

We now apply importance sampling to a realistic engineering model. We use the analytical cantilever beam with \(\mathrm{Beta}(2, 5)\) priors on the two Young’s moduli, the same model as the QMC tutorial.

from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox.benchmarks.instances.analytic.cantilever_beam import (
    cantilever_beam_1d_analytical,
)

bkd = NumpyBkd()

benchmark = cantilever_beam_1d_analytical(bkd)
model = benchmark.function()
variable = benchmark.prior()

print(f"Number of uncertain parameters: {model.nvars()}")
print(f"Number of QoIs: {model.nqoi()}")

Number of uncertain parameters: 2
Number of QoIs: 3

We define “failure” as tip deflection exceeding a threshold \(c\). We choose \(c\) so that the true failure probability is small but still computable by quadrature — around \(p \approx 0.01\) to \(0.001\), depending on the model’s output range.

# Choose a threshold in the upper tail of the deflection distribution
# using quadrature-weighted CDF evaluation.
# Search for a threshold that gives true p ≈ 0.005
# by evaluating P(QoI > c) via quadrature at many candidate thresholds.
qoi_np = bkd.to_numpy(qoi_quad)
wts_np = bkd.to_numpy(quad_wts)
candidates = np.linspace(np.percentile(qoi_np, 95), np.max(qoi_np), 200)
probs = np.array([
    float(np.dot(wts_np, (qoi_np > c).astype(float))) for c in candidates
])

# Pick the threshold closest to the target failure probability
target_p = 0.005
threshold = float(candidates[np.argmin(np.abs(probs - target_p))])

# True failure probability via quadrature
indicator = bkd.asarray((qoi_np > threshold).astype(float))
true_p = bkd.to_float(bkd.sum(quad_wts * indicator))

print(f"Failure threshold c:         {threshold:.6f}")
print(f"True failure probability p:  {true_p:.6f}")
print(f"Samples needed for 10% relative SE (MC): {int(100 / true_p):,}")

Failure threshold c:         4.470276
True failure probability p:  0.004884
Samples needed for 10% relative SE (MC): 20,476

Visualizing the Failure Region

Figure 4 shows the beam model’s response surface in \((E_1, E_2)\) input space. The failure region — where the tip deflection exceeds the threshold — is a thin sliver in the corner of the input domain. Most samples drawn from the prior land far from this region and contribute nothing to the failure probability estimate.

Figure 4: The beam’s input space with tip deflection shown as a contour map. The failure region \(\delta_{\text{tip}} > c\) (red shading) occupies a small fraction of the prior’s support. Of \(M = 2{,}000\) prior samples, only a handful (red diamonds) land in the failure region — the rest (blue dots) contribute zero to the failure probability estimate.

Choosing a Biasing Distribution

For the beam model with bounded Beta priors on \((E_1, E_2)\), a practical strategy is to shift the sampling distribution toward the failure region. If we know approximately where failures occur — for example, from the response surface in Figure 4 — we can center \(q\) near the failure boundary.

A simple approach: use a truncated normal (or Beta with adjusted parameters) centered near the region of input space where \(f\) is close to the threshold. For this tutorial, we approximate the failure region’s center from the response surface and define \(q\) as a product of independent shifted Beta distributions.

# Identify the approximate center of the failure region
# using the prior samples that exceeded the threshold
fail_samples = samples_vis_np[:, fail_mask]
if fail_samples.shape[1] > 0:
    shift_center = np.mean(fail_samples, axis=1)
else:
    # Fallback: use the point with the highest QoI value
    shift_center = samples_vis_np[:, np.argmax(qoi_vis)]

print(f"Approximate failure region center: {shift_center}")

Approximate failure region center: [18161.47802503  4556.64103994]

# Build the biasing distribution q as a product of independent
# truncated normals centered near the failure region, truncated
# to the prior's support bounds.

# Extract prior support bounds from the marginals
marginal_list = variable.marginals()
prior_bounds = [m.bounds() for m in marginal_list]

# Define biasing distribution: truncated normal per dimension
# Spread chosen to cover the failure region without being too narrow
is_spread_factor = 0.15  # fraction of support width

def log_q_pdf(theta, shift_center, bounds, spread_factor):
    """Log-pdf of the biasing distribution (product of truncated normals)."""
    log_pdf = 0.0
    for d in range(len(shift_center)):
        lo, hi = bounds[d]
        sigma = spread_factor * (hi - lo)
        a, b = (lo - shift_center[d]) / sigma, (hi - shift_center[d]) / sigma
        log_pdf += stats.truncnorm.logpdf(
            theta[d], a, b, loc=shift_center[d], scale=sigma
        )
    return log_pdf

def sample_q(M, shift_center, bounds, spread_factor, rng):
    """Draw M samples from the biasing distribution."""
    nvars = len(shift_center)
    samples = np.empty((nvars, M))
    for d in range(nvars):
        lo, hi = bounds[d]
        sigma = spread_factor * (hi - lo)
        a, b = (lo - shift_center[d]) / sigma, (hi - shift_center[d]) / sigma
        samples[d] = stats.truncnorm.rvs(
            a, b, loc=shift_center[d], scale=sigma, size=M,
            random_state=rng,
        )
    return samples

def log_prior_pdf(theta, variable, bkd):
    """Log-pdf of the prior (IndependentJoint) evaluated at theta."""
    # variable.logpdf expects shape (nvars, nsamples) and returns (1, nsamples)
    log_pdf_arr = variable.logpdf(bkd.asarray(theta))
    return bkd.to_numpy(log_pdf_arr).ravel()

Visualizing the Biasing Distribution

Figure 5 compares prior sampling (left) with importance sampling (right). The IS samples cluster in the failure region; the marker sizes show the importance weights that correct for the distributional shift.

Figure 5: Left: \(M = 500\) samples from the prior \(\pi\) — few reach the failure region. Right: \(M = 500\) samples from the biasing distribution \(q\) (truncated normal centered near the failure boundary). Many more samples land in the failure region. Marker size is proportional to the importance weight \(w_i = \pi(\boldsymbol{\theta}^{(i)}) / q(\boldsymbol{\theta}^{(i)})\); samples far from the prior center receive smaller weights.

IS in Action: Variance Reduction

We now run the same 500-replicate experiment as before, but using importance sampling with \(M = 500\) samples per replicate.

def is_failure_estimates(model, variable, threshold, M, n_reps,
                         shift_center, bounds, spread_factor, bkd):
    """Run n_reps IS experiments, each with M samples."""
    estimates = np.empty(n_reps)
    for i in range(n_reps):
        rng = np.random.RandomState(i + 100000)
        s = sample_q(M, shift_center, bounds, spread_factor, rng)
        # Evaluate model using bkd arrays
        s_bkd = bkd.asarray(s)
        q_vals = bkd.to_numpy(model(s_bkd)[0])
        indicator = (q_vals > threshold).astype(float)

        # Importance weights: pi(theta) / q(theta)
        log_w = log_prior_pdf(s, variable, bkd) - log_q_pdf(
            s, shift_center, bounds, spread_factor
        )
        w = np.exp(log_w)
        estimates[i] = np.mean(indicator * w)
    return estimates

is_estimates_500 = is_failure_estimates(
    model, variable, threshold, 500, n_reps,
    shift_center, prior_bounds, is_spread_factor, bkd,
)

Figure 6 places the MC and IS histograms side by side.

Figure 6: Failure probability estimates from 500 experiments with \(M = 500\) samples each. Left: standard MC — wide spread with many zero estimates. Right: importance sampling with a shifted truncated normal — dramatically narrower, no zeros, centered on the true probability.

The IS histogram is dramatically narrower and has no zeros. The variance reduction factor quantifies the improvement: IS with \(M = 500\) is as precise as MC with many thousands of samples.

Convergence Comparison

Figure 7 shows the relative RMSE of both estimators as a function of \(M\). Both converge at rate \(O(M^{-1/2})\) — importance sampling does not change the rate. It reduces the constant, which on a log-log plot appears as a vertical shift between two parallel lines.

Figure 7: Relative RMSE of the failure probability estimate vs. sample size. Both MC (blue) and IS (green) converge at \(O(M^{-1/2})\) (parallel lines on the log-log plot), but the IS constant is orders of magnitude smaller. IS with \(M = 100\) can achieve the same accuracy as MC with \(M \approx 50{,}000\).

The Danger of Bad IS: Weight Degeneracy

Importance sampling is not a free lunch. If \(q\) is poorly chosen — shifted too far into the tail, or made too narrow — a small number of samples can acquire enormous importance weights while the rest contribute almost nothing. This is weight degeneracy: the effective sample size collapses even though \(M\) is large.

The effective sample size diagnoses this:

\[ \text{ESS} = \frac{\bigl(\sum_{i=1}^{M} w_i\bigr)^2}{\sum_{i=1}^{M} w_i^2} \]

When all weights are equal, \(\text{ESS} = M\). When one weight dominates, \(\text{ESS} \approx 1\). A ratio \(\text{ESS}/M \ll 1\) signals that the IS estimate is unreliable.

Figure 8 compares the distribution of log-importance-weights for three choices of the biasing distribution: a well-chosen shift, an overshifted \(q\), and an overly narrow \(q\).

Figure 8: Histograms of log importance weights \(\log w_i\) for three biasing distributions, each with \(M = 500\). Left: well-chosen \(q\) near the failure boundary — weights are moderate, ESS is high. Center: \(q\) shifted too far past the failure region — a few extreme weights dominate, ESS collapses. Right: \(q\) with too-small spread — similar degeneracy. The effective sample size (ESS) is shown in each panel.

The well-chosen \(q\) produces log-weights that are concentrated in a narrow range — the effective sample size is a healthy fraction of \(M\). The overshifted and too-narrow distributions produce log-weight histograms with heavy tails, indicating that a handful of extreme weights dominate the estimator.

Finding the Failure Region in Practice

The IS estimator above assumed we already knew approximately where failures occur. For the 2D beam model, we could read this from the response surface plot. In realistic problems with expensive models and high-dimensional input spaces, this is not straightforward.

NoteLooking ahead: surrogates and optimization

A common practical workflow is to build a cheap surrogate model (e.g., a polynomial chaos expansion or Gaussian process emulator) from a modest pilot sample, and then use the surrogate to identify the failure region. Specifically:

1. Optimization on the surrogate. Solve \(\min_{\boldsymbol{\theta}} -f_{\text{surr}}(\boldsymbol{\theta})\) subject to \(\boldsymbol{\theta}\) being in the prior’s support, or equivalently find the most probable failure point (MPP) by minimizing the distance from the prior mode subject to \(f_{\text{surr}}(\boldsymbol{\theta}) \ge c\). This is a standard reliability analysis technique (FORM/SORM) and locates the region of input space where failures are most likely under the prior.

2. Hessian-based proposal covariance. The Hessian of \(\log \pi(\boldsymbol{\theta})\) (or of the limit-state function \(g = f - c\)) at the MPP provides curvature information about the failure boundary. Its inverse gives a natural covariance matrix for the biasing distribution \(q\) — a multivariate normal centered at the MPP with covariance proportional to the inverse Hessian. This adapts the shape and orientation of \(q\) to the local geometry of the failure surface, without requiring manual tuning.

3. Iterative refinement. Run IS with the surrogate-based \(q\), use the resulting weighted samples to update the surrogate in the failure region, and re-optimize. This focuses surrogate accuracy where it matters most for the failure probability estimate.

These techniques make IS practical for expensive, high-dimensional models where visual inspection of the failure region is impossible. They are the subject of a later tutorial.

Exercises

1. Vary the threshold \(c\) to produce true failure probabilities ranging from \(10^{-2}\) to \(10^{-4}\). How does the IS variance reduction factor (ratio of MC variance to IS variance at fixed \(M\)) change as \(p\) decreases? Does IS become more or less valuable for rarer events?

2. For the beam model, use a pilot MC run of \(M = 200\) to estimate the center of the failure region (mean of the failing samples). Use this as shift_center for IS. Compare the resulting variance reduction to using an arbitrary shift. How sensitive is the result to the accuracy of the center estimate?

3. Implement the effective sample size \(\text{ESS} = (\sum w_i)^2 / \sum w_i^2\) and monitor it as you vary the is_spread_factor from 0.05 to 0.5. Plot \(\text{ESS}/M\) vs. spread factor. Where does ESS collapse? What is the practical “sweet spot”?

4. (Challenge) Combine importance sampling with QMC: use a Sobol sequence in \([0,1]^2\) and transform it through \(q\)’s inverse CDF (instead of the prior’s). Compare the variance of this QMC-IS estimator to both plain IS and plain MC. Does the QMC structure help even though the integrand contains a discontinuous indicator?

Key Takeaways

• MC needs \(M \approx 100/p\) samples for 10% relative error on a failure probability \(p\) — prohibitive for rare events
• QMC cannot help much: the indicator integrand is discontinuous, so the Koksma–Hlawka advantage disappears
• Importance sampling reweights samples from a biasing distribution \(q\) that concentrates mass in the failure region, correcting via the likelihood ratio \(\pi/q\)
• IS does not change the convergence rate (\(1/\sqrt{M}\)) but can reduce the variance constant by orders of magnitude
• Poor choice of \(q\) leads to weight degeneracy — a few extreme weights dominate and the estimate becomes unreliable; the effective sample size ESS diagnoses this
• In practice, surrogates and optimization (FORM/SORM, Hessian-based proposals) provide systematic ways to construct good biasing distributions for expensive, high-dimensional models

Next Steps

Continue with:

• Control Variates — Using correlated auxiliary quantities to reduce estimator variance