1. Why and how to resample

SMC relies on importance sampling. So let’s start quickly recapping that, and defining some notation. We approximate a target distribution \(\target(\cdot)\) with a set of samples from some surrogate proposal distribution, by assigning each sample an importance weight \(\impwt^\idx\) (proportional to the density ratio of target to proposal). Then the set of weighted particles \((\state^\idx,\impwt^\idx)_{\idx=1}^\np\) defines an empirical approximation to the current target \(\target(\cdot)\):

\[\widehat{\target}(\cdot) \triangleq \sum_{\idx=1}^\np \normwt^\idx \delta_{\state^\idx}(\cdot)\]

where \(\normwt^\idx = \impwt^\idx / \sum_j \impwt^j\) are the normalized weights.Throughout the rest of this post, “weights” and weights in code refer to the normalized weights \(\normwt^\idx\). For any test function \(f\), taking a weighted average of the function applied over the particles gives an estimator of its expectation under the target: \(\widehat{\E_{\target}[f]}\triangleq\E_{\widehat{\target}}[f]=\sum_\idx \normwt^\idx f(\state^\idx)\).

In sequential importance sampling (SIS), we apply importance sampling sequentially, to approximate a sequence of target distributions by maintaining an evolving population of particles, iteratively propagating and reweighting. One main issue with SIS is that it can suffer from weight degeneracy: When the weights become concentrated and the budget of \(\np\) particles effectively behaves as if it were just one sample, defeating the purpose of having multiple particles.

This weight degeneracy issue is one key motivation for SMC, which generalizes from SIS by the addition of resampling steps. Resampling replaces a current set of particles (potentially with degenerate weights) with a new set of samples with weights all set to be equal, duplicating high-weight particles and dropping low-weight ones. This addresses weight degeneracy, but introduces a different problem: path degeneracy.I’m getting this terminology for the two kinds of degeneracy from the excellent Naesseth et al. (2019, Chapter 2). They note that adaptive resampling can be a partial mitigation for path degeneracy. Only resample when the effective sample size \(\text{ESS} = 1/\sum_\idx (\normwt^\idx)^2\) (which ranges from 1 when one particle has all the weight to \(\np\) when weights are uniform) drops below a threshold (e.g., \(\np/2\)), rather than at every step. In this post I want to focus just on what happens when we do resample, rather than on when to resample, but the ESS values in the illustration can give a sense of when resampling would be triggered in an adaptive resampling method. After enough steps, all current particles may trace back to a single ancestor at earlier timesteps [click on a timestep label in the illustration to see ancestry].

Adding resampling can keep our particle approximation useful (fixing SIS’s weight degeneracy problem), but the method we use to resample can affect both the variance of our estimates and how quickly path diversity is lost. The rest of this post explores several resampling methods, visualizes their differences, to help build intuition for why some would perform better than others.

2. Inverse transform sampling

The first methods we’ll look at share the same core idea of inverse transform sampling. We can think of partitioning the unit interval into segments of width \(\normwt^\idx\), and mapping from ‘probe’ locations on the unit interval to determine the resampled particles. More precisely, the cumulative distribution function (CDF) of the discrete distribution defined by the weights, \(\cdf(\idx) = \sum_{j=1}^{\idx} \normwt^j\), defines this partition, and its inverse (the quantile function) it will map a value \(\probe\in(0,1]\) to a selected (resampled) particle \(\invcdf(\probe)\).

If you haven’t heard of inverse transform sampling before, this is actually a nice way to build the intuition. Think about how to randomly place \(\np\) probes so that any particle \(i\) gets selected \(\np\normwt^\idx\) times on average. Try clicking on the CDF plot at right to place probes, and consider what strategy would you use for placing \(\np\) probes, in order to have the distribution of resampled particles recreate the weight histogram on average.

Drag bars on left to change weights, or choose preset from dropdown. Click on plot to place ‘probes.’

3. Multinomial, stratified, and systematic resampling

You may have found that the most natural idea is to use \(\np\) independent draws from the uniform distribution.That is, inverse transform sampling. A little while ago, I made another post exploring density transformations in the continuous case. Here we are doing this in a discrete setting: Each uniform-distributed independent probe \(\probe_k\) is transformed through the discrete quantile function \(\invcdf\) to produce a sample from the particle-weight distribution, just as passing a uniform draw through a continuous inverse CDF yields a draw from the corresponding distribution. This works, and leads to the first standard algorithm, multinomial resampling. We will then look at two other common strategies that reduce variance by spreading probes more evenly for lower variance.

# Step 1: build the inverse CDF (shared by all three CDF-based methods)
cumulative_sum = np.cumsum(weights)
# Step 2: choose N independent uniform probes
positions = random(N)
# Step 3: map probes through the inverse CDF → ancestor indices
indices = np.searchsorted(cumulative_sum, positions)

# Steps 1 & 3 as above; only step 2 changes:
# one uniform draw per stratum instead of N independent
positions = (random(N) + range(N)) / N

# Steps 1 & 3 as above; only step 2 changes:
# a single random offset → N equally-spaced probes
positions = (random() + np.arange(N)) / N

4. Residual resampling

We could also consider a deterministic-stochastic hybrid, where some particles are deterministically set, and others are allocated randomly. In residual resampling, we give particle \(\idx\) exactly \(\lfloor \np\normwt^\idx \rfloor\) copies, then fill the remaining \(R = \np - \sum_\idx \lfloor \np\normwt^\idx \rfloor\) slots by resampling on the residual weights \(\resid^\idx = (\np\normwt^\idx - \lfloor \np\normwt^\idx \rfloor)/R\). This nondeterministic part of the algorithm could be done using any of the three CDF methods (select below).⊕\(\Var_{\text{resid}} \leq \Var_{\text{mult}}\) when phase 2 uses multinomial or stratified resampling (Douc et al., 2005). The deterministic phase removes variance for the integer parts; the phase-2 choice affects only the residual variance. Note that residual-systematic does not have this guarantee, since systematic resampling can hit the same counterexample on the residual weights. The right plot shows the residual CDF (solid) overlaid on the original weights CDF (dotted). For highly skewed weights, you can see that first allocating the deterministic weights makes the residual CDF much more even than the original was.

# Phase 1 (deterministic): guaranteed copies from the integer part
num_copies = np.floor(N * weights) # ⌊Nwⁱ⌋
R = N - sum(num_copies)  # remaining slots
# Phase 2 (stochastic): multinomial on the fractional remainders
residual = weights * N - num_copies
residual /= sum(residual)# normalize residuals
positions = random(R)    # multinomial: R independent probes
indexes[k:N] = np.searchsorted(cumsum(residual), positions)

Phase 2 method: Multinomial Stratified Systematic

5. Comparison

Comparing all four methods on the same weights: Colored error bars show mean count \(\pm\) 1 std over \(K\) trials.

6. Other resampling schemes

The four methods above are the most widely used, but they are not the only options.

# For each particle independently:
num_copies = np.floor(N * weights)
p_bonus = N * weights - num_copies # fractional part
u = np.random.rand(N)    # one uniform draw per particle
bonus = (u >= 1 - p_bonus)    # inverse CDF: right of step → bonus
num_copies += bonus # total may differ from N

7. Resampling methods in action

To see how the choice of resampling method affects particle diversity in practice, here is the same random walk model from the introduction, now with a method selector.