Tutorial on ABC rejection and ABC SMC for parameter estimation and model selection

TL;DR

This tutorial explains four algorithms based on Approximate Bayesian Computation (ABC), including rejection and Sequential Monte Carlo (SMC) methods, for parameter estimation and model selection in complex models.

Contribution

It provides a clear schematic illustration of four ABC algorithms, enhancing understanding of their implementation for parameter estimation and model selection.

Findings

Clarifies the differences between rejection and SMC ABC methods.

Provides schematic guidance for implementing ABC algorithms.

Facilitates understanding of ABC for model selection.

Abstract

In this tutorial we schematically illustrate four algorithms: (1) ABC rejection for parameter estimation (2) ABC SMC for parameter estimation (3) ABC rejection for model selection on the joint space (4) ABC SMC for model selection on the joint space.

Figures

Full text

Tutorial on ABC rejection and ABC SMC for parameter estimation and model selection

Tina Toni, Michael P. H. Stumpf

ttoni@imperial.ac.uk, m.stumpf@imperial.ac.uk

In this tutorial we schematically illustrate four algorithms:

ABC rejection for parameter estimation [1, 2], 2. 2.

ABC SMC for parameter estimation [3, 4, 5], 3. 3.

ABC rejection for model selection on the joint space [6], 4. 4.

ABC SMC for model selection on the joint space [7].

We suggest to read this tutorial from the beginning. We start with a detailed explanation of the ABC rejection algorithm, which later helps to understand ABC SMC as it is based on the same concepts. Also, both model selection algorithms are closely related to parameter estimation algorithms and it is therefore helpful to understand those first.

This tutorial forms a part of the supplementary material of the paper ”Simulation-based model selection for dynamical systems in systems and population biology, Bioinformatics, 26 (1), 104-110, 2010” (T. Toni, M. P. H. Stumpf).

**ABC rejection

(a) We define a prior distribution $P(\theta)$ and we would like to approximate the posterior distribution $P(\theta|D_{0})$. We start by sampling a parameter $\theta^{*}$ from the prior distribution. We call this sampled parameter a particle.

(b) We simulate a data set $D^{*}$ according to some simulation framework $f(D|\theta^{*})$. In our examples we use different simulation frameworks. If we simulate a deterministic dynamical model, we add some noise at the time points of interest. If we simulate a stochastic dynamical model, we do not add any additional noise to the trajectories. We compare the simulated data set $D^{*}$ (circles) to the experimental data $D_{0}$ (crosses) using a distance function, $d$, and tolerance $\epsilon$; if $d(D_{0},D^{*})\leq\epsilon$, we accept $\theta^{*}$. The tolerance $\epsilon\geq 0$ is the desired level of agreement between $D_{0}$ and $D^{*}$.

(c) The particle $\theta^{*}$ is accepted because $D^{*}$ and $D_{0}$ are sufficiently close.

(d) We sample another parameter $\theta^{*}$ from the prior distribution and simulate a corresponding dataset $D^{*}$. In this case $D^{*}$ and $D_{0}$ are very different and we reject the particle (we ”throw it away”).

(e) We repeat the whole procedure until N particles have been accepted. They represent a sample from $P(\theta|d(D_{0},D^{*})\leq\epsilon)$, which approximates the posterior distribution. If $\epsilon$ is sufficiently small then the distribution $P(\theta|d(D_{0},D^{*})\leq\epsilon)$ will be a good approximation for the “true” posterior distribution, $P(\theta|D_{0})$.

(f) Many particles were rejected in the procedure, for which we have spent a lot of computational effort for simulation. ABC rejection is therefore computationally inefficient. We can use ABC SMC to reduce the computational cost.**

**ABC SMC (Toni et al., 2009)

(a) As in ABC rejection, we define a prior distribution $P(\theta)$ and we would like to approximate a posterior distribution $P(\theta|D_{0})$. In ABC SMC we do this sequentially by constructing intermediate distributions, which converge to the posterior distribution. We define a tolerance schedule $\epsilon_{1}>\epsilon_{2}>\ldots\epsilon_{T}\geq 0$.

(b) We sample particles from a prior distribution until $N$ particles have been accepted (have reached the distance smaller than $\epsilon_{1}$). For all accepted particles we calculate weights (see [4] for formulas and derivation). We call the sample of all accepted particles ”Population 1”.

(c) We then sample a particle $\theta^{*}$ from population 1 and perturb it to obtain a perturbed particle $\theta^{**}\sim K(\theta|\theta^{*})$, where K is a perturbation kernel (for example a Gaussian random walk). We then simulate a dataset $D^{*}\sim f(D|\theta^{**})$ and accept the particle $\theta^{**}$ if $d(D_{0},D^{**})\leq\epsilon_{2}$. We repeat this until we have accepted $N$ particles in population 2. We calculate weights for all accepted particles.

(d) We repeat the same procedure for the following populations, until we have accepted $N$ particles of the last population $T$ and calculated their weights. Population $T$ is a sample of particles that approximates the posterior distribution.

ABC SMC is computationally much more efficient than ABC rejection (see [4] for comparison).**

**ABC rejection for model selection

**

**ABC SMC for model selection

**