Tutorials › Stochastic Differential Equations in Pharmacokinetics › Multinomial, Stratified, Systematic, and Residual Resampling

Particle Filters and Bayesian Estimation · Tutorial 69 of 100

Multinomial, Stratified, Systematic, and Residual Resampling

Learn to turn normalized particle weights into ancestor indices using four standard resampling schemes, and understand what each scheme guarantees and where it can mislead.

Advanced 12 min read

What You'll Learn

  • Define offspring counts and ancestor indices from normalized particle weights.
  • Implement multinomial resampling using cumulative weights and independent uniforms.
  • Construct stratified and systematic resampling points and tally their offspring.
  • Apply residual resampling by allocating integer copies before sampling the remainder.
  • Derive why each scheme has expected offspring count equal to the particle’s weight times the sample size.
  • Compare variability, ordering sensitivity, and practical diagnostics across the four schemes.

From Weights to a New Particle Set

After a filtering update, a particle filter represents the current state distribution by particles \(x^{(1)},\ldots,x^{(N)}\) and normalized weights \(w_1,\ldots,w_N\), where \(w_i\geq 0\) and \(\sum_{i=1}^N w_i=1\). Resampling selects \(N\) ancestor indices from \(\{1,\ldots,N\}\), with selection probabilities determined by these weights. The new particle set consists of the selected particles, often with equal weights \(1/N\).

Resampling changes multiplicities, not particle locations: if particle \(i\) is selected three times, it contributes three copies to the resampled set. It is useful to record the result as either a list of ancestor indices or as offspring counts \(O_1,\ldots,O_N\), where \(O_i\) is the number of times particle \(i\) is selected and \(\sum_i O_i=N\).

Definition: Given normalized weights \(w_1,\ldots,w_N\), a resampling scheme produces integer counts \(O_i\geq 0\) with \(\sum_i O_i=N\). A scheme is unbiased in offspring counts if \(E[O_i]=Nw_i\) for every \(i\). This condition concerns expected counts; it does not say that every resampled set reproduces the weights exactly.

The four schemes below all satisfy this expected-count property, but they differ in how much randomness remains after the weights are given. To make their mechanics concrete, the examples use the same four weights and \(N=10\).

Four Ways to Allocate Offspring

Write cumulative weights as \(C_0=0\) and \(C_i=\sum_{k=1}^i w_k\). Associate particle \(i\) with the interval \((C_{i-1},C_i]\) on the unit interval. A uniform point in this interval selects particle \(i\). Whether an endpoint is included is immaterial for continuous uniforms, since any specified endpoint has probability zero.

Multinomial resampling: Draw \(N\) independent uniform random variables on \((0,1)\), and map each point to the particle interval containing it. The resulting ancestor indices are independent draws with probabilities \(w_1,\ldots,w_N\).

Worked Example: Multinomial Offspring Counts

Let the weights be \((0.05,0.15,0.40,0.40)\), with \(N=10\). Then the cumulative cut points are \(0.05,0.20,0.60,1.00\). Suppose the ten generated uniforms, in order, are \(0.07,0.12,0.18,0.25,0.32,0.48,0.55,0.66,0.81,0.97\).

The first three points lie in \((0.05,0.20]\), so they select particle 2. The next four lie in \((0.20,0.60]\), so they select particle 3. The final three lie in \((0.60,1]\), so they select particle 4; none lies in particle 1’s interval. Thus the ancestor list is \((2,2,2,3,3,3,3,4,4,4)\), and the offspring counts are \((0,3,4,3)\). They sum to 10, but they need not equal the expected counts \((0.5,1.5,4,4)\).

Stratified resampling: Divide \((0,1)\) into \(N\) equal strata, \(((j-1)/N,j/N)\). Draw one independent uniform point from each stratum, then map all \(N\) points to particles using the cumulative-weight intervals.

Worked Example: Stratified Offspring Counts

Use weights \((0.05,0.15,0.40,0.40)\) and \(N=10\). In stratum \(j\), choose the point \(u_j=(j-1+0.6)/10\). These points are \(0.06,0.16,0.26,0.36,0.46,0.56,0.66,0.76,0.86,0.96\).

The points \(0.06\) and \(0.16\) select particle 2; \(0.26,0.36,0.46,0.56\) select particle 3; and \(0.66,0.76,0.86,0.96\) select particle 4. Particle 1 receives none. The offspring counts are therefore \((0,2,4,4)\). The calculation illustrates the stratum-by-stratum method; in an implementation, each point is generated independently within its own stratum rather than fixed at the same relative position.

Systematic resampling: Draw one uniform offset \(U\) on \((0,1/N)\), form the ordered points \(U+(j-1)/N\) for \(j=1,\ldots,N\), and map them to particles. Unlike stratified resampling, the points are linked by the single shared offset.

Worked Example: Systematic Offspring Counts

For the same weights and \(N=10\), take the particular offset \(U=0.03\). The systematic points are \(0.03,0.13,0.23,0.33,0.43,0.53,0.63,0.73,0.83,0.93\).

The point \(0.03\) selects particle 1; \(0.13\) selects particle 2; \(0.23,0.33,0.43,0.53\) select particle 3; and \(0.63,0.73,0.83,0.93\) select particle 4. The resulting counts are \((1,1,4,4)\). This is one possible outcome, not a claim that the counts are fixed by the weights.

Residual resampling: First assign particle \(i\) the deterministic number \(\lfloor Nw_i\rfloor\) of offspring. Let \(R=N-\sum_i\lfloor Nw_i\rfloor\). If \(R>0\), draw the remaining \(R\) offspring multinomially with probabilities proportional to the residual masses \(Nw_i-\lfloor Nw_i\rfloor\). If \(R=0\), no random draws are needed.

Worked Example: Residual Offspring Counts

With weights \((0.05,0.15,0.40,0.40)\) and \(N=10\), the expected counts are \((0.5,1.5,4,4)\). Taking floors gives deterministic allocations \((0,1,4,4)\), which total 9, so \(R=1\) offspring remains.

The residual masses are \((0.5,0.5,0,0)\), whose sum is 1. Their normalized probabilities are therefore \((0.5,0.5,0,0)\). If the single residual draw uses a uniform value \(0.7\), the inverse cumulative probabilities select particle 2. Adding this draw to the deterministic allocations gives counts \((0,2,4,4)\). Particle 3 or 4 cannot receive the residual draw here because their expected counts were already integers.

Why the Expected Counts Agree

The schemes differ in dependence between selections, but each preserves the expected number of copies of every particle. This matters because for any function \(h\) of the particle state, the resampled average has the original weighted average as its expectation.

Theorem: Conditional on the particles and normalized weights, multinomial, stratified, systematic, and residual resampling each satisfy \(E[O_i]=Nw_i\) for every particle \(i\). Consequently, for any real-valued \(h\) with finite values on the particle set, \(E[\frac{1}{N}\sum_i O_i h(x^{(i)})]=\sum_i w_i h(x^{(i)})\).
1
Multinomial.
Write \(O_i=\sum_{j=1}^N I_{ji}\), where \(I_{ji}\) is 1 if draw \(j\) selects particle \(i\), and 0 otherwise. Since \(E[I_{ji}]=w_i\), linearity gives \(E[O_i]=\sum_j E[I_{ji}]=Nw_i\).
2
Stratified.
Let \(I_{ji}\) indicate that the point from stratum \(j\) falls in particle \(i\)’s interval. The probability is \(N\) times the length of the intersection of that stratum with \((C_{i-1},C_i]\), because the point is uniform on an interval of length \(1/N\). Summing over strata partitions the particle interval, whose length is \(w_i\). Thus \(E[O_i]=\sum_j E[I_{ji}]=Nw_i\).
3
Systematic.
For offset \(U\) uniform on \((0,1/N)\), the points are \(U+k/N\), \(k=0,\ldots,N-1\). For any interval \((a,b]\subseteq(0,1)\), the expected number of these points in it is \(N(b-a)\): integrating the number of lattice points in \((a,b]\) over one offset period gives the interval length times \(N\). Applying this to \((C_{i-1},C_i]\) yields \(E[O_i]=Nw_i\).
4
Residual.
Set \(a_i=\lfloor Nw_i\rfloor\), \(r_i=Nw_i-a_i\), and \(R=N-\sum_i a_i=\sum_i r_i\). If \(R>0\), the residual multinomial allocation has expected count \(R(r_i/R)=r_i\) for particle \(i\). Its total expected count is therefore \(a_i+r_i=Nw_i\). If \(R=0\), every \(r_i=0\), so the deterministic count \(a_i\) already equals \(Nw_i\).

For the consequence, write the resampled average as \(\frac{1}{N}\sum_i O_i h(x^{(i)})\). The particle values are fixed conditional on the filtering update, so linearity of expectation gives \(\frac{1}{N}\sum_i E[O_i]h(x^{(i)})=\frac{1}{N}\sum_i Nw_i h(x^{(i)})=\sum_i w_i h(x^{(i)})\). This proves the stated expectation identity. The result does not assert that resampling reduces error in every realized update. \(\square\)

Choosing a Scheme and Checking an Implementation

Multinomial resampling is simple and its draws are independent, but it can produce substantial count variation. For a fixed particle \(i\), \(O_i\) has a binomial distribution with parameters \(N,w_i\), so \(\operatorname{Var}(O_i)=Nw_i(1-w_i)\). Stratification reduces the randomness associated with placing points in separate parts of the cumulative-weight interval: each stratum contributes one point, rather than allowing several independent points to cluster in the same stratum. Systematic resampling uses only one random number and is often computationally convenient, but its points are dependent. Residual resampling removes the integer part of each expected count before random allocation, often reducing avoidable count variability.

These are not interchangeable guarantees. All four schemes preserve expected offspring counts, but a particular outcome can differ, and systematic resampling’s result depends on the particle ordering. If particles are arranged in a meaningful order, systematic points can interact with that order; if they are arbitrarily reordered, the realized allocation can change even though the weights and expected counts do not. Stratified and residual schemes also depend on the supplied weights and their numerical implementation, but do not use the same single-offset lattice.

In a particle filter, resampling can duplicate high-weight particles and eliminate low-weight particles. It does not create new state diversity. Frequent resampling may therefore cause genealogical concentration, while skipping resampling when weights are already even may add unnecessary Monte Carlo variation. Resampling decisions and later particle estimates are separate issues; the expectation identity above is a conditional property of the resampling step, not a guarantee about every downstream estimate.

Useful implementation checks include confirming that weights are nonnegative and sum to one within numerical tolerance, that the output has exactly \(N\) indices, that every index is in \(\{1,\ldots,N\}\), and that offspring counts sum to \(N\). For residual resampling, also verify that the residual masses are nonnegative and sum to \(R\). For systematic resampling, generate \(U\) on \((0,1/N)\), not on \((0,1)\); the latter would not place exactly one point in each of the intended equal-spacing positions. Cumulative sums should be monotone, and a search routine should handle a uniform point exactly at a cumulative boundary consistently. Such exact ties have probability zero in the mathematical algorithms but can arise through finite-precision arithmetic.

Common pitfall: It is incorrect to say that resampling “sets particle \(i\)’s count to \(Nw_i\).” Counts must be integers, so this is generally impossible. The correction is that each scheme targets \(Nw_i\) in expectation; a particular count is random except when residual allocation’s deterministic part, or the entire allocation, fixes it.

Practice and Solutions

Check Your Understanding

Use the definitions and calculations above to answer each question. Brief solutions follow the questions.

  1. For weights \((0.2,0.3,0.5)\) and \(N=5\), state the expected offspring count for each particle.
  2. In multinomial resampling with \(N=8\), what is the distribution and variance of the offspring count for a particle of weight \(0.25\)?
  3. For weights \((0.1,0.4,0.5)\) and \(N=6\), calculate the deterministic allocations and number of residual draws.
  4. Why does systematic resampling use one offset in \((0,1/N)\), rather than \(N\) independent uniforms on \((0,1)\)?

Solutions. (1) Multiply each weight by 5: the expected counts are \((1,1.5,2.5)\). (2) The count is binomial with parameters \(8\) and \(0.25\), so its variance is \(8(0.25)(0.75)=1.5\). (3) The expected counts are \((0.6,2.4,3)\), their floors are \((0,2,3)\), and the total deterministic allocation is 5; therefore one residual draw remains. Its residual masses are \((0.6,0.4,0)\), which sum to 1. (4) The one offset generates a lattice with spacing \(1/N\) and exactly \(N\) points, one corresponding to each systematic position. Independent uniforms instead define multinomial resampling and can cluster without that spacing.

Clinical Biostats