Multilevel rejection sampling for approximate Bayesian computation

The most elementary implementation of ABC is ABC rejection sampling (Pritchard et al., 1999; Sunnåker et al., 2013), see Algorithm 1. This method generates $N$ independent and identically distributed samples $\bm{\theta}^{1},\ldots,\bm{\theta}^{N}$ from the posterior distribution by accepting proposals, $\bm{\theta}^{*}~\sim~p(\bm{\theta})$, when the data generated by the model simulation process $f(\mathcal{D}\mid\nobreak\bm{\theta}^{*})$ is within $\epsilon$ of the observed data, $\mathcal{D}$, under the discrepancy metric $d(\mathcal{D},\cdot)$. ABC rejection sampling is simple to implement, and samples are independent and identically distributed. Therefore ABC rejection sampling is widely used in many applications (Browning et al., 2018; Grelaud et al., 2009; Navascués et al., 2017; Ross et al., 2017; Vo et al., 2015). However, ABC rejection sampling can be computationally prohibitive in practice (Barber et al., 2015; Fearnhead and Prangle, 2012). This is especially true when the prior density is highly diffuse compared with the target posterior density (Marin et al., 2012), as most proposals are rejected.

To improve the efficiency of ABC rejection sampling, one can consider a likelihood-free modification of Markov chain Monte Carlo (MCMC) (Beaumont et al., 2002; Marjoram et al., 2003; Tanaka et al., 2006) in which a Markov chain is constructed with a stationary distribution identical to the desired posterior. Given the Markov chain is in state $\bm{\theta}^{i}$, a state transition is proposed via a proposal kernel, $K(\bm{\theta}\mid\nobreak\bm{\theta}^{i})$.

The Metropolis-Hastings (Hastings, 1970; Metropolis et al., 1953) state transition probability, $h$, can be modified within an ABC framework to yield

The stationary distribution of such a Markov chain is the desired approximate posterior (Marjoram et al., 2003). Algorithm 2 provides a method for computing $N_{T}$ iterations of this Markov chain.

While MCMC-ABC sampling can be highly efficient, the samples in the sequence, $\bm{\theta}^{1},\ldots,\bm{\theta}^{N_{T}}$, are not independent. This can be problematic as it is possible for the Markov chain to take long excursions into regions of low posterior probability. This incurs additional bias that is potentially significant (Sisson et al., 2007). A poor choice of proposal kernel can also have considerable impact upon the efficiency of MCMC-ABC (Green et al., 2015). The question of how to choose the proposal kernel is non-trivial. Typically proposal kernels are determined heuristically. However, automatic and adaptive schemes are available to assist in obtaining near optimal proposals in some cases (Cabras et al., 2015; Roberts and Rosenthal, 2009). Another additional complication is that of determining when the Markov Chain has converged; this is a challenging problem to solve in practice (Roberts and Rosenthal, 2004).

Sequential Monte Carlo (SMC) sampling was introduced to address these potential inefficiencies (Del Moral et al., 2006) and later extended within an ABC context (Sisson et al., 2007; Drovandi and Pettitt, 2011; Toni et al., 2009). A set of samples, referred to as particles, is evolved through a sequence of ABC posteriors defined through a sequence of $T$ acceptance thresholds, $\epsilon_{1},\ldots,\epsilon_{T}$ (Sisson et al., 2007; Beaumont et al., 2009). At each step, $t\in[0,T]$, the current ABC posterior, $p(\bm{\theta}\mid\nobreak d(\mathcal{D},\mathcal{D}_{s})<\epsilon_{t})$, is approximated by a discrete distribution constructed from a set of $N_{P}$ particles $\bm{\theta}_{t}^{1},\ldots,\bm{\theta}_{t}^{N_{P}}$ with importance weights $W_{t}^{1},\ldots,W_{t}^{N_{P}}$ that is, $\mathbb{P}(\bm{\theta}=\bm{\theta}_{t}^{i})=W_{t}^{i}$ for $i=1,\ldots,N_{P}$. The collection is updated to represent the ABC posterior of the next step, $p(\bm{\theta}\mid\nobreak d(\mathcal{D},\mathcal{D}_{s})<\epsilon_{t+1})$, through application of rejection sampling on particles perturbed by a proposal kernel, $K(\bm{\theta}\mid\nobreak\bm{\theta}^{i-1})$. If $p(\bm{\theta}\mid\nobreak d(\mathcal{D},\mathcal{D}_{s})<\epsilon_{t})$ is similar to $p(\bm{\theta}\mid\nobreak d(\mathcal{D},\mathcal{D}_{s})<\epsilon_{t+1})$, the acceptance rate should be high. The importance weights of the new family of particles are updated using an optimal backwards kernel (Sisson et al., 2007; Beaumont et al., 2009). Algorithm 3 outlines the process.

Through the use of independent weighted particles, SMC-ABC avoids long excursions into the distribution tails that are possible with MCMC-ABC. However, SMC-based methods can be affected by particle degeneracy, and the efficiency of each step is still dependent on the choice of proposal kernel (Green et al., 2015; Filippi et al., 2013). Just as with MCMC-ABC, adaptive schemes are available to guide proposal kernel selection (Beaumont et al., 2009; Del Moral et al., 2012). The choice of the sequence of acceptance thresholds is also important for the efficiency of SMC-ABC. However, there are good solutions to generate these sequences in an adaptive manner (Drovandi and Pettitt, 2011; Silk et al., 2013).

Multilevel Monte Carlo (MLMC) is a recent development that can significantly reduce the computational burden in the estimation of expectations (Giles, 2008). To demonstrate the basic idea of MLMC, consider computing the expectation of a continuous-time stochastic process $X_{t}$ at time $T$. Let $Z_{t}^{\tau}$ denote a discrete-time approximation to $X_{t}$ with time step $\tau$: the expectations of $X_{T}$ and $Z_{T}^{\tau}$ are related according to

That is, an estimate of $\mathbb{E}\left[Z_{T}^{\tau}\right]$ can be treated as a biased estimate of $\mathbb{E}\left[X_{T}\right]$. By taking a sequence of time steps $\tau_{1}>\cdots>\tau_{L}$, the indices of which are referred to as levels, we can arrive at a telescoping sum,

Computing this form of the expectation returns the same bias as that returned when computing $\mathbb{E}\left[Z_{T}^{\tau_{L}}\right]$ directly. However, Giles (2008) demonstrates that a Monte Carlo estimator for the telescoping sum can be computed more efficiently than directly estimating $\mathbb{E}\left[Z_{T}^{\tau_{L}}\right]$ in the context of stochastic differential equations (SDEs). This efficiency comes from exploiting the fact that the bias correction terms, $\mathbb{E}\left[Z_{T}^{\tau_{\ell}}-Z_{T}^{\tau_{\ell-1}}\right]$, measure the expected difference between the estimates on levels $\ell$ and $\ell-1$. Therefore, sample paths of $Z_{T}^{\tau_{\ell-1}}$ need not be independent of sample paths of $Z_{T}^{\tau_{\ell}}$. In the case of SDEs, samples may be generated in pairs driven by the same underlying Brownian motion, that is, the pair is coupled. By the strong convergence properties of numerical schemes for SDEs, Giles (2008) shows that this coupling is sufficient to reduce the variance of the Monte Carlo estimator. This reduction in variance is achieved through optimally trading off statistical error and computational cost across all levels. Through this trade-off, an estimator is obtained with the same accuracy in mean-square to standard Monte Carlo, but at a reduced computational cost. This saving of computational cost is achieved since fewer samples of the most accurate discretisation are required.

Recently, several examples of MLMC applications to Bayesian inference problems have appeared in the literature. One of the biggest challenges in the application of MLMC to inverse problems is the introduction of a strong coupling between levels. That is, the construction of a coupling mechanism that reduces the variances of the bias correction terms enough to enable the MLMC estimator to be computed more efficiently than standard Monte Carlo. Dodwell et al. (2015) demonstrate a MLMC scheme for MCMC sampling applicable to high-dimensional Bayesian inverse problems with closed-form likelihood expressions. The coupling of Dodwell et al. (2015) is based on correlating Markov chains defined on a hierarchy in parameter space. A similar approach is also employed by Efendiev et al. (2015). A multilevel method for ensemble particle filtering is proposed by Gregory et al. (2016) that employs an optimal transport problem to correlate a sequence of particle filters of different levels of accuracy. Due to the computational cost of the transport problem, a local approximation scheme is introduced for multivariate parameters (Gregory et al., 2016). Beskos et al. (2017) look more generally at the case of applying MLMC methods when independent and identically distributed sampling of the distributions on some levels is infeasible, the result is an extension of MLMC in which coupling is replaced with sequential importance sampling, that is, a multilevel variant of SMC (MLSMC).

MLMC has also recently been considered in an ABC context. Guha and Tan (2017) extend the work of Efendiev et al. (2015) by replacing the Metropolis-Hasting acceptance probability in a similar way to the MCMC-ABC method (Marjoram et al., 2003). The MLSMC method (Beskos et al., 2017) is exploited to achieve coupling in an ABC context by Jasra et al. (2017).

ABC samplers based on MCMC and SMC are generally more computationally efficient than ABC rejection sampling (Marjoram et al., 2003; Sisson et al., 2007; Toni et al., 2009). However, there are many advantages to using ABC rejection sampling. Specific advantages are: (i) its simple implementation; (ii) it produces truly independent and identically distributed samples, that is, there is no need to re-weight samples; and (iii) there are no algorithm parameters that affect the computational efficiency. In particular, no proposal kernels need be heuristically defined for ABC rejection sampling.

The aim of this work is to design a new algorithm for ABC inference that retains the aforementioned advantages of ABC rejection sampling, while still being efficient and accurate. The aim is not to develop a method that is always superior to MCMC and SMC methods, but rather provide a method that requires little user-defined configuration, and is still computationally reasonable. To this end, we investigate MLMC methods applied directly to ABC rejection sampling.

Various applications of MLMC to ABC inference have been demonstrated very recently in the literature (Guha and Tan, 2017; Jasra et al., 2017), with each implementation contributing different ideas for constructing effective control variates in an ABC context. It is not yet clear when one approach will be superior to another. Thus the question of how best to apply MLMC in an ABC context remains an open problem. Since there are no examples of an application of MLMC methods to the most fundamental of ABC samplers, that is rejection sampling, such an application is a significant contribution to this field.

We describe a new algorithm for ABC rejection sampling, based on MLMC. Our new algorithm, called MLMC-ABC, is as general as standard ABC rejection sampling but is more computationally efficient through the use of variance reduction techniques that employ a novel construction for the coupling problem. We describe the algorithm, its implementation, and validate the method using a tractable Bayesian inverse problem. We also compare the performance of the new method against MCMC-ABC and SMC-ABC using a standard benchmark problem from epidemiology.

Our method benefits from the simplicity of ABC rejection sampling. We require only the discrepancy metric and a sequence of acceptance thresholds to be defined for a given inverse problem. Our approach is also efficient, and achieves comparable or superior performance to MCMC-ABC and SMC-ABC methods, at least for the examples considered in this paper. Therefore, we demonstrate that our algorithm is a promising method that could be extended to design viable alternatives to current state-of-the-art approaches. This work, along with that of Guha and Tan (2017) and Jasra et al. (2017), provides an additional set of computational tools to further enhance the utility of ABC methods in practice.

We first reformulate the Bayesian inference problem (Equation (1)) as an expectation. To this end, note that, given a $k$-dimensional parameter space, $\bm{\Theta}$, and the random parameter vector, $\bm{\theta}$, if the joint posterior CDF, $F(\mathbf{s})~=~\mathbb{P}(\theta_{1}\leq s_{1},\ldots,\theta_{k}\leq s_{k})$ is differentiable, that is, its probability density function (PDF) exists, then

where $A_{\mathbf{s}}=\{\bm{\theta}\in\bm{\Theta}:\theta_{1}\leq s_{1},\ldots,\theta_{k}\leq s_{k}\}$. This can be expressed as an expectation by noting

where $\mathds{1}_{A_{\mathbf{s}}}(\bm{\theta})$ is the indicator function: $\mathds{1}_{A_{\mathbf{s}}}(\bm{\theta})=1$ whenever $\bm{\theta}\in A_{\mathbf{s}}$ and $\mathds{1}_{A_{\mathbf{s}}}(\bm{\theta})=0$ otherwise. Now, consider the ABC approximation given in Equation (2) with discrepancy metric $d(\mathcal{D},\mathcal{D}_{s})$ and acceptance threshold $\epsilon$. The ABC posterior CDF, denoted by $F_{\epsilon}(\mathbf{s})$, will be

The marginal ABC posterior CDFs, $F_{\epsilon,1}(s),\ldots,F_{\epsilon,k}(s)$, are

We now introduce some notation that will simplify further derivations. We define $\bm{\theta}_{\epsilon}$ to be a random vector distributed according to the ABC posterior CDF, $F_{\epsilon}(\mathbf{s})$, with acceptance threshold, $\epsilon$, as given in Equation (2.1). This provides us with simplification in notation for the ABC posterior PDF, $p(\bm{\theta}_{\epsilon})=p(\bm{\theta}\mid\nobreak d(\mathcal{D},\mathcal{D}_{s})<\epsilon)$, and the conditional expectation, $\mathbb{E}\left[\mathds{1}_{A_{\mathbf{s}}}(\bm{\theta}_{\epsilon})\right]=\mathbb{E}\left[\mathds{1}_{A_{\mathbf{s}}}(\bm{\theta})\mid\nobreak d(\mathcal{D},\mathcal{D}_{s})<\epsilon\right]$. For any expectation, $P$, we use $\hat{P}^{N}$ to denote the Monte Carlo estimate of this expectation using $N$ samples.

The standard Monte Carlo integration approach is to generate $N$ samples $\bm{\theta}_{\epsilon}^{1},\ldots,\bm{\theta}_{\epsilon}^{N}$ from the ABC posterior, $p(\bm{\theta}_{\epsilon})$, then evaluate the empirical CDF (eCDF),

for $\mathbf{s}\in\mathcal{S}$, where $\mathcal{S}$ is a discretisation of the parameter space $\bm{\Theta}$. For simplicity, we will consider $\mathcal{S}$ to be a $k$-dimensional regular lattice. In general, a regular lattice will not be well suited for high dimensional problems. However, since this is the first time that this MLMC approach has been presented, it is most natural to begin the exposition with a regular lattice, and then discuss other more computationally efficient approaches later. More attention is given to other possibilities in Section 4.

The eCDF is not, however, the only Monte Carlo approximation to the CDF one may consider. In particular, Giles et al. (2015) demonstrate the application of MLMC to a univariate CDF approximation. We now present a multivariate equivalent of the MLMC CDF of Giles et al. (2015) in the context of ABC posterior CDF estimation. Consider a sequence of $L$ acceptance thresholds, $\{\epsilon_{\ell}\}_{\ell=1}^{\ell=L}$, that is strictly decreasing, that is, $\epsilon_{\ell}>\epsilon_{\ell+1}$. In this work, the problem of constructing optimal sequences is not considered, as the focus is the initial development of the method. More discussion around this problem is given in Section 4. Given such a sequence, $\{\epsilon_{\ell}\}_{\ell=1}^{\ell=L}$, we can represent the CDF (Equation (2.1)) using the telescoping sum

where

Using our notation, the MLMC estimator for Equation (6) and Equation (7) is given by

where

and $g_{\mathbf{s}}(\bm{\theta})$ is a Lipschitz continuous approximation to the indicator function; this approximation is computed using a tensor product of cubic polynomials,

where $\delta_{j}$ is the lattice spacing in the $j$th dimension, and $\xi(x)$ is a piece-wise continuous polynomial,

This expression is based on the rigorous treatment of smoothing required in the univariate case given by Giles et al. (2015). While other polynomials that satisfy certain conditions are possible (Giles et al., 2015; Reiss, 1981), here we restrict ourselves to this relatively simple form. Application of the smoothing function improves the quality of the final CDF estimate, just as using smoothing kernels improves the quality of PDF estimators (Silverman, 1986). Such a smoothing is also necessary to avoid convergence issues with MLMC caused by the discontinuity of the indicator function (Avikainen, 2009; Giles et al., 2015).

To compute the $\hat{Y}^{N_{1}}_{1}(\mathbf{s})$ term (Equation (9)), we generate $N_{1}$ samples $\bm{\theta}_{\epsilon_{1}}^{1},\ldots,\bm{\theta}_{\epsilon_{1}}^{N_{1}}$ from $p(\bm{\theta}_{\epsilon_{1}})$; this represents a biased estimate for $F_{\epsilon_{L}}(\mathbf{s})$. To compensate for this bias, correction terms, $\hat{Y}^{N_{\ell}}_{\ell}(\mathbf{s})$, are evaluated for $\ell>2$, each requiring the generation of $N_{\ell}$ samples $\bm{\theta}_{\epsilon_{\ell}}^{1},\ldots,\bm{\theta}_{\epsilon_{\ell}}^{N_{\ell}}$ from $p(\bm{\theta}_{\epsilon_{\ell}})$ and $N_{\ell}$ samples $\bm{\theta}_{\epsilon_{\ell-1}}^{1},\ldots,\bm{\theta}_{\epsilon_{\ell-1}}^{N_{\ell}}$ from $p(\bm{\theta}_{\epsilon_{\ell-1}})$, as given in Equation (9). It is important to note that the samples, $\bm{\theta}_{\epsilon_{\ell-1}}^{1},\ldots,\bm{\theta}_{\epsilon_{\ell-1}}^{N_{\ell}}$, used to compute $\hat{Y}^{N_{\ell}}_{\ell}(\mathbf{s})$ are independent of the samples, $\bm{\theta}_{\epsilon_{\ell-1}}^{1},\ldots,\bm{\theta}_{\epsilon_{\ell-1}}^{N_{\ell-1}}$, used to compute $\hat{Y}^{N_{\ell-1}}_{\ell-1}(\mathbf{s})$.

The goal is to introduce a coupling between levels that controls the variance of the bias correction terms. With an effective coupling, the result is an estimator with lower variance, hence the number of samples required to obtain an accurate estimate is reduced. Denote $v_{\ell}$ as the variance of the estimator $\hat{Y}^{N_{\ell}}_{\ell}(\mathbf{s})$. For $\ell\geq 2$ this can be expressed as

where $\mathbb{V}\left[\cdot\right]$ and $\mathbb{C}\left[\cdot,\cdot\right]$ denote the variance and covariance, respectively. Introducing a positive correlation between the random variables $\bm{\theta}_{\epsilon_{\ell}}$ and $\bm{\theta}_{\epsilon_{\ell-1}}$ will have the desired effect of reducing the variance of $\hat{Y}^{N_{\ell}}_{\ell}(\mathbf{s})$.

In many applications of MLMC, a positive correlation is introduced through driving samplers at both the $\ell$ and $\ell-1$ level with the same randomness. Properties of Brownian motion or Poisson processes are typically used for the estimation of expectations involving SDEs or Markov processes (Giles, 2008; Anderson and Higham, 2012; Lester et al., 2016). In the context of ABC methods, however, simulation of the quantity of interest is necessarily based on rejection sampling. The reliance on rejection sampling makes a strong coupling, in the true sense of MLMC, a difficult, if not impossible task. Rather, here we introduce a weaker form of coupling through exploiting the fact that our MLMC estimator is performing the task of computing an estimate of the ABC posterior CDF. We combine this with a property of nested ABC rejection samplers to arrive at an efficient algorithm for computing $\hat{F}_{\epsilon_{L}}(\mathbf{s})$.

We proceed to establish a correlation between levels as follows. Assume we have computed, for some $\ell<L$, the terms $\hat{Y}^{N_{1}}_{1}(\mathbf{s}),\ldots,\hat{Y}^{N_{\ell}}_{\ell}(\mathbf{s})$ in Equation (8). That is, we have an estimator to the CDF at level $\ell$ by taking the sum

with marginal distributions $\hat{F}^{N_{1},\ldots,N_{\ell}}_{\epsilon_{\ell},j}(s_{j})$ for $j=1,\ldots,k$. We can use this to determine a coupling based on matching marginal probabilities when computing $\hat{Y}_{\ell+1}(\mathbf{s})$. After generating $N_{\ell+1}$ samples $\bm{\theta}_{\epsilon_{\ell+1}}^{1},\ldots,\bm{\theta}_{\epsilon_{\ell+1}}^{N_{\ell+1}}$ from $p(\bm{\theta}_{\epsilon_{\ell+1}})$, we compute the eCDF, $\hat{F}^{N_{\ell+1}}_{\epsilon_{\ell+1}}(\mathbf{s})$, using Equation (5) and obtain the marginal distributions, $\hat{F}^{N_{\ell+1}}_{\epsilon_{\ell+1},j}(s_{j})$ for $j=1,\ldots,k$. We can thus generate $N_{\ell+1}$ coupled pairs $\{\bm{\theta}^{i}_{\epsilon_{\ell+1}},\bm{\theta}^{i}_{\epsilon_{\ell}}\}$ by choosing the $\bm{\theta}^{i}_{\epsilon_{\ell}}$ with the same marginal probabilities as the empirical probability of $\bm{\theta}^{i}_{\epsilon_{\ell+1}}$. That is, the $j$th component of $\bm{\theta}^{i}_{\epsilon_{\ell}}$ is given by

where $\theta^{i}_{\epsilon_{\ell+1}}$ is the $j$th component of $\bm{\theta}^{i}_{\epsilon_{\ell+1}}$ and $\hat{G}^{N_{1},\ldots,N_{\ell}}_{\epsilon_{\ell},j}(s)$ is the inverse of the $j$th marginal distribution of $\hat{F}^{N_{1},\ldots,N_{\ell}}_{\epsilon_{\ell}}(\mathbf{s})$. This introduces a positive correlation between the sample pairs, $\{\bm{\theta}_{\epsilon_{\ell+1}}^{i},\bm{\theta}_{\epsilon_{\ell}}^{i}\}$, since an increase in any of the components of $\bm{\theta}_{\epsilon_{\ell+1}}^{i}$ will cause an increase in the same component $\bm{\theta}_{\epsilon_{\ell}}^{i}$. This correlation reduces the variance in the bias correction estimator $\hat{Y}^{N_{\ell+1}}_{\epsilon_{\ell+1}}(\mathbf{s})$ computed according to Equation (9). We can then update the MLMC CDF to get an improved estimator by using

and apply an adjustment that ensures monotonicity. We continue this process iteratively to obtain $\hat{F}^{N_{1},\dots,N_{L}}_{\epsilon_{L}}(\mathbf{s})$.

It must be noted here that this coupling mechanism introduces an approximation for the general inference problem; therefore some additional bias can be introduced. This is made clear when one considers the process in terms of the copula distributions of $\bm{\theta}_{\epsilon_{\ell+1}}$ and $\bm{\theta}_{\epsilon_{\ell}}$. If these copula distributions are the same, then the coupling is exact and there is no additional bias. The coupling is also exact for the univariate case ($k=1$). Therefore, under the assumption that the correlation structure does not change significantly between levels, then the bias should be low. In practice, this requirement affects the choice of the acceptance threshold sequence, $\epsilon_{1},\ldots,\epsilon_{L}$; we discuss this in more detail in Section 4. In Section 3.1, we demonstrate that for sensible choices of this sequence, the introduced bias is small compared with the bias that is inherent in the ABC approximation.

We now require the sample numbers $N_{1},\ldots,N_{L}$ that are the optimal trade-off between accuracy and efficiency. Denote $d_{\ell}$ as the number of data generation steps required during the computation of $\hat{Y}^{N_{\ell}}_{\ell}(\mathbf{s})$ and let $c_{\ell}=d_{\ell}/N_{\ell}$ be the average number of data generation steps per accepted ABC posterior sample using acceptance threshold $\epsilon_{\ell}$.

Given $v_{1}~=~\mathbb{V}\left[g_{\mathbf{s}}(\bm{\theta}_{\epsilon_{1}})\right]$ and $v_{\ell}~=~\mathbb{V}\left[g_{\mathbf{s}}(\bm{\theta}_{\epsilon_{\ell}})-g_{\mathbf{s}}(\bm{\theta}_{\epsilon_{\ell-1})}\right]$, for $\ell~>~1$, one can construct the optimal $N_{\ell}$ under the constraint $\mathbb{V}\left[\hat{F}^{N_{1},\ldots,N_{L}}_{\epsilon_{L}}(\mathbf{s})\right]=\mathcal{O}(h^{2})$, where $h^{2}$ is the target variance of the MLMC CDF estimator. As shown by Giles (2008), using a Lagrange multiplier method, the optimal $N_{1},\ldots,N_{L}$ are given by

In practice, the values for $v_{1},\ldots,v_{L}$ and $c_{1},\ldots,c_{L}$ will not have analytic expressions available; rather, we perform a low accuracy trial simulation with all $N_{1}=\cdots=N_{L}=c$, for some comparatively small constant, $c$, to obtain the relative scaling of variances and data generation requirements.

A MLMC method based on the estimator in Equation (8) and the variance reduction strategy given in Section 2.3 would depend on standard ABC rejection sampling (Algorithm 1) for the final bias correction term $\hat{Y}^{N_{L}}_{\epsilon_{L}}$. For many ABC applications of interest, the computation of this final term will dominate the computational costs. Therefore, the potential computational gains depend entirely on the size of $N_{L}$ compared to the number of samples, $N$, required for the equivalent standard Monte Carlo approach (Equation (5)). While this approach is often an improvement over rejection sampling, we can achieve further computational gains by exploiting the iterative computation of the bias correction terms.

Let $\operatorname*{supp}(f(x))$ denote the support of a function $f(x)$, and note that, for any $\ell\in[2,L]$, $\operatorname*{supp}(p(\bm{\theta}_{\epsilon_{\ell}}))\subseteq\operatorname*{supp}(p(\bm{\theta}_{\epsilon_{\ell-1}}))$. This follows from the fact that if, for any $\bm{\theta}$, $\mathbb{P}(d(\mathcal{D},\mathcal{D}_{s})<\epsilon_{\ell}\mid\nobreak\bm{\theta})>0$, then $\mathbb{P}(d(\mathcal{D},\mathcal{D}_{s})<\epsilon_{\ell-1}\mid\nobreak\bm{\theta})>0$ since $\epsilon_{\ell}<\epsilon_{\ell-1}$. That is, any simulated data generated using parameter values taken from outside $\operatorname*{supp}(p(\bm{\theta}_{\epsilon_{\ell-1}}))$ cannot be accepted on level $\ell$ since $d(\mathcal{D},\mathcal{D}_{s})>\epsilon_{\ell-1}$ almost surely. Therefore, we can truncate the prior to the support of $p(\bm{\theta}_{\epsilon_{\ell-1}})$ when computing $\hat{Y}^{N_{\ell}}_{\ell}(\mathbf{s})$, thus increasing the acceptance rate of level $\ell$ samples. In practice, we need to approximate the support regions through sampling. For simplicity, in this work we restrict sampling of the prior at level $\ell$ to within the bounding box that contains all the samples generated at level $\ell-1$. However, more sophisticated approaches could be considered, and may result in further computational improvements.

We now have all the components to construct our MLMC-ABC algorithm. We compute the MLMC CDF estimator (Equation (8)) using the coupling technique for the bias correction terms (Section 2.3) and prior truncation for improved acceptance rates (Section 2.5).

Optimal $N_{1},\ldots,N_{L}$ are estimated as per Equation (10) and Section 2.4. Once $N_{1},\ldots,N_{L}$ have been estimated, computation of the MLMC-ABC posterior CDF $\hat{F}_{\epsilon_{L}}(\mathbf{s})$ proceeds according to Algorithm 4.

The computational complexity of MLMC-ABC (Algorithm 4) is roughly $\mathcal{O}(c_{S}+N_{L}(c_{L}+c_{G}))$ where $c_{L}$ is the expected cost of generating a single sample of the ABC posterior with threshold $\epsilon_{L}$, $c_{S}$ is the cost of updating the CDF estimate and $c_{G}$ is the cost of the coupling which involves the marginal CDF inverses. From Algorithm 4, we have $c_{S}=\mathcal{O}(N_{1}|\mathcal{S}|)$ and from Barber et al. (2015) we have $c_{L}=\mathcal{O}(\epsilon^{-n})$ where $n$ is the dimensionality of $\mathcal{D}$. The marginal inverse operations require only two steps:

1. 1.

   find $\mathbf{s}\in\mathcal{S}$ such that, $s_{j}\leq\theta_{\epsilon_{\ell},j}^{i}<s_{j}+\delta_{j}$. Such an operation is, at most, $\mathcal{O}(\log_{k}|\mathcal{S}|)$;

2. 2.

   invert the interpolating cubic spline, which can be done in $\mathcal{O}(1)$.

It follows that $c_{G}=\mathcal{O}(\log_{k}|\mathcal{S}|)$. For any practical application, $\epsilon_{L}$ will need to be sufficiently small, that is, $(N_{1}|\mathcal{S}|/N_{L}-\log_{k}|\mathcal{S}|)\ll c_{L}$, in order for the cost of generating posterior samples at level $L$ dominates the cost of the marginal CDF inverse operations and the lattice updating. Computational gains over ABC rejection sampling are achieved through decreasing $N_{L}$ and $c_{L}$, via variance reduction and prior truncation.

Our primary focus in Algorithm 4 is on using MLMC to estimate the posterior CDF. The coupling mechanism is more readily communicated in this case. However, the MLMC-ABC method is more general and can be used to estimate $\mathbb{E}\left[U(\bm{\theta}_{\epsilon_{\ell}})\right]$ where $U(\bm{\theta}_{\epsilon_{\ell}})$ is any Lipschitz continuous function. In this more general case, only the marginal CDFs need be accumulated to facilitate the coupling mechanism. For more details see A.

The SIS model is a common model from epidemiology that describes the spread of a disease or infection for which no significant immunity is obtained after recovery; the common cold, for example. The model is given by

with parameters $\bm{\theta}=\{\beta,\gamma\}$ and hazard functions for infection and recovery given by

respectively. This process defines a discrete-state continuous-time Markov process with a forward transitional density function that is computationally feasible to evaluate exactly for small populations $N_{pop}=S+I$. That is, for $t>s$, the probability of $S(t)=x$ given $S(s)=y$, with density denoted by $p(x,t\mid\nobreak y,s;\beta,\gamma)$, has a solution obtained through the $x$th element of the vector

where the $y$th element of column vector, $\mathbf{y}$, is one and all other elements zero, $\exp(\cdot)$ denotes the matrix exponential and $Q(\beta,\gamma)$ is the infinitesimal generator matrix of the Markov process; for the SIS model, $Q(\beta,\gamma)$ is a tri-diagonal matrix dependent only on the parameters of the model.

Let $S_{obs}(t)$ be an observation at time $t$ of the number of susceptible individuals in the population. We generate observations using a single realisation of the SIS model with parameters $\beta=0.003$ and $\gamma=0.1$, population size $N_{pop}=101$, and initial conditions $S(0)=100$ and $I(0)=1$; Observations are taken at times $t_{1}=4,\,t_{2}=8,\ldots,\,t_{10}=40$. Using the analytic solution to the SIS transitional density (using Equation (11)), we arrive at the likelihood function

where $S_{obs}(t_{0})=s_{0}$ almost surely. Hence, we can obtain an exact solution to the SIS posterior $p(\beta,\gamma\mid\nobreak\mathcal{D})$ given the priors $\beta\sim\mathcal{U}(0,0.06)$ and $\gamma\sim\mathcal{U}(0,2)$. Given this exact posterior, quadrature can be applied to compute the posterior CDF, given by

For the ABC approximation we use a discrepancy metric based on the sum of squared errors,

where $S^{*}(t)$ is a realisation of the model generated using the Gillespie algorithm (Gillespie, 1977) for a given set of parameters, and $\mathcal{D}_{s}=[S^{*}(t_{1}),S^{*}(t_{2}),\ldots,S^{*}(t_{10})]$. An appropriate acceptance threshold sequence for this metric is $\epsilon_{1},\ldots,\epsilon_{L}$, with $\epsilon_{\ell}=\epsilon_{1}m^{1-\ell}$, $m=2$ and $\epsilon_{1}=75$ (Toni et al., 2009).

We can use the exact likelihood (Equation (12)) to evaluate the convergence properties of our MLMC-ABC method. Based on the analysis of ABC rejection sampling by Barber et al. (2015), the rate of decay of the root mean-squared error (RMSE) as the computational cost is increased is slower for ABC than standard Monte Carlo. We find experimentally for ABC rejection sampling, that the decay of RMSE of the SIS model is approximately $\mathcal{O}(C^{-1/4})$, where $C$ is the expected computational cost of generating $N$ ABC posterior samples. The $95\%$ confidence interval of this experimental rate is $[-0.27,-0.23]$, which is consistent with the expected theoretical result of $-0.25$ (Barber et al., 2015). This is slower than the expected $\mathcal{O}(C^{-1/2})$ decay typically achieved with standard Monte Carlo. We compare the ABC rejection sampler decay rate with that achieved by our MLMC-ABC methods.

We use the $L_{\infty}$ norm for the RMSE, that is,

where $F$ is the exact posterior CDF, evaluated directly from the likelihood (Equation (12)), and $\hat{F}_{\epsilon_{L}}$ is the Monte Carlo estimate of the CDF. A regular lattice that consists of $300\times 300$ points is used for this estimate. The computation required for Gillespie simulations completely dominates the added computation of updating the lattice. The RMSE is computed using $20$ independent MLMC-ABC CDF estimations for $L=1,\ldots,3$. The sequence of sample numbers, $N_{1},\ldots,N_{L}$, is computed using Equation (10) with target variance of $\mathcal{O}(\epsilon^{2})$ and $100$ trial samples. Figure 1 demonstrates the improved convergence rate over the ABC rejection sampler convergence rate. Based on the least-squares fit, the RMSE decay is approximately $\mathcal{O}(C^{-1/3})$. The $95\%$ confidence interval for the rate is $[-0.26,-0.34]$ which is consistent with theory from Giles et al. (2015) for the univariate situation. For smaller target RMSE, this results in an order of magnitude reduction in the computational cost.

We also consider the effect of the additional bias introduced through the coupling mechanism as presented in Section 2.3. We set the sample numbers at all levels to be $10^{4}$ to ensure the Monte Carlo error is negligible and compare the bias for different values of acceptance threshold scaling factor $m$. The bias is computed according to the $L_{\infty}$ norm, that is,

where $\hat{F}_{\epsilon_{L}}^{c}$ denotes the estimator computed according to Algorithm 4 and $\hat{F}_{\epsilon_{L}}^{u}$ denotes the estimator computed without any coupling. That is, $\hat{F}_{\epsilon_{L}}^{u}$ is computed using standard Monte Carlo to evaluate each term in the MLMC telescoping sum (Equation (6)). Note that, computationally, $\hat{F}_{\epsilon_{L}}^{u}$ will always be inferior to a standard Monte Carlo estimate. Figure 2 shows that, not only does the bias decay as $m$ decreases, but also the additional bias is well within the order of bias expected from the ABC approximation. This result, along with Figure 1, demonstrates that the reduction in estimator variance can dominate the increase in bias. Thus, compared with standard Monte Carlo, a significantly lower RMSE for the same computational effort is achieved with MLMC using this coupling strategy.

We now evaluate the performance of MLMC-ABC using the model developed by Tanaka et al. (2006) in the study of tuberculosis transmission rate parameters using DNA fingerprint data (Small et al., 1994). This model has been selected due to the availability of published comparative performance evaluations of MCMC-ABC and SMC-ABC (Tanaka et al., 2006; Sisson et al., 2007).

The model proposed by Tanaka et al. (2006) describes the occurrence of tuberculosis infections over time and the mutation of the bacterium responsible, Myobacterium tuberculosis. The number of infections, $I$, at time $t$ is

where $G_{t}$ is the number of distinct genotypes and $X_{i,t}$ is the number of infections caused by the $i$th genotype at time $t$. For each genotype, new infections occur with rate $\alpha$, infections terminate with rate $\delta$, and mutation occurs with rate $\mu$; causing an increase in the number of genotypes. This process, as with the SIS model, can be described by a discrete-state continuous-time Markov process. In this case, however, the likelihood is intractable, but the model can still be simulated using the Gillespie algorithm (Gillespie, 1977). After a realisation of the model is completed, either by extinction or when a maximum infection count is reached, a sub-sample of $473$ cases is collected and compared against the IS$6110$ DNA fingerprint data of tuberculosis bacteria samples (Small et al., 1994). The dataset consists of $326$ distinct genotypes; the infection cases are clustered according to the genotype responsible for the infection. The collection of clusters can be summarised succinctly as $30^{1}\,23^{1}\,15^{1}\,10^{1}\,8^{1}\,5^{2}\,4^{4}\,3^{13}\,2^{20}\,1^{282},$ where $n^{j}$ denotes there are $j$ clusters of size $n$. The discrepancy metric used is

where $n$ is the size of the population sub-sample (e.g., $n=473$), $g(\mathcal{D})$ denotes the number of distinct genotypes in the dataset (e.g., $g(\mathcal{D})=326$) and the genetic diversity is

where $n_{i}(\mathcal{D})$ is the cluster size of the $i$th genotype in the dataset.

We perform likelihood-free inference on the tuberculosis model for the parameters $\bm{\theta}=\left\{\alpha,\delta,\mu\right\}$ with the goal of evaluating the efficiencies of MLMC-ABC, MCMC-ABC and SMC-ABC. We use a target posterior distribution of $p(\bm{\theta}\mid\nobreak d(\mathcal{D},\mathcal{D}_{s})<\epsilon)$ with $d(\mathcal{D},\mathcal{D}_{s})$ as defined in Equation (13) and $\epsilon=0.0025$. The acceptance threshold sequence, $\epsilon_{1},\ldots,\epsilon_{10}$, used for both SMC-ABC and MLMC-ABC is $\epsilon_{i}=\epsilon_{10}+(\epsilon_{i-1}-\epsilon_{10})/2$ with $\epsilon_{1}=1$ and $\epsilon_{10}=0.0025$. The improper prior is given by $\alpha\sim\mathcal{U}(0,5)$, $\delta\sim\mathcal{U}(0,\alpha)$ and $\mu\sim\mathcal{N}(0.198,0.06735^{2})$ (Sisson et al., 2007; Tanaka et al., 2006). For the MCMC-ABC and SMC-ABC algorithms we apply a typical Gaussian proposal kernel,