Sequential Monte Carlo ABC for Mechanistic Models

Why Use a Sequential ABC Population?

In a mechanistic model, simulating data at a proposed parameter value may be feasible even when evaluating the likelihood is not. Approximate Bayesian computation (ABC) uses the simulation and a discrepancy between observed and simulated summaries to retain parameter values compatible with the observations. Rejection ABC can become inefficient when the tolerance is small: most prior draws may be simulated and discarded. Sequential Monte Carlo ABC (ABC-SMC) addresses this by moving a population of particles through a sequence of tolerances, using the population at one stage to propose candidates for the next.

The sequence does not make the final result exact. Its final population represents an ABC approximation whose accuracy depends on the final tolerance, summary statistics, discrepancy, and Monte Carlo quality. The purpose of the sequence is to reach a useful tolerance more efficiently while keeping track of how the distribution of particles changes.

Definition: Let $\theta$ denote a mechanistic parameter vector, with prior density $\pi(\theta)$. Let $s_{\mathrm{obs}}$ be the observed summary vector and $S_{\mathrm{sim}}$ the summary vector from a model simulation at $\theta$. Given a discrepancy $d(S_{\mathrm{sim}},s_{\mathrm{obs}})$ and tolerance $\epsilon$, the ABC likelihood factor is $$ L_\epsilon(\theta)=P\!\left(d(S_{\mathrm{sim}},s_{\mathrm{obs}})\leq\epsilon\mid\theta\right). $$ The tolerance-specific ABC target has density proportional to $\pi(\theta)L_\epsilon(\theta)$. It is an approximation to the likelihood-based posterior, not generally that exact posterior.

The discrepancy and summaries must be defined before running the sequence. In particular, summaries on different scales should be scaled consistently, as discussed in the earlier material on choosing and scaling ABC summary statistics. The tolerance is meaningful only relative to that discrepancy. A common schedule is non-increasing, $\epsilon_1\geq\epsilon_2\geq\cdots\geq\epsilon_T$, so later populations impose at least as strict an acceptance condition as earlier ones.

The ABC-SMC Sequence

At the first stage, draw parameters from the prior, simulate data, and retain draws whose discrepancy is at most $\epsilon_1$. These accepted draws approximate the first ABC target. At each later stage, select an ancestor from the preceding weighted population, perturb its parameter, and simulate again. Keep a candidate only if it lies in the prior support and its simulated discrepancy is at most the new tolerance. The accepted candidates are then assigned importance weights to account for the fact that they were proposed from the previous population rather than directly from the prior.

Initialize.
Choose $\epsilon_1$. Sample from the prior, simulate, and retain $N$ particles meeting the ABC condition. Give them equal weights $1/N$.

Set the next tolerance and kernel.
Choose $\epsilon_t\leq\epsilon_{t-1}$ and a perturbation kernel $K_t(\theta\mid\theta')$, often adapted from the preceding weighted particles.

Propose and accept.
Select an ancestor according to its normalized weight, perturb it, and simulate until a candidate passes the prior-support and discrepancy checks.

Reweight and assess.
Assign each accepted particle its importance weight, normalize the weights, and check acceptance, weight concentration, and particle diversity before proceeding.

This procedure assumes that the prior density and perturbation-kernel density can be evaluated, that the simulator can generate draws at proposed parameters, and that the recorded discrepancy is calculated consistently. Any constraints on parameters must be respected. If a transformed or truncated proposal is used, its actual proposal density—not an untruncated density that ignores the transformation or truncation—must be used in the weight calculation.

Why the Importance Weight Has This Form

At stage $t$, the proposal density for a candidate parameter is the mixture of perturbation kernels centered on the previous particles:

q_t(\theta)=\sum_{j=1}^{N}W_{t-1,j}K_t(\theta\mid\theta_{t-1,j}), \qquad \sum_{j=1}^{N}W_{t-1,j}=1.

Before the ABC acceptance check, candidates therefore come from $q_t$. Conditioning on passing the check changes their density in proportion to $q_t(\theta)L_{\epsilon_t}(\theta)$. The desired tolerance-specific target is proportional to $\pi(\theta)L_{\epsilon_t}(\theta)$. Taking the ratio cancels the common ABC likelihood factor, giving an importance ratio proportional to $\pi(\theta)/q_t(\theta)$. Thus an accepted particle $\theta_{t,i}$ receives the unnormalized weight

\widetilde{W}_{t,i} = \frac{\pi(\theta_{t,i})} {\displaystyle\sum_{j=1}^{N}W_{t-1,j} K_t(\theta_{t,i}\mid\theta_{t-1,j})}, \qquad W_{t,i}= \frac{\widetilde{W}_{t,i}} {\displaystyle\sum_{r=1}^{N}\widetilde{W}_{t,r}}.

This derivation requires that the proposal density cover the relevant target support: if a region with positive target density cannot be proposed, weighting cannot recover it. It also assumes the same ABC acceptance rule is used in the target and in candidate selection. The weights correct the proposal distribution; they do not correct a poor choice of summaries or a tolerance that is too large.

Worked Example: Reweighting a Candidate From a Kernel Mixture

Suppose the prior is uniform on $[0,2]$, so its density at an interior point is $0.5$. At one stage, three previous particles have weights $0.2$, $0.5$, and $0.3$. At the accepted candidate $\theta=0.4$, their perturbation-kernel densities are respectively $0.5794$, $0.3884$, and $0.0048$. Calculate the mixture density and the candidate's unnormalized importance weight.

The mixture proposal density at $0.4$ is

q_t(0.4) =0.2(0.5794)+0.5(0.3884)+0.3(0.0048) =0.11588+0.19420+0.00144 =0.31152 \approx0.3115.

The importance weight before normalization is the prior density divided by this mixture density:

\widetilde{W} =\frac{0.5}{0.31152} \approx1.605.

This is an unnormalized weight, not a probability. To obtain normalized weights, calculate the corresponding unnormalized weights for all accepted particles and divide each by their sum. The candidate's weight is relatively large because its prior density is appreciable compared with the density under which the mixture proposed it.

Adapting the Perturbation Kernel

A kernel that is too narrow tends to reproduce nearly identical particles and may fail to move between separated regions. A kernel that is too broad can propose many parameters with low ABC acceptance. One standard adaptive choice is a multivariate Gaussian random-walk kernel centered on each ancestor, with covariance based on the weighted covariance of the previous population. A common rule multiplies that covariance by two; this is a practical scale choice, not a guarantee of optimal acceptance or exploration.

K_t(\theta\mid\theta_{t-1,j}) = N\!\left(\theta_{t-1,j},\,2\widehat{\Sigma}_{t-1}\right).

Here $\widehat{\Sigma}_{t-1}$ is the weighted sample covariance. With normalized weights $W_j$, weighted mean $\bar{\theta}_w=\sum_j W_j\theta_j$, and $\sum_jW_j=1$, one common unbiased weighted estimate is

\widehat{\Sigma} = \frac{\displaystyle\sum_j W_j(\theta_j-\bar{\theta}_w)(\theta_j-\bar{\theta}_w)^\mathsf{T}} {\displaystyle 1-\sum_j W_j^2}.

The covariance should be estimated on a parameter scale where Gaussian perturbations are reasonable. For positive parameters, for example, proposing on a log scale can avoid invalid negative candidates, but the transformation and density must be handled consistently. If the covariance is singular or nearly singular, adding a small diagonal regularization term can make it usable; this does not fix a population that has already lost important regions.

Worked Example: Calculating a One-Parameter Adaptive Kernel

Suppose the previous stage has three equally weighted particles for a parameter $\theta$: $0.2$, $0.6$, and $1.0$. Calculate the weighted mean, the unbiased weighted variance, and the variance of a Gaussian perturbation using twice that variance.

The weighted mean is

\bar{\theta}_w=\frac{0.2+0.6+1.0}{3}=0.6.

The sum of weighted squared deviations is

\sum_jW_j(\theta_j-\bar{\theta}_w)^2 =\frac{1}{3}\left[(-0.4)^2+0^2+(0.4)^2\right] =\frac{0.32}{3} \approx0.1067.

Because $\sum_jW_j^2=3(1/3)^2=1/3$, the denominator in the weighted covariance estimate is $1-1/3=2/3$. Therefore,

\widehat{\Sigma} =\frac{0.1067}{2/3}=0.16, \qquad 2\widehat{\Sigma}=0.32.

The perturbation is thus centered on the selected ancestor with variance $0.32$, or standard deviation $\sqrt{0.32}\approx0.5657$, in the parameter's units. This calculation describes the proposal spread; the ABC filter and prior support still determine which candidates can be retained.

Choosing a Sequence and Checking Its Behavior

Tolerances may be set in advance or chosen adaptively, for example by a chosen quantile of discrepancies among accepted candidates. An adaptive rule should still enforce a non-increasing sequence and a defensible final tolerance. A quantile rule is a computational choice, not evidence by itself that the resulting approximation is scientifically adequate.

Worked Example: Applying a Decreasing Tolerance

At the first stage, suppose $\epsilon_1=2.0$. Five candidate simulations have scaled discrepancies $0.7$, $1.4$, $1.8$, $2.1$, and $2.6$. At the next stage, use $\epsilon_2=1.5$. Which first-stage candidates would pass the second-stage discrepancy condition if simulated again, and why is the first-stage accepted set not automatically the second-stage population?

The stage-two condition is $d\leq1.5$. Of the listed discrepancy values, $0.7$ and $1.4$ satisfy it; $1.8$, $2.1$, and $2.6$ do not. This identifies which values would pass if they were the discrepancies of stage-two simulations. In an actual ABC-SMC run, the new stage uses perturbed parameters and fresh simulations, so a particle's earlier discrepancy does not determine its new discrepancy. Each proposed candidate must be simulated and checked against $\epsilon_2$. Accepted candidates then receive the stage-two importance weights.

Diagnostics should be checked at every stage, not only at the end. The acceptance rate—the number accepted divided by the number proposed—can reveal a tolerance that is too stringent for the current proposal. The effective sample size of normalized weights, $1/\sum_i W_{t,i}^2$, indicates weight concentration: a value near one means a few particles dominate, while a value near $N$ indicates more even weights. Also inspect the spread and distinctness of particles, parameter boundaries, and whether plausible modes have disappeared. Low acceptance and low effective sample size together often signal that a stage is struggling; high acceptance alone does not establish that the approximation is accurate.

Alternatives, Limitations, and Common Mistakes

Approach	How it proceeds	Main trade-off
Rejection ABC	Draw from the prior and retain candidates passing one tolerance.	Simple, but can waste many simulations when the prior covers a large region and the tolerance is small.
ABC-SMC	Move a weighted population through a tolerance sequence.	Can reuse information from earlier stages, but requires careful proposal densities, weights, and population diagnostics.
ABC-based MCMC	Use a Markov chain whose proposals are accepted according to parameter and simulation-based criteria.	Can avoid storing a full sequential population, but depends on chain mixing and tuning rather than a sequence of weighted populations.

A central limitation is that ABC-SMC targets the posterior conditional on the chosen summaries and tolerance, not necessarily the posterior that would result from the full data likelihood. A decreasing tolerance can reduce the tolerance component of the approximation, but cannot restore information discarded by summaries. Nor can sequential reweighting reliably recover a mode that has been lost from the particle population. For stochastic simulators, variation in simulated summaries also affects acceptance, so repeated runs can be needed to assess Monte Carlo stability.

Mistake: treating accepted particles as equally weighted at every stage. At later stages, candidates come from a mixture of kernels, not directly from the prior. Use the prior density divided by the mixture proposal density, then normalize across accepted particles.
Mistake: calculating the mixture density using only the chosen ancestor. The candidate could have been generated from any previous particle. The denominator must sum kernel densities over the full preceding population, weighted by its normalized weights.
Mistake: assuming a smaller tolerance guarantees a better overall calibration. It tightens the discrepancy filter, but can cause very low acceptance and does not repair inadequate summaries, a poor discrepancy scale, or particle impoverishment.
Mistake: using a kernel density inconsistent with the actual proposal. Truncation, parameter transformations, and other modifications affect the proposal density and must be reflected in the importance weight.

Key takeaway: ABC-SMC combines a sequence of tolerance-specific ABC targets with importance-corrected proposals from the preceding weighted population. A reliable run needs a meaningful discrepancy, valid proposal densities, a non-increasing tolerance schedule, and diagnostics for acceptance, weight concentration, and particle diversity.

Check Your Understanding

Answer the questions, then compare your reasoning with the solutions below.

If the ABC tolerance decreases, does the ABC target necessarily become the exact likelihood-based posterior? Explain.
At a later stage, why is the importance-weight denominator a weighted sum over all preceding particles?
For normalized particle weights $0.5$, $0.3$, and $0.2$, calculate the effective sample size $1/\sum_iW_i^2$.
Give two signs that a perturbation kernel or tolerance schedule may be causing difficulties.
Why can a decreasing tolerance fail to recover a parameter mode that was missed by an earlier population?

Solutions. (1) No. It can reduce tolerance-related approximation, but the result still depends on the summaries and discrepancy; exact recovery is not guaranteed. (2) A candidate can arise from perturbation of any ancestor, so its proposal density is the full mixture density, with each kernel weighted by its ancestor's probability. (3) $1/(0.5^2+0.3^2+0.2^2)=1/0.38\approx2.632$. (4) Examples include very low acceptance, a low effective sample size, near-duplicate particles, loss of spread, or failure to retain plausible regions. (5) Later candidates are generated by perturbing particles already present. If no particle represents that mode and the kernels do not reach it, the sequence has no mechanism to recreate it.