On the use of approximate Bayesian computation Markov chain Monte Carlo with inflated tolerance and post-correction

Approximate Bayesian computation is a form of likelihood-free inference (see, e.g., the reviews Marin et al., 2012; Sunnåker et al., 2013) which is used when exact Bayesian inference of a parameter $\theta\in\mathsf{T}$ with posterior density $\pi(\theta)\propto\mathrm{pr}(\theta)L(\theta)$ is impossible, where $\mathrm{pr}(\theta)$ is the prior density and $L(\theta)=g(y^{*}\mid\theta)$ is an intractable likelihood with data $y^{*}\in\mathsf{Y}$. More specifically, when the generative model of observations $g(\,\cdot\,\mid\theta)$ cannot be evaluated, but allows for simulations, we may perform relatively straightforward approximate inference based on the following (pseudo-)posterior:

where $\epsilon>0$ is a ‘tolerance’ parameter, and $K_{\epsilon}:\mathsf{Y}^{2}\to[0,\infty)$ is a ‘kernel’ function, which is often taken as a simple cut-off $K_{\epsilon}(y,y^{*})=1\left(\|s(y)-s(y^{*})\|\leq\epsilon\right)$, where $s:\mathsf{Y}\to\mathbb{R}^{d}$ extracts a vector of summary statistics from the (pseudo) observations.

The summary statistics are often chosen based on the application at hand, and reflect what is relevant for the inference task; see also (Fearnhead and Prangle, 2012; Raynal et al., to appear). Because $L_{\epsilon}(\theta)$ may be regarded as a smoothed version of the true likelihood $g(y^{*}\mid\theta)$ using the kernel $K_{\epsilon}$, it is intuitive that using a too large $\epsilon$ may blur the likelihood and bias the inference. Therefore, it is generally desirable to use as small a tolerance $\epsilon>0$ as possible, but because the computational methods suffer from inefficiency with small $\epsilon$, the choice of tolerance level is difficult (cf. Bortot et al., 2007; Sisson and Fan, 2018; Tanaka et al., 2006).

We discuss simple post-processing procedure which allows for consideration of a range of values for the tolerance $\epsilon\leq\delta$, based on a single run of approximate Bayesian computation Markov chain Monte Carlo (Marjoram et al., 2003) with tolerance $\delta$. Such post-processing was suggested in (Wegmann et al., 2009) (in case of simple cut-off), and similar post-processing has been suggested also with regression adjustment (Beaumont et al., 2002) (in a rejection sampling context). The method, discussed further in Section 2, can be useful for two reasons: A range of tolerances $\epsilon\leq\delta$ may be routinely inspected, which can reveal excess bias in the pseudo-posterior $\pi_{\delta}$; and the Markov chain Monte Carlo inference may be implemented with sufficiently large $\delta$ to allow for good mixing.

Our contribution is two-fold. We suggest straightforward-to-calculate approximate confidence intervals for the posterior mean estimates calculated from the post-processing output, and discuss some theoretical properties related to it. We also introduce an adaptive approximate Bayesian computation Markov chain Monte Carlo which finds a balanced $\delta$ during burn-in, using acceptance rate as a proxy, and detail a convergence result for it.

For the rest of the paper, we assume that the kernel function in (1) has the form

where $d:\mathsf{Y}^{2}\to[0,\infty)$ is any ‘dissimilarity’ function and $\phi:[0,\infty)\to[0,1]$ is a non-increasing ‘cut-off’ function. Typically $d(y,y^{*})=\|s(y)-s(y^{*})\|$, where $s:\mathsf{Y}^{2}\to\mathbb{R}^{d}$ are the chosen summaries, and in case of the simple cut-off discussed in Section 1, $\phi(t)=\phi_{\mathrm{simple}}(t)=1\left(t\leq 1\right)$. We will implicitly assume that the pseudo-posterior $\pi_{\epsilon}$ given in (1) is well-defined for all $\epsilon>0$ of interest, that is, $c_{\epsilon}=\int\mathrm{pr}(\theta)L_{\epsilon}(\theta)\mathrm{d}\theta>0$.

The following summarises the approximate Bayesian computation Markov chain Monte Carlo algorithm of Marjoram et al. (2003), with proposal $q$ and tolerance $\delta>0$:

Algorithm 1 may be implemented by storing only $\Theta_{k}$ and the related distances $T_{k}=d(Y_{k},y^{*})$, and in what follows, we regard either $(\Theta_{k},Y_{k})_{k\geq 1}$ or $(\Theta_{k},T_{k})_{k\geq 1}$ as the output of Algorithm 1. In practice, the initial values $(\Theta_{0},Y_{0})$ should be taken as the state of the Algorithm 1 run for a number of initial ‘burn-in’ iterations. We also introduce an adaptive algorithm for parameter tuning later (Section 4).

It is possible to consider a variant of Algorithm 1 where many (possibly dependent) observations $\tilde{Y}_{k}^{(1)},\ldots,\tilde{Y}_{k}^{(m)}\sim g(\,\cdot\,\mid\tilde{\Theta}_{k})$ are simulated in each iteration, and an average of their kernel values is used in the accept-reject step (cf. Andrieu et al., 2018). We focus here in the case of single pseudo-observation per iteration, following the asymptotic efficiency result of Bornn et al. (2017), but remark that our method may be applied in a straightforward manner also with multiple observations.

Algorithm 11 in Appendix details how $E_{\delta,\epsilon}(f)$ and $S_{\delta,\epsilon}(f)$ can be calculated in $O(n\log n)$ time simultaneously for all $\epsilon\leq\delta$ in case of simple cut-off. The estimator $E_{\delta,\epsilon}(f)$ approximates $\mathbb{E}_{\pi_{\epsilon}}[f(\Theta)]$ and $S_{\delta,\epsilon}(f)$ may be used to construct a confidence interval; see Algorithm 5 below. Theorem 4 details consistency of $E_{\delta,\epsilon}(f)$, and relates $S_{\delta,\epsilon}(f)$ to the limiting variance, in case the following (well-known) condition ensuring a central limit theorem holds:

Proof of Theorem 4 is given in Appendix. Inspired by Theorem 4, we suggest to report the following approximate confidence intervals for the suggested estimators:

The confidence interval in Algorithm 5 is straightforward application of Theorem 4, except for using a common integrated autocorrelation estimate $\hat{\tau}_{\delta}(f)$ for all $\tau_{\delta,\epsilon}(f)$. This relies on the approximation $\tau_{\delta,\epsilon}(f)\lessapprox\tau_{\delta}(f)$, which may not always be entirely accurate, but likely to be reasonable, as illustrated by Theorem 6 in Section 3 below. We suggest using a common $\hat{\tau}_{\delta}(f)$ for all tolerances because direct estimation of integrated autocorrelation is computationally demanding, and likely to be unstable for small $\epsilon$.

The classical choice for $\hat{\tau}_{\delta}(f)$ in Algorithm 5(ii) is windowed autocorrelation, $\hat{\tau}_{\delta}(f)=\sum_{k=-\infty}^{\infty}\omega(k)\hat{\rho}_{k}$, with some $0\leq\omega(k)\leq 1$, where $\hat{\rho}_{k}$ is the sample autocorrelation of $\big{(}f(\Theta_{k})\big{)}$ (cf. Geyer, 1992). We employ this approach in our experiments with $\omega(k)=1\left(|k|\leq M\right)$ where the cut-off lag $M$ is chosen adaptively as the smallest integer such that $M\geq 5\big{(}1+2\sum_{i=1}^{M}\hat{\rho}_{k}\big{)}$ (Sokal, 1996). Also more sophisticated techniques for the calculation of the asymptotic variance have been suggested (e.g. Flegal and Jones, 2010).

We remark that, although we focus here on the case of using a common cut-off $\phi$ for both the abc-mcmc($\delta$) and the post-correction, one could also use a different cut-off $\phi_{s}$ in the simulation phase, as considered by Beaumont et al. (2002) in the regression context. The extension to Definition 2 is straightforward, setting $U_{k}^{(\delta,\epsilon)}=\phi(T_{k}/\epsilon)\big{/}\phi_{s}(T_{k}/\delta)$, and Theorem 4 remains valid under a support condition.

The following result, whose proof is given in Appendix, gives an expression for the integrated autocorrelation in case of simple cut-off.

We next discuss how this loosely suggests that $\tau_{\delta,\epsilon}(f)\lessapprox\tau_{\delta,\delta}(f)=\tau_{\delta}(f)$. The weight $\bar{w}_{\delta,\delta}\equiv 1$, and under suitable regularity conditions both $\bar{w}_{\delta,\epsilon}(\theta)$ and $\check{\tau}_{\delta,\epsilon}(f)$ are continuous with respect to $\epsilon$, and $\bar{w}_{\delta,\epsilon}(\theta)\to 0$ as $\epsilon\to 0$. Then, for $\epsilon\approx\delta$, we have $\bar{w}_{\delta,\epsilon}\approx 1$ and therefore $\tau_{\delta,\delta}(f)\approx\tau_{\delta,\epsilon}(f)$. For small $\epsilon$, the terms with $\mathrm{var}_{\pi_{\delta}}(f_{\delta,\epsilon})$ are of order $O(\bar{w}_{\delta,\epsilon}^{2})$, and are dominated by the other terms of order $O(\bar{w}_{\delta,\epsilon})$. The remaining ratio may be written as

where $\bar{g}_{\delta,\epsilon}\propto\{\bar{w}_{\delta,\epsilon}(1-\bar{w}_{\delta,\epsilon})\}^{1/2}f$ with $\pi_{\delta}(\bar{g}_{\delta,\epsilon}^{2})=1$. If $r_{\delta}(\theta)\leq r_{*}<1$, then the term is upper bounded by $2r_{*}(1-r_{*})^{-1}$, and we believe it to be often less than $\tau_{\delta,\delta}(f)$, because the latter expression is similar to the contribution of rejections to the integrated autocorrelation; see the proof of Theorem 6.

For general $\phi$, it appears to be hard to obtain similar theoretical result, but we expect the approximation to be still sensible. Theorem 6 relies on $Y_{k}^{(s)}$ being independent of $(\Theta_{k}^{(0)},Y_{k}^{(0)})$ conditional on $\Theta_{k}^{(s)}$, assuming at least single acceptance. This is not true with other cut-offs, but we believe that the dependence of $Y_{k}^{(s)}$ from $(\Theta_{0}^{(s)},Y_{0}^{(s)})$ given $\Theta_{k}^{(s)}$ is generally weaker than dependence of $\Theta_{k}^{(s)}$ and $\Theta_{0}^{(s)}$, suggesting similar behaviour.

We conclude the section with a general (albeit pessimistic) upper bound for the asymptotic variance of the post-corrected estimators.

Theorem 7 follows directly from (Franks and Vihola, 2017, Corollary 4). The upper bound guarantees that a moderate correction, that is, $\epsilon$ close to $\delta$ and $c_{\delta}$ close to $c_{\epsilon}$, is nearly as efficient as direct abc-mcmc($\delta$). Indeed, typically $w_{\delta,\epsilon}\to 1$ and $c_{\epsilon}\to c_{\delta}$ as $\epsilon\to\delta$, in which case Theorem 7 implies $\limsup_{\epsilon\to\delta}\sigma_{\delta,\epsilon}^{2}(f)\leq\sigma_{\epsilon}^{2}(f)$. However, as $\epsilon\to 0$, the bound becomes less informative.

We propose Algorithm 8 below to adapt the tolerance $\delta$ in abc-mcmc($\delta$) during a burn-in of length $n_{b}$, in order to obtain a user-specified overall acceptance rate $\alpha^{*}\in(0,1)$. Tolerance optimisation has been suggested earlier based on quantiles of distances, with parameters simulated from the prior (e.g. Beaumont et al., 2002; Wegmann et al., 2009). This heuristic might not be satisfactory in the Markov chain Monte Carlo context, if the prior is relatively uninformative. We believe that acceptance rate optimisation is a more natural alternative, and Sisson and Fan (2018) suggested this as well.

Our method requires also a sequence of decreasing positive step sizes $(\gamma_{k})_{k\geq 1}$. We used $\alpha^{*}=0.1$ and $\gamma_{k}=k^{-2/3}$ in our experiments, and discuss these choices later.

In practice, we use Algorithm 8 with a Gaussian symmetric random walk proposal $q_{\Sigma_{k}}$, where the covariance parameter $\Sigma_{k}$ is adapted simultaneously (Haario et al., 2001; Andrieu and Moulines, 2006) (see Algorithm 23 of Supplement D). We only detail theory for Algorithm 8, but note that similar simultaneous adaptation has been discussed earlier (cf. Andrieu and Thoms, 2008), and expect that our results could be elaborated accordingly.

The following conditions suffice for convergence of the adaptation:

Proof of Theorem 10 will follow from the more general Theorem 14 of Supplement B.

Polynomially decaying step size sequences as in Assumption 9 (i) are common in adaptation which is of the stochastic approximation type as our approach (Andrieu and Thoms, 2008). Slower decaying step sizes such as $n^{-2/3}$ often behave better with acceptance rate adaptation (cf. Vihola, 2012, Remark 3).

Simple random walk Metropolis with covariance adaptation (Haario et al., 2001) typically leads to a limiting acceptance rate around $0.234$ (Roberts et al., 1997). In case of a pseudo-marginal algorithm such as abc-mcmc($\delta$), the acceptance rate is lower than this, and decreases when $\delta$ is decreased (see Lemma 16 of Supplement C). Markov chain Monte Carlo would typically be necessary when rejection sampling is not possible, that is, when the prior is far from the posterior. In such a case, the likelihood approximation must be accurate enough to provide reasonable approximation $\pi_{\delta}\approx\pi_{\epsilon}$. This suggests that the desired acceptance rate should be taken substantially lower than $0.234$.

The choice of the desired acceptance rate $\alpha^{*}$ could also be motivated by theory developed for pseudo-marginal Markov chain Monte Carlo algorithms. Doucet et al. (2015) rely on log-normality of the likelihood estimators, which is problematic in our context, because the likelihood estimators take value zero. Sherlock et al. (2015) find the acceptance rate $0.07$ to be optimal under certain conditions, but also in a quite dissimilar context. Indeed, in our context, the $0.07$ guideline assumes a fixed tolerance, and informs about choosing the number of pseudo-data per iteration. As we stick with single pseudo-data per iteration following (Bornn et al., 2017), the $0.07$ guideline cannot be taken too informative. We recommend slightly higher $\alpha^{*}$ such as $0.1$ to ensure sufficient mixing.

Beaumont et al. (2002) suggested similar post-processing as in Section 2, applying a further regression correction. Namely, in the context of Section 2, we may consider a function $\tilde{f}^{(\epsilon)}(\theta,y)=f(\theta)-\bar{s}(y)^{\mathrm{\scriptscriptstyle T}}b_{\epsilon}$ where $\bar{s}(y)=s(y)-s(y^{*})$ and $b_{\epsilon}$ is a solution of

where $\tilde{\pi}_{\delta}$ is the stationary distribution of abc-mcmc($\delta$), with marginal $\pi_{\delta}$, given in Appendix. When the latter expectation is replaced by its empirical version, the solution coincides with weighted least squares $(\hat{a}_{\epsilon},\hat{b}_{\epsilon}^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}=(\mathrm{M}^{\mathrm{\scriptscriptstyle T}}\mathrm{W}_{\epsilon}\mathrm{M})^{-1}\mathrm{M}^{\mathrm{\scriptscriptstyle T}}\mathrm{W}_{\epsilon}v$, with $v_{k}=f(\Theta_{k})$, $\mathrm{W}_{\epsilon}=\mathrm{diag}(W_{1}^{(\delta,\epsilon)},\ldots,W_{n}^{(\delta,\epsilon)})$ and with matrix $\mathrm{M}$ having rows $[M]_{k,\,\cdot\,}=(1,\bar{s}(Y_{k})^{\mathrm{\scriptscriptstyle T}})$.

We suggest the following confidence interval for $a_{\epsilon}=E_{\tilde{\pi}_{\epsilon}}[\tilde{f}^{(\epsilon)}(\Theta,Y)]$ in the spirit of Algorithm 5:

where $\hat{\tau}^{\mathrm{reg}}_{\delta}$ is the integrated autocorrelation estimate for $(\hat{F}_{k}^{(\delta)})$ where $\hat{F}_{k}^{(\delta)}=f(\Theta_{k})-\bar{s}^{T}\hat{b}_{\delta}$ and $S_{\delta,\epsilon}^{\mathrm{reg}}=[(\mathrm{M}^{\mathrm{\scriptscriptstyle T}}\mathrm{W}_{\epsilon}M)^{-1}]_{1,1}\sum_{k=1}^{n}(W_{k}^{(\delta,\epsilon)})^{2}(\hat{F}_{k}^{(\epsilon)}-\hat{a}_{\epsilon})^{2}$, where the first term is included as an attempt to account for the increased uncertainty due to estimated $\hat{b}_{\epsilon}$, analogous to weighted least squares. Experimental results show some promise for this confidence interval, but we stress that we do not have better theoretical backing for it, and leave further elaboration of the confidence interval for future research.

We experiment with our methods on two models, a lightweight Gaussian toy example, and a Lotka-Volterra model. Our experiments focus on three aspects: can abc-mcmc($\delta$) with larger tolerance $\delta$ and post-correction to a desired tolerance $\epsilon<\delta$ deliver more accurate results than direct abc-mcmc($\epsilon$); does the approximate confidence interval appear reliable; how well does the tolerance adaptation work in practice. All the experiments are implemented in Julia (Bezanson et al., 2017), and the codes are available in https://bitbucket.org/mvihola/abc-mcmc.

Because we believe that Markov chain Monte Carlo is most useful when little is known about the posterior, we apply covariance adaptation (Haario et al., 2001; Andrieu and Moulines, 2006) throughout the simulation in all our experiments, using an identity covariance initially. When running the covariance adaptation alone, we employ the step size $n^{-1}$ as in the original method of Haario et al. (2001), and in case of tolerance adaptation, we use step size $n^{-2/3}$.

Regarding our first question, we investigate running abc-mcmc($\delta$) starting near the posterior mode with different pre-selected tolerances $\delta$. We first attempted to perform the experiments by initialising the chains from independent samples of the prior distribution, but in this case, most of the chains failed to accept a single move during the entire run. In contrast, our experiments with tolerance adaptation are initialised from the prior, and both the tolerances and the covariances are adjusted fully automatically by our algorithm.

We believe that approximate Bayesian computation inference with Markov chain Monte Carlo is a useful approach, when the chosen simulation tolerance allows for good mixing. Our confidence intervals for post-processing and automatic tuning of simulation tolerance may make this approach more appealing in practice.

A related approach by Bortot et al. (2007) makes tolerance an auxiliary variable with a user-specified prior. This approach avoids explicit tolerance selection, but the inference is based on a pseudo-posterior $\check{\pi}(\theta,\delta)$ not directly related to $\pi_{\delta}(\theta)$ in (1). Bortot et al. (2007) also provide tolerance-dependent analysis, showing parameter means and variances with respect to conditional distributions of $\check{\pi}(\theta,\delta)$ given $\delta\leq\epsilon$. We believe that our approach, where the effect of tolerance in the expectations with respect $\pi_{\epsilon}$ can be investigated explicitly, can be more immediate to interpret. Our confidence interval only shows the Monte Carlo uncertainty related to the posterior mean, and we are currently investigating how the overall parameter uncertainty could be summarised in a useful manner.

The convergence rates of approximate Bayesian computation has been investigated by Barber et al. (2015) in terms of cost and bias with respect to true posterior, and recently by Li and Fearnhead (2018a, b) in the large data limit, the latter in the context of regression. It would be interesting to consider extensions of these results in the Markov chain Monte Carlo context. In fact, Li and Fearnhead (2018a) already suggest that the acceptance rate must be lower bounded, which is in line with our adaptation rule.

Automatic selection of tolerance has been considered earlier in Ratmann et al. (2007), who propose an algorithm based on tempering and a cooling schedule. Based on our experiments, the tolerance adaptation we present in this paper appears to perform well in practice, and provides reliable results with post-correction. For the adaptation to work efficiently, the Markov chains must be taken relatively long, rendering the approach difficult for the most computationally demanding models.

We conclude with a brief discussion of certain extensions of the suggested post-correction method; more details are given in Supplement E. First, in case of non-simple cut-off, the rejected samples may be ‘recycled’ by using the acceptance probability as weight (Ceperley et al., 1977). The accuracy of the post-corrected estimator could be enhanced with smaller values of $\epsilon$ by performing further independent simulations from $g(\,\cdot\,\mid\Theta_{k})$ (which may be calculated in parallel). The estimator is rather straightforward, but requires some care because the estimators of the pseudo-likelihood take value zero. The latter extension, which involves additional simulations as post-processing, is similar to the ‘lazy’ version of Prangle (2016, 2015) incorporating a randomised stopping rule for simulation, and to debiased ‘exact’ approach of Tran and Kohn (2015), which may lead to estimators which get rid of $\epsilon$-bias entirely.

This work was supported by Academy of Finland (grants 274740, 284513 and 312605). The authors wish to acknowledge CSC, IT Center for Science, Finland, for computational resources, and thank Christophe Andrieu for useful discussions.