Scientific Computing · Simulation

Rare-Event Estimation with Splitting and the Cross-Entropy Shift

When the event of interest has probability a million to one, plain sampling sees nothing, and the standard tools bias the proposal or multiply the paths to reach the tail in bounded steps.

An event of probability $10^{-6}$ needs on the order of a million replications before a single occurrence is seen, and the standard error of the estimate is then the probability itself. Rare-event simulation replaces the hopeless brute-force count with constructions that reach the tail deliberately: importance splitting multiplies the paths that are already rare, and the cross-entropy method learns a proposal that makes the event typical. Both attack the same failure of the square-root law, and both have failure modes of their own that are as instructive as the methods.

The problem is the coefficient of variation, not the mean

For a Bernoulli estimate of $p = P(X \in A)$, the relative error after $n$ draws is $\sqrt{(1-p)/(np)} \approx 1/\sqrt{np}$, so the number of draws needed for a given relative accuracy scales as $1/p$. At $p = 10^{-6}$ a ten-percent relative error needs around $10^8$ draws, and the estimate is zero for the overwhelming majority of runs that see no event at all. The estimator is not just slow; it is degenerate, because the event is invisible until the rare path is actually generated.

The fix has to reach the event, and reaching it requires changing the sampling distribution so the event becomes probable. This is the same likelihood-ratio idea as importance sampling, but the target is a probability rather than a mean, and the optimal proposal is known in closed form: sample from the distribution conditioned on the event. The entire field is about approximating that conditional distribution without being able to sample from it directly.

Splitting multiplies the survivors

Importance splitting works on a nested sequence of rare events. Choose levels $A_0 \supset A_1 \supset \dots \supset A_m = A$ where $A_0$ is certain, and each step is only moderately rare, so that $P(A_{i+1} \mid A_i)$ is around, say, $1/10$. The estimator runs paths, and whenever a path reaches level $i$, it is cloned into $R$ copies that continue independently. The probability estimate is the product of the observed level-crossing ratios, and the cloning keeps the population alive through levels that plain sampling would almost never reach.

The estimator is unbiased if the branching is done correctly. Each path reaching a level spawns a fixed number of offspring, and the final estimate combines the fraction of paths reaching each level with the offspring counts. The variance is reduced because the rare part of the trajectory, the climb from $A_{m-1}$ to $A_m$, is sampled many times instead of once per million runs.

The difficulty is the choice of levels. If the levels are too far apart, each crossing is itself rare and the estimator is noisy at that step; if they are too close, the number of levels is large and the bookkeeping of the branching tree becomes the cost. The standard guidance is to choose levels so that each conditional crossing probability is roughly constant, which makes every step contribute equal information and keeps the branching population from either exploding or dying out.

import random

def split_estimate(level_func, R, n_levels, trials_per_level):
    """Estimate P(end in target) by multilevel splitting; level_func(k)
    reports the highest level reached by a path started from level k."""
    n = trials_per_level
    for k in range(n_levels):
        reached = 0
        for _ in range(n):
            if level_func(k) > k:
                reached += 1
        # Each survivor spawns R offspring for the next level.
        n = reached * R
        if n == 0:
            return 0.0
    return n / (trials_per_level * R ** n_levels)

The arithmetic exposes the bookkeeping. The final probability is the product of the per-level crossing ratios, and each ratio is estimated by the surviving fraction times the branching factor. If a level is too hard, reached is zero, the population collapses, and the estimate returns zero even though the true probability is positive. The branching factor and the level spacing have to be balanced so that the population stays near a target size.

The cross-entropy method learns the proposal

The cross-entropy method attacks the same problem by fitting a parametric family to the event. Starting from an initial proposal $g_0$, it draws a batch, keeps the small fraction $\rho$ with the largest scores (the elite), and refits the parameters to the elite by maximizing their likelihood. The process repeats, and the fitted proposal migrates toward the conditional distribution of the event, so that after a few iterations the event is no longer rare under the proposal and the importance-sampling weights are moderate.

The name comes from the optimization that drives the refit. Maximizing the likelihood of the elite samples under the new parameters is equivalent to minimizing the Kullback-Leibler divergence from the ideal proposal to the parametric family, and the cross-entropy of the two distributions is the quantity that the elite-likelihood maximization minimizes. The method does not need the event to be reachable in one step, only that a monotone score can be defined whose elite are closer to the event.

The failure mode is a family that cannot represent the event. If the true conditional distribution is multimodal and the family is a single Gaussian, the refit collapses to a compromise between the modes that is rare under every one of them, and the proposal never concentrates. The remedy is a richer family or a mixture, and the choice of family is the modeling decision that the method cannot make for itself. The same limitation appears in the rejection envelope, where the envelope constant is set by how well the proposal can represent the target.

The bias that a learned proposal introduces

Both methods are exact in principle and biased in practice, because the parameters they learn are estimated from finite samples. Splitting estimates the level-crossing ratios from a finite number of paths, and the per-level ratio estimates are biased slightly upward, which compounds across levels. The cross-entropy method fits the proposal to a finite elite, and the fitted proposal is not exactly the conditional distribution, so the importance-sampling estimator carries a bias that shrinks with the batch size but never vanishes.

The more serious bias is the one the practitioner introduces by looking. If the levels, the branching factor, or the elite fraction are tuned by inspecting the output and rerunning, the estimator inherits a selection effect: the reported variance is computed as if the design were fixed, but the design was chosen because it looked good, and the effective variance is larger. The standard defence is to fix the design on a pilot, then estimate on a fresh stream, exactly as common random numbers requires the paired runs to consume the stream in lockstep.

There is also a bias in the interpretation of the weights. A rare-event estimator that produces a handful of paths with large weights and thousands with negligible weights has an effective sample size far below the path count, and the effective sample size formula is the check that the weight distribution is not dominated by one lucky path.

What the methods can and cannot certify

Neither method can certify that an event was not missed. A rare-event estimator is only as good as the proposal or the levels, and if the event sits in a region the proposal never visits, the estimate is zero with a small reported variance. The only genuine check is a lower bound from an independent argument: a union bound, a comparison against a known special case, or a second method that reaches the event by a different route.

The cross-entropy method has one property that makes it attractive in practice: it produces the fitted proposal as a byproduct, and the proposal is a usable summary of where the tail is. A network model whose rare event is a buffer overflow produces, after cross-entropy fitting, a shifted traffic distribution that concentrates on the overload, and that distribution is the answer to the inverse question of what conditions cause the overflow, not just an intermediate in the estimate of its probability.

The methods compose rather than replace. Splitting is the tool when the rare event is a path property of a dynamic system, cross-entropy when a parametric shift can be fitted to a score, and plain importance sampling when the conditional distribution is available in closed form. All three are applications of the same principle, that the sampling distribution should be moved to the tail, and all three inherit the weight variance as their shared diagnostic of whether the move went too far.

References

  1. Rare Event Simulation using Monte Carlo Methods — Gerardo Rubino and Bruno Tuffin (eds.); Wiley
  2. A Tutorial on the Cross-Entropy Method — Pieter-Tjerk de Boer, Dirk P. Kroese, Shie Mannor and Reuven Y. Rubinstein; Annals of Operations Research, 2005
  3. Monte Carlo: Concepts, Algorithms, and Applications — George S. Fishman; Springer
  4. Simulation and the Monte Carlo Method, 3rd edition — Reuven Y. Rubinstein and Dirk P. Kroese; Wiley