In many mathematical biology modelling, we eventually encounter the problem of coupling experimental data with simulation data. In many of these cases, we also care about uncertainty quantification of our model parameters in prediciton experimental results. Once approach for parameter inference is Bayesian inference. However, in many modelling senarios the likelihood is intractable meaning the likelihood is unavailable in closed form, or where evaluation of the likelihood is infeasible. Hence, we turn to approximating of the posterior distribution using a technique called approximate Bayesian computation (ABC). This blog is primiarly based off two papers (Sisson, Fan, and Tanaka 2007) and (Marjoram et al. 2003)

A brief introduction

This blog is dedicated to my findings and implementations of simple approximate Bayesian computation (ABC) algorithms. The motivation is to infer parameters of an agent-based model from experimental results in my PhD project. However, in this blog I will only note down the basic concepts and a implementation of a Metropolis-Hastings (MH) Monte Carlo Markov Chain (MCMC) method will be shown in MATLAB.

Algorithms

In ABC inferencing the first appoach is rejection sampling where we accept or reject samples drawn from a distibution. A template ABC rejection sampling algorithm is given below.

Rejection sampling

  1. Generate a candidate value \(\theta^* \sim \pi(\theta)\) from the prior.
  2. Generate a data set \(x^*\sim f(x|\theta^*).\)
  3. Accept \(\theta^*\) if \(\rho(S(x^*),S(x_0))\leq \epsilon\).
  4. If rejected, go to step 1.

However, the samples are uncorrelated meaning the acceptance rate is relatively low and therefore more computational expensive. This motivates us to find a better sampling algorithm that has a higher acceptance rate. The next best thing is MCMC as given below.

Markov chain Monte Carlo (MCMC)

  1. Initialise \(\theta_1 ,i=1\).
  2. Generate a candidate value \(\theta^*\sim q(\theta|\theta_i)\), where \(q\) is some proposal density.
  3. Generate a data set \(x^* \sim f(x|\theta^*).\)
  4. Set \(\theta_{i+1}=\theta^*\) with probability

\[ \alpha = \min \left\{ 1, \frac{\pi(\theta^*)q(\theta_i|\theta^*)}{\pi(\theta_i)q(\theta^*|\theta_i)} \boldsymbol{1}(\rho(S(x^*),S(x_0)) \leq \epsilon) \right\}, \] otherwise set \(\theta_{i+1}=\theta_i\).

  1. If \(i<N\), increment \(i=i+1\) and go to 2.

An example

An example for simulating yeast colonies from an agent-based model will be shown. TBA. See diagram below for an outline of the procedure.

References

Marjoram, Paul, John Molitor, Vincent Plagnol, and Simon Tavaré. 2003. “Markov Chain Monte Carlo Without Likelihoods.” Proceedings of the National Academy of Sciences 100 (26): 15324–28.
Sisson, Scott A, Yanan Fan, and Mark M Tanaka. 2007. “Sequential Monte Carlo Without Likelihoods.” Proceedings of the National Academy of Sciences 104 (6): 1760–65.