Asymptotics of ABC

TL;DR

This paper reviews recent theoretical insights into the asymptotic behavior of Approximate Bayesian Computation (ABC), highlighting its strengths in point estimation and challenges in uncertainty quantification, with practical implementation guidance.

Contribution

It provides an informal review of recent asymptotic results for ABC, emphasizing the impact of different implementations on uncertainty estimation.

Findings

ABC performs well in point estimation with more data.

Standard ABC overestimates parameter uncertainty.

Regression correction improves uncertainty quantification.

Abstract

We present an informal review of recent work on the asymptotics of Approximate Bayesian Computation (ABC). In particular we focus on how does the ABC posterior, or point estimates obtained by ABC, behave in the limit as we have more data? The results we review show that ABC can perform well in terms of point estimation, but standard implementations will over-estimate the uncertainty about the parameters. If we use the regression correction of Beaumont et al. then ABC can also accurately quantify this uncertainty. The theoretical results also have practical implications for how to implement ABC.

Equations

$${\pi }_{ABC}(\mathbf{\theta })\propto p(\mathbf{\theta })\int f_n(\mathbf{s}|\mathbf{\theta })K\left(\frac{{\mathbf{s}}_{obs}-\mathbf{s}}{ϵ}\right)\text{mbox}d\mathbf{s}.$$

$$h_{ABC}=\int h(\mathbf{\theta }){\pi }_{ABC}(\mathbf{\theta })\text{mbox}d\mathbf{\theta }.$$

$$\text{mbox}{Pr}_{ABC}(\Vert \mathbf{\theta }-{\mathbf{\theta }}_0\Vert <\delta )={\int }_{\mathbf{\theta }:\Vert \mathbf{\theta }-{\mathbf{\theta }}_0\Vert <\delta }{\pi }_{ABC}(\mathbf{\theta })\text{mbox}d\mathbf{\theta },$$

$$\text{mbox}{Pr}_{ABC}(\Vert \mathbf{\theta }-{\mathbf{\theta }}_0\Vert <\delta )\to 1$$

$${\mathbf{S}}_{n,\mathbf{\theta }}\to \mathbf{b}(\mathbf{\theta })$$

$$\text{mbox}{Pr}_{ABC}(\Vert \mathbf{\theta }-{\mathbf{\theta }}_0\Vert <{\lambda }_n)\to 1$$

$$\sqrt{n}\left\{{\mathbf{S}}_{n,\mathbf{\theta }}-\mathbf{b}(\mathbf{\theta })\right\}\to \text{mbox}N\left\{0,A(\mathbf{\theta })\right\},$$

$${\mathbf{S}}_{n,\mathbf{\theta }}=\mathbf{b}({\mathbf{\theta }}_0)+D_0(\mathbf{\theta }-{\mathbf{\theta }}_0)+\frac{1}{\sqrt{n}}\mathbf{Z},$$

$${\mathbf{S}}_{n,\mathbf{\theta }}=\mathbf{b}({\mathbf{\theta }}_0)+D_0(\mathbf{\theta }-{\mathbf{\theta }}_0)+\frac{1}{\sqrt{n}}\mathbf{Z}+{ϵ}_n\mathbf{V}.$$

$${\tilde{\mathbf{S}}}_{n,\mathbf{\theta }}=D_0\mathbf{t}+\mathbf{V}+\frac{1}{{ϵ}_n\sqrt{n}}\mathbf{Z}\to D_0\mathbf{t}+\mathbf{V}.$$

$${\tilde{\mathbf{S}}}_{n,\mathbf{\theta }}=D_0\mathbf{t}+\mathbf{Z}+{ϵ}_n\sqrt{n}\mathbf{V}\to D_0\mathbf{t}+\mathbf{Z}+c\mathbf{V},$$

$$\left\{D_0^{T}A({\theta }_0{)}^{-1}{D}_0\right){\}}^{-1}{D}_0^{T}A({\theta }_0{)}^{-1}{\tilde{\mathbf{s}}}_{obs},$$

$${\left\{D_0^T(A({\theta }_0)+c^2\Sigma {)}^{-1}{D}_0\right\}}^{-1}{D}_0^T\{A({\theta }_0)+c^2\Sigma {\}}^{-1}{\tilde{\mathbf{s}}}_{obs}$$

$$\tilde{I}=D_0^T\{A({\theta }_0)+c^2\Sigma {\}}^{-1}{D}_0.$$

$${\left(D_0^{T}A({\theta }_0{)}^{-1}{D}_0\right)}^{-1}{D}_0^{T}A({\theta }_0{)}^{-1}{\tilde{\mathbf{s}}}_{obs}.$$

$${\hat{\mathbf{\theta }}}_N=\frac{1}{{\sum }_{k=1}^{N_{acc}}{w}_k}\sum _{k=1}^{N_{acc}}{w}_k{\mathbf{\theta }}^k.$$

$$\Vert {\mathbf{s}}_{obs}-{\mathbf{s}}_i\Vert =\Vert \{{\mathbf{s}}_{obs}-\mathbf{b}({\mathbf{\theta }}_0)\}+\{\mathbf{b}({\mathbf{\theta }}_0)-\mathbf{b}({\mathbf{\theta }}_i)\}+\{\mathbf{b}({\mathbf{\theta }}_i)-{\mathbf{s}}_i\}\Vert .$$

$${\tilde{\mathbf{\theta }}}^k={\mathbf{\theta }}^k-\hat{B}({\mathbf{s}}^k-{\mathbf{s}}_{obs}).$$

$$\sum _{i=1}^p{\left({\mathbf{\theta }}_i-{\mathbf{\alpha }}_i-\sum _{j=1}^d{B}_{ij}{\mathbf{S}}_j\right)}^2$$

Full text

Asymptotics of ABC

Paul Fearnhead

Department of Mathematics and Statistics, Lancaster University

Abstract

This document is due to appear as a chapter of the forthcoming Handbook of Approximate Bayesian Computation (ABC) edited by S. Sisson, Y. Fan, and M. Beaumont.

We present an informal review of recent work on the asymptotics of Approximate Bayesian Computation (ABC). In particular we focus on how does the ABC posterior, or point estimates obtained by ABC, behave in the limit as we have more data? The results we review show that ABC can perform well in terms of point estimation, but standard implementations will over-estimate the uncertainty about the parameters. If we use the regression correction of Beaumont et al. then ABC can also accurately quantify this uncertainty. The theoretical results also have practical implications for how to implement ABC.

1 Introduction

This chapter aims to give an overview of recent work on the asymptotics of Approximate Bayesian Computation (ABC). By asymptotics here we mean how does the ABC posterior, or point estimates obtained by ABC, behave in the limit as we have more data? The chapter summarises results from three papers, Li and Fearnhead (2015), Frazier et al. (2016) and Li and Fearnhead (2016). The presentation in this chapter is deliberately informal, with the hope of conveying both the intuition behind the theoretical results from these papers and the practical consequences of this theory. As such we will not present all the technical conditions for the results we give: the interested reader should consult the relevant papers for these, and the results we state should be interpreted as holding under appropriate regularity conditions.

We will focus on ABC for a $p$-dimensional parameter, $\boldsymbol{\theta}$, from a prior $p(\boldsymbol{\theta})$ (we use the common convention of denoting vectors in bold, and we will assume these are column vectors). We assume we have data of size $n$ that is summarised through a $d$-dimensional summary statistic. The asymptotic results we review consider the limit $n\rightarrow\infty$, but assume that the summary statistic is of fixed dimension. Furthermore all results assume that the dimension of the summary statistic is at least as large as the dimension of the parameters, $d\geq p$ – this is implicit in the identifiability conditions that we will introduce later. Examples of such a setting are where the summaries are sample means of functions of individual data points, quantiles of the data, or, for time-series data, are empirical auto-correlations of the data. It also includes summaries based on fixed-dimensional auxillary models Drovandi et al. (2015) or on composite likelihood score functions Ruli et al. (2016).

To distinguish the summary statistic for the observed data from the summary statistic of data simulated within ABC, we will denote the former by $\boldsymbol{s}_{obs}$, and the latter by $\boldsymbol{s}$. Our model for the data will define a probability model for the summary. We assume that this in turn specifies a probability density function, or likelihood, for the summary, $f_{n}(\boldsymbol{s};\boldsymbol{\theta})$, which depends on the parameter. In some situations we will want to refer to the random variable for the summary statistic, and this will be $\boldsymbol{S}_{n,\boldsymbol{\theta}}$. As is standard with ABC, we assume that we can simulate from the model but cannot calculate $f_{n}(\boldsymbol{s};\boldsymbol{\theta})$.

The most basic ABC algorithm is a rejection sampler Pritchard et al. (1999) , which iterates the following three steps:

• (RS1)

Simulate a parameter from the prior: $\theta_{i}\sim p(\boldsymbol{\theta})$.

• (RS2)

Simulate a summary statistic from the model given $\boldsymbol{\theta}_{i}$: $\boldsymbol{s}_{i}\sim f_{n}(\boldsymbol{s}|\boldsymbol{\theta}_{i})$.

• (RS3)

Accept $\boldsymbol{\theta}_{i}$ if $\lVert\boldsymbol{s}_{obs}-\boldsymbol{s}_{i}\rVert<\epsilon$.

Here $\lVert\boldsymbol{s}_{obs}-\boldsymbol{s}_{i}\rVert$ is a suitably chosen distance between the observed and simulated summary statistics, and $\epsilon$ is a suitably chosen bandwidth. In the following we will assume that $||\boldsymbol{x}||$ is either Euclidean distance, $||\boldsymbol{x}||^{2}=\boldsymbol{x}^{T}\boldsymbol{x}$, or a Mahalanobis distance, $||\boldsymbol{x}||^{2}=\boldsymbol{x}^{T}\Gamma\boldsymbol{x}$ for some chosen positive-definite $d\times d$ matrix $\Gamma$.

If we define a (uniform) kernel function, $K(\boldsymbol{x})$, to be 1 if $\lVert\boldsymbol{x}\rVert<1$ and 0 otherwise, then this rejection sampler is drawing from the following distribution

[TABLE]

We call this the ABC posterior . If we are interested in estimating a function of the parameter $h(\boldsymbol{\theta})$ we can use the ABC posterior mean

[TABLE]

In practice we cannot calculate this posterior mean analytically, but would have to estimate it based on the sample mean of $h(\boldsymbol{\theta}_{i})$ for parameter values $\boldsymbol{\theta}_{i}$ simulated using the above rejection sampler.

In this chapter we review results on the behaviour of the ABC posterior, the ABC posterior mean, and Monte Carlo estimates of this mean as $n\rightarrow\infty$. In particular we consider whether the ABC posterior concentrates around the true parameter value in Section 2. We then consider the limiting form of the ABC posterior and the frequentist asymptotic distribution of the ABC posterior mean in Section 3. For the latter two results we compare these asymptotic distributions with those of the true posterior given the summary – which is the best we can hope for once we have chosen our summary statistics.

The results in these two sections ignore any Monte Carlo error. The impact of Monte Carlo error on the asymptotic variance of our ABC posterior mean estimate is the focus of Section 4. This impact depends on the choice of algorithm we use to sample from the ABC posterior (whereas the choice of algorithm has no effect on the actual ABC posterior or posterior mean that are analysed in the earlier sections). The rejection sampling algorithm above is inefficient in the limit as $n\rightarrow\infty$ and thus we consider more efficient importance sampling and MCMC generalisations in this section.

We then review results that show how post-processing the output of ABC can lead to substantially stronger asymptotic results. The chapter then finishes with a discussion that aims to draw out the key practical insights from the theory.

Before we review these results, it is worth mentioning that we can generalise the definition of the ABC posterior, and the associate posterior mean, given above. Namely we can use a more general form of kernel than the uniform kernel. Most of the results we review apply if we replace the uniform kernel by a different kernel, $K(\boldsymbol{x})$, that is monotonically decreasing in $\lVert\boldsymbol{x}\rVert$. Furthermore the specific form of the kernel has little affect on the asymptotic results – what matters most is how we choose the bandwidth and, in some cases, the choice of distance. The fact that most of the theoretical results do not depend on the choice of kernel means that, for concreteness, we will primarily assume a uniform kernel in our presentation below. The exceptions being in Section 3 where it is easier to get an intuition for the results if we use a Gaussian kernel. By focussing on these two choices we do not mean to suggest that they are necessarily better than other choices, it is just that they simplify the exposition. We will return to the choice of kernel in the Discussion.

2 Posterior Concentration

The results we present in this section are from Frazier et al. (2016) (though see also Martin et al., 2016), and consider the question of whether the ABC posterior will place increasing probability mass around the true parameter value as $n\rightarrow\infty$. It is the most basic convergence result we would wish for, requires weaker conditions than results we give in Section 3, and is thus easier to apply to other ABC settings (see for example Marin et al., 2014; Bernton et al., 2017).

We will denote the true parameter value by $\boldsymbol{\theta}_{0}$. If we define

[TABLE]

the ABC posterior probability that $\boldsymbol{\theta}$ is within some distance $\delta$ of the true parameter value, then for posterior concentration we want that for any $\delta>0$

[TABLE]

as $n\rightarrow\infty$. That is, for any strictly positive choice of distance, $\delta$, regardless of how small it is, as $n\rightarrow\infty$ we need the ABC posterior to place all its probability on the event that $\boldsymbol{\theta}$ is within $\delta$ of the true parameter value.

To obtain posterior concentration for ABC we will need to let the bandwidth depend on $n$, and henceforth we denote the bandwidth by $\epsilon_{n}$.

2.1 ABC Posterior Concentration

The posterior concentration result of Frazier et al. (2016) is based upon assuming a law of large numbers for the summary statistics. Specifically we need the existence of a binding function, $\boldsymbol{b}(\boldsymbol{\theta})$, such that for any $\boldsymbol{\theta}$

[TABLE]

in probability as $n\rightarrow\infty$. If this holds, and the binding function satisfies an identifiability condition: that $\boldsymbol{b}(\boldsymbol{\theta})=\boldsymbol{b}(\boldsymbol{\theta}_{0})$ implies $\boldsymbol{\theta}=\boldsymbol{\theta}_{0}$, then we have posterior concentration providing the bandwidth tends to zero, $\epsilon_{n}\rightarrow 0$.

To gain some insight into this result and the assumptions behind it, we present an example. To be able to visuallise what is happening we will assume that the parameter and summary statistic are both 1-dimensional. Figure 1 shows an example binding function, a value of $\theta_{0}$ and $s_{obs}$, and output from the ABC rejection sampler.

As $n$ increases we can see the plotted points, that show proposed parameter and summary statistic values, converge towards the line that shows the binding function. This stems from our assumption of a law of large numbers for the summaries, so that for each $\theta$ value the summaries should tend to $b(\boldsymbol{\theta})$ as $n$ increases.

We also have that the observed summary statistic, $s_{obs}$, converges towards $b(\theta_{0})$. Furthermore we are decreasing the bandwidth as we increase $n$, which corresponds to narrower acceptance regions for the summaries, which means that the accepted summary statistics converge towards $b(\theta_{0})$. Asymptotically, only parameter values close to $\theta_{0}$, which have values $b(\theta)$ which are close to $b(\theta_{0})$, will simulate summaries close to $b(\theta_{0})$. Hence the only accepted parameter values will be close to, and asymptotically will concentrate on, $\theta_{0}$. This can be seen in practice from the plots in the bottom row of Figure 1.

The identifiability condition on the binding function is used to ensure that concentration of accepted summaries around $\boldsymbol{b}(\boldsymbol{\theta}_{0})$ results in ABC posterior concentration around $\boldsymbol{\theta}_{0}$. What happens when this identifiability condition does not hold is discussed in Section 2.3.

2.2 Rate of Concentration

We can obtain stronger results by looking at the rate at which concentration occurs. Informally we can think of this as the supremum of rates, $\lambda_{n}\rightarrow 0$, such that

[TABLE]

as $n\rightarrow\infty$. For parametric Bayesian inference with independent and identically distributed data this rate would be $1/\sqrt{n}$.

Assuming the binding function is continuous at $\boldsymbol{\theta}_{0}$, then the rate of concentration will be determined by the rate at which accepted summaries concentrate on $\boldsymbol{b}(\boldsymbol{\theta}_{0})$. As described above, this depends on the variability (or ‘noise’) of the simulated summaries around the binding function and on the bandwidth, $\epsilon_{n}$. The rate of concentration will be the slower of the rate at which the noise in the summary statistics and the rate at which $\epsilon_{n}$ tend to 0.

We can see this from the example in Figure 2, where we show output from the ABC rejection sampler for different values of $n$, but with $\epsilon_{n}$ tending to 0 at either a faster or slower rate than that of the noise in the summaries. For each regime the rate of concentration of both the accepted summaries and of the accepted parameter values is determined by the slower of the two rates.

2.3 Effect of Binding Function

The shape of the binding function for values of $\boldsymbol{\theta}$ for which $\boldsymbol{b}(\boldsymbol{\theta})$ is close to $\boldsymbol{b}(\boldsymbol{\theta}_{0})$ affects the ABC posterior as it affects the range of $\boldsymbol{\theta}$ values that will have a reasonable chance of producing summary statistic values that would be accepted by the ABC rejection sampler.

If the identifiability condition holds and the binding function is differentiable at $\boldsymbol{\theta}_{0}$ then the value of this gradient will directly impact the ABC posterior variance. This is shown in the top row of Figure 3. If this gradient is large (top-left plot) then even quite large differences in summary statistics would correspond to small differences in the parameter, and hence a small ABC posterior variance. By comparison if the gradient is small (top-right plot) then large differences in parameters may mean only small differences in summary statistics. In this case we expect a much larger ABC posterior variance for the same width of the region in which the summary statistics are accepted.

The bottom row of Figure 3 shows what can happen if the identifiability condition does not hold. The bottom-left plot gives an example where there are two distinct parameter values for which the binding function is equal to $b(\theta_{0})$. In this case we have a bi-modal ABC posterior that concentrates on these two values. The bottom-right plot shows an example where there is a range of parameter values whose binding function value is equal to $b(\theta_{0})$, and in this case the ABC posterior will concentrate on this range of parameter values.

It can be difficult in practice to know whether the identifiability condition holds. In large data settings, observing a multi-modal posterior as in the bottom-left plot of Figure 3 would suggest that it does not hold. In such cases it may be possible to obtain identifiability by adding extra summaries. The wish to ensure identifiability is one reason for choosing a higher dimensional summary than parameter. However this does not come without potential cost, as we show in Section 3.

2.4 Model Error

One of the implicit assumptions behind the result on posterior concentration is that our model is correct. This manifests itself within the assumption that as we get more data the observed summary statistic will converge to the value $\boldsymbol{b}(\boldsymbol{\theta}_{0})$. If the model we assume in ABC is incorrect then this may not be the case (see Frazier et al., 2017, for a fuller discussion of the impact of model error). There are then two possibilities, the first is that the observed summary statistic will converge to a value $\boldsymbol{b}(\tilde{\boldsymbol{\theta}})$ for some parameter value $\tilde{\boldsymbol{\theta}}\neq\boldsymbol{\theta}_{0}$. In this case, by the arguments above, we can still expect posterior concentration but to $\tilde{\boldsymbol{\theta}}$ and not $\boldsymbol{\theta}_{0}$.

The other possibility is that the observed summary statistic converges to a value that is not equal to $\boldsymbol{b}(\boldsymbol{\theta})$ for any $\boldsymbol{\theta}$. This is most likely to occur when the dimension of the summary statistic is greater than the dimension of the parameter. To give some insight into this scenario, we give in an example in Figure 4, where we have independent identically distributed data from a Gaussian distribution with mean $\theta$ and variance $\theta^{2}+2$, but our model assumes the mean and variance are $\theta$ and $\theta^{2}+1$ respectively. This corresponds to a wrong assumption about the variance. We then apply ABC with summary statistics that are the sample mean and variance.

As shown in the figure, we still can get posterior concentration in this setting. If we denote the limiting value of the binding function for the true model as $b_{0}$, then the posterior concentrates on parameter value, or values, whose binding function value is closest, according to the distance we use for deciding whether to accept simulated summaries, to $b_{0}$.

In this second scenario it may be possible to detect the model error by monitoring the closeness of the accepted summaries to the observed summaries. If the model is correct, then the distance between accepted and observed summaries tends to 0 with increasing $n$. Whereas in this second model error scenario, these distances will tend towards some non-zero constant.

3 ABC Posterior and Posterior Mean

We now consider stronger asymptotic results for ABC. To obtain these results we need extra assumptions in addition to those required for posterior concentration (see Frazier et al., 2016; Li and Fearnhead, 2015, for full details). The most important of these is that the summary statistics obey a central limit theorem

[TABLE]

for some $d\times d$ positive definite matrix $A(\boldsymbol{\theta})$. In the above central limit theorem we have assumed a $1/\sqrt{n}$ rate of convergence, but it is trivial to generalise this Li and Fearnhead (2015).

3.1 ABC Posterior

Under this central limit assumption we first consider convergence of the ABC posterior. Formal results can be found in Frazier et al. (2016) (but see also Li and Fearnhead, 2016). Here we give an informal presentation of these results.

To gain intuition about the limiting form of the ABC posterior, we can use the fact from the previous section that there is posterior concentration around $\boldsymbol{\theta}_{0}$. Thus asymptotically we need only consider the behaviour of the model for $\boldsymbol{\theta}$ close to $\boldsymbol{\theta}_{0}$. Also asymptotically the noise in the summaries is Gaussian. So if we make a linear approximation to $\boldsymbol{b}(\boldsymbol{\theta})$ for $\boldsymbol{\theta}$ close to $\boldsymbol{\theta}_{0}$, our model will be well approximated by

[TABLE]

where $D_{0}$ is the $d\times p$ matrix of first derivatives of $\boldsymbol{b}(\boldsymbol{\theta})$ with respect to $\boldsymbol{\theta}$, with these derivatives evaluated at $\boldsymbol{\theta}_{0}$; and $\boldsymbol{Z}$ is a $d$-dimensional Gaussian random variable with covariance matrix $A(\boldsymbol{\theta}_{0})$. Furthermore, for $\boldsymbol{\theta}$ close to $\boldsymbol{\theta}_{0}$ the prior will be well approximated by a uniform prior. For the following we assume that $D_{0}$ is of rank $p$.

Wilkinson (2013) shows that the effect of the approximation in ABC, whereby we accept simulated summaries which are similar, but not identical, to the observed summary, is equivalent to performing exact Bayesian inference under a different model. This different model has additional additive noise, where the distribution of the noise is given by the kernel, $K(\cdot)$, we use in ABC. So if $\boldsymbol{V}$ is a $d$-dimensional random variable with density $K(\cdot)$, independent of $\boldsymbol{Z}$, then our ABC posterior will behave like the true posterior for the model

[TABLE]

From Section 2.2, we know that the rate of concentration is the slower of the rate of the noise in the summaries, $1/\sqrt{n}$ under our central limit theorem, and the bandwidth $\epsilon_{n}$. This means that we get different limiting results depending on whether $\epsilon_{n}=O(1/\sqrt{n})$ or not. This can be seen from (1), as whether $\epsilon_{n}=O(1/\sqrt{n})$ or not will affect whether the $\epsilon_{n}\boldsymbol{V}$ noise term dominates or not.

If $\sqrt{n}\epsilon_{n}\rightarrow\infty$, so $\epsilon_{n}$ is the slower rate, then to get convergence of the ABC posterior we need to consider the re-scaled variable $\boldsymbol{t}=(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})/\epsilon_{n}$. If we further define $\tilde{\boldsymbol{S}}_{n,\boldsymbol{\theta}}=\{\boldsymbol{S}_{n,\boldsymbol{\theta}}-\boldsymbol{b}(\boldsymbol{\theta}_{0})\}/\epsilon_{n}$ then we can re-write (1) as

[TABLE]

Thus the limiting form of the ABC posterior is equivalent to the true posterior for this model, given observation $\tilde{\boldsymbol{s}}_{obs}=\{\boldsymbol{s}_{obs}-\boldsymbol{b}(\boldsymbol{\theta}_{0})\}/\epsilon_{n}$, with a uniform prior for $\boldsymbol{t}$. The shape of this posterior will be determined by the ABC kernel. If we use the standard uniform kernel, then the ABC posterior will asymptotically be uniform. By converting from $\boldsymbol{t}$ to $\boldsymbol{\theta}$ we see that the asymptotic variance for $\boldsymbol{\theta}$ is $O(1/\epsilon_{n}^{2})$ in this case.

The other case is that $\sqrt{n}\epsilon_{n}\rightarrow c$ for some positive, finite constant $c$. In this case we consider the re-scaled variable $\boldsymbol{t}=\sqrt{n}(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})$, and re-scaled observation $\tilde{\boldsymbol{S}}_{n,\boldsymbol{\theta}}=\sqrt{n}\{\boldsymbol{S}_{n,\boldsymbol{\theta}}-\boldsymbol{b}(\boldsymbol{\theta}_{0})\}$. The ABC posterior will asymptotically be equivalent to the true posterior for $\boldsymbol{t}$ under a uniform prior, for a model

[TABLE]

and given an observation $\tilde{\boldsymbol{s}}_{obs}=\sqrt{n}\{\boldsymbol{s}_{obs}-\boldsymbol{b}(\boldsymbol{\theta}_{0})\}$.

We make three observations from this. First if $\epsilon_{n}=o(1/\sqrt{n})$, so $c=0$, then using standard results for the posterior distribution of a linear model, the ABC posterior for $\boldsymbol{t}$ will converge to a Gaussian with mean

[TABLE]

and variance $I^{-1}$ where $I=D_{0}^{T}A(\boldsymbol{\theta}_{0})^{-1}D_{0}$. This is the same limiting form as the true posterior given the summaries. The matrix $I$ can be viewed as an information matrix, and note that this is larger if the derivatives of the binding function, $D_{0}$, are larger; in line with the intuition we presented in Section 2.3.

Second if $c\neq 0$, the ABC posterior will have a larger variance than the posterior given summaries. This inflation of the ABC posterior variance will increase as $c$ increases. In general it is hard to say the form of the posterior, as it will depend on the distribution of noise in our limiting model, $\boldsymbol{Z}+c\boldsymbol{V}$, which is a convolution of the limiting Gaussian noise of the summaries and a random variable drawn from the ABC kernel.

Our final observation is that we can get some insight into the behaviour of the ABC posterior when $c\neq 0$ if we assume a Gaussian kernel, as again the limiting ABC posterior will be the true posterior for a linear a model with Gaussian noise. If the Gaussian kernel has variance $\Sigma$, which corresponds to measuring distances between summary statistics using the scaled distance $\lVert\boldsymbol{x}\rVert=\boldsymbol{x}^{T}\Sigma^{-1}\boldsymbol{x}$, then the ABC posterior for $\boldsymbol{t}$ will converge to a Gaussian with mean

[TABLE]

and variance, $\tilde{I}^{-1}$, where

[TABLE]

3.2 ABC Posterior Mean

We now consider the asymptotic distribution of the ABC posterior mean . By this we mean the frequentist distribution, whereby we view the posterior mean as a function of the data, and look at the distribution of this under repeated sampling of the data. Formal results appear in Li and Fearnhead (2015), but we will give informal results, building on the results we gave for the ABC posterior. We will focus on the case where $\epsilon_{n}=O(1/\sqrt{n})$, but note that results hold for the situation where $\epsilon_{n}$ decays more slowly; in fact Li and Fearnhead (2015) show that if $\epsilon_{n}=o(n^{-3/10})$ then the ABC posterior mean will have the same asymptotic distribution as for the case we consider, where $\epsilon_{n}=O(1/\sqrt{n})$.

The results we stated for the ABC posterior in section 3.1 for the case $\epsilon_{n}=O(1/\sqrt{n})$ included expressions for the posterior mean; see (2) and (3). The latter expression was under the assumption of a Gaussian kernel in ABC, but most of the exposition we give below holds for a general kernel (see Li and Fearnhead, 2015, for more details).

The first of these, (2), is the true posterior mean given the summaries. Asymptotically our re-scaled observation $\tilde{\boldsymbol{s}}_{obs}$ has a Gaussian distribution with mean 0 and variance $A(\boldsymbol{\theta}_{0})$ due to the central limit theorem assumption, and the posterior mean for $\boldsymbol{t}$ is a linear transformation of $\tilde{\boldsymbol{s}}_{obs}$. This immediately gives that the asymptotic distribution of the ABC posterior mean of $\boldsymbol{t}$ is Gaussian with mean 0 and variance $I^{-1}$. Equivalently, for large $n$, the ABC posterior mean for $\boldsymbol{\theta}$ will be approximately normally distributed with mean $\boldsymbol{\theta}_{0}$ and variance $I^{-1}/n$.

The case where $\sqrt{n}\epsilon_{n}\rightarrow c$ for some $c>0$ is more interesting. If we have $d=p$, so we have the same number of summaries as we have parameters, then $D_{0}$ is a square matrix. Assuming this matrix is invertible, we see that the ABC posterior mean simplifies to $D_{0}^{-1}\tilde{\boldsymbol{s}}_{obs}$. Alternatively if $d>p$ but $\Sigma=\gamma A(\boldsymbol{\theta}_{0})$ for some scalar $\gamma>0$, so that the variance of our ABC kernel is proportional to the asymptotic variance of the noise in our summary statistics, then the ABC posterior mean again simplifies; this time to

[TABLE]

In both cases the expressions for the ABC posterior mean are the same as for the $c=0$ case, and are identical to the true posterior mean given the summaries. Thus the ABC posterior mean has the same limiting Gaussian distribution as the true posterior mean in these cases.

More generally for the $c>0$ case, the ABC posterior mean will be different from the true posterior mean given the summaries. In particular the asymptotic variance of the ABC posterior mean can be greater than the asymptotic variance of the true posterior mean given the summaries. Li and Fearnhead (2015) show that it is always possible to project a $d>p$ dimensional summary to a $p$ dimensional summary such that the asymptotic variance of the true posterior mean is not changed. This suggests using such a $p$ dimensional summary statistic for ABC (see Fearnhead and Prangle, 2012, for a different argument for choosing $d=p$). An alternative conclusion from these results is to scale the distance used when deciding whether to accept or reject summaries to be proportional an estimate of the variance of the noise in the summaries.

It is interesting to compare the asymptotic variance of the ABC posterior mean to the limiting value of the ABC posterior variance. Ideally these would be the same, as that implies that the ABC posterior is correctly quantifying uncertainty. We do get equality when $\epsilon_{n}=o(1/\sqrt{n})$; but in other cases we can see that the ABC posterior variance is larger than the asymptotic variance of the ABC posterior mean, and thus ABC over-estimates uncertainty. We will return to this in Section 5.

4 Monte Carlo Error

The previous section included results on the asymptotic variance of the ABC posterior mean – which gives a measure of accuracy of using the ABC posterior mean as a point estimate for the parameter. In practice we cannot calculate the ABC posterior mean analytically and we need to use output from a Monte Carlo algorithm, such as the rejection sampler described in the introduction. A natural question is what effect does the resulting Monte Carlo error have? And can we implement ABC in such a way that, for a fixed Monte Carlo sample size, the Monte Carlo estimate of the ABC posterior mean is an accurate point estimate? Or do we necessarily require the Monte Carlo sample size to increase as $n$ increases.

Li and Fearnhead (2015) explore these questions. To do so they consider an importance sampling version of the rejection sampling algorithm we previously introduced. This algorithm requires the specification of a proposal distribution for the parameter, $q(\boldsymbol{\theta})$, and involves iterating the following $N$ times

• (IS1)

Simulate a parameter from the proposal distribution: $\boldsymbol{\theta}_{i}\sim q(\boldsymbol{\theta})$.

• (IS2)

Simulate a summary statistic from the model given $\boldsymbol{\theta}_{i}$: $\boldsymbol{s}_{i}\sim f_{n}(\boldsymbol{s}|\boldsymbol{\theta}_{i})$.

• (IS3)

If $\lVert\boldsymbol{s}_{obs}-\boldsymbol{s}_{i}\rVert<\epsilon_{n}$ accept $\boldsymbol{\theta}_{i}$ and assign it a weight proportional to $\pi(\boldsymbol{\theta}_{i})/q(\boldsymbol{\theta}_{i})$.

The output is a set of, $N_{acc}$ say, weighted parameter values which can be used to estimate, for example, posterior means. With a slight abuse of notation, if the accepted parameter values are denoted $\boldsymbol{\theta}^{k}$ and their weights $w_{k}$ for $k=1,\ldots,N_{acc}$ then we would estimate the posterior mean of $\boldsymbol{\theta}$ by

[TABLE]

The use of this Monte Carlo estimator will inflate the error in our point estimate of the parameter by $\mbox{Var}(\hat{\boldsymbol{\theta}}_{N})$, where we calculate variance with respect to randomness of the Monte Carlo algorithm.

If the asymptotic variance of the ABC posterior mean is $O(1/n)$ we would want the Monte Carlo variance to be $O(1/(nN))$. This would mean that the overall impact of the Monte Carlo error is to inflate the mean square error of our estimator of the parameter by a factor $1+O(1/N)$ (similar to other likelihood free methods; e.g. Gourieroux et al., 1993; Heggland and Frigessi, 2004).

Now the best we can hope for with a rejection or importance sampler would be equally weighted, independent samples from the ABC posterior. The Monte Carlo variance of such an algorithm would be proportional to the ABC posterior variance. Thus if we want the Monte Carlo variance to be $O(1/n)$ then we need $\epsilon_{n}=O(1/\sqrt{n})$, as for slower rates the ABC posterior variance will decay more slowly than $O(1/n)$.

Thus we will focus on $\epsilon_{n}=O(1/\sqrt{n})$. The key limiting factor in terms of the Monte Carlo error of our rejection or importance sampler is the acceptance probability. To have a Monte Carlo variance that is $O(1/n)$ we will need an implementation whereby the acceptance probability is bounded away from 0 as $n$ increases. To see whether and how this is possible we can examine the acceptance criteria in step, (RS3) or (IS3):

[TABLE]

We need this distance to have a non-negligible probability of being less than $\epsilon_{n}$. Now the first and third bracketed terms on the right-hand side will be $O_{p}(1/\sqrt{n})$ under our assumption for the central limit theorem for the summaries. Thus this distance is at best $O_{p}(1/\sqrt{n})$, and if $\epsilon_{n}=o(1/\sqrt{n})$ the probability of the distance being less than $\epsilon_{n}$ should tend to 0 as $n$ increases.

This suggests we need $\sqrt{n}\epsilon_{n}\rightarrow c$ for some $c>0$. For this choice, if we have a proposal which has a reasonable probability of simulating $\boldsymbol{\theta}$ values within $O(1/\sqrt{n})$ of $\boldsymbol{\theta}_{0}$, then we could expect the distance to have a non-zero probability of being less than $\epsilon_{n}$ as $n$ increases. This rules out the rejection sampler, or any importance sampler with a pre-chosen proposal distribution. But an adaptive importance sampler that learns a good proposal distribution (e.g. Sisson et al., 2007; Beaumont et al., 2009; Peters et al., 2012) can have this property.

Note that such an importance sampler would need a proposal distribution for which the importance sampling weights are also well-behaved. Li and Fearnhead (2015) give a family a proposal distributions that have both an acceptance probability that is non-zero as $n\rightarrow\infty$ and have well-behaved importance sampling weights.

Whilst Li and Fearnhead (2015) did not consider MCMC based implementations of ABC Marjoram et al. (2003); Bortot et al. (2007), the intuition behind the results for the importance sampler suggest that we can implement such algorithms in a way that the Monte Carlo variance will be $O(1/(nN))$. For example if we use a random walk proposal distribution with a variance that is $O(1/n)$ then after convergence the proposed $\boldsymbol{\theta}$ values will be a distance $O_{p}(1/\sqrt{n})$ away from $\boldsymbol{\theta}_{0}$ as required. Thus the acceptance probability should be bounded away from 0 as $n$ increases. Furthermore such a scaling is appropriate for a random walk proposal to efficiently explore a target whose variance is $O(1/n)$ Roberts et al. (2001). Note that care would be needed whilst the MCMC algorithm is converging to stationarity as the proposed parameter values at this stage will be far away from $\boldsymbol{\theta}_{0}$.

5 The Benefits of Regression Adjustment

We finish this chapter by briefly reviewing asymptotic results for a popular version of ABC which post-processes the output of ABC using regression adjustment. This idea was first proposed by Beaumont et al. (2002) (see Nott et al., 2014, for links to Bayes linear methods). We will start with a brief description, then show how using regression adjustment can enable the adjusted ABC posterior to have the same asymptotic properties as the true posterior given the summaries, even if $\epsilon_{n}$ decays slightly slower than $1/\sqrt{n}$.

Figure 5 provides an example of the ABC adjustment. The idea is to run an ABC algorithm that accepts pairs of parameters and summaries. Denote these by $(\boldsymbol{\theta}^{k},\boldsymbol{s}^{k})$ for $k=1\ldots,N_{acc}$. These are shown in the top-left plot of Figure 5. We then fit $p$ linear models that, in turn, aim to predict each component of the parameter vector from the summaries. The output of this fitting procedure is a $p$-dimensional vector $\hat{\boldsymbol{\alpha}}$, the intercepts in the $p$ linear models, and a $p\times d$ matrix $\hat{B}$, whose $ij$th entry is the coefficient of the $j$ summary statistic in the linear model for estimating the $i$th component of $\boldsymbol{\theta}$.

An example of such fit is shown in the top-left hand plot of Figure 5. This fit is indicative of biases in our accepted $\boldsymbol{\theta}$ which correspond to different values of the summaries. In our example, the fit suggests that $\boldsymbol{\theta}$ values accepted for smaller, or larger, values of the summary statistic will, on average, be less then, or greater than, the true parameter value. We can then use the fit to correct for this bias. In particular we can adjust each of the accepted parameter values, to $\tilde{\boldsymbol{\theta}}^{k}$ for $k=1,\ldots,N_{acc}$ where

[TABLE]

The adjusted parameter values are shown in the bottom-left plot of Figure 5, and a comparison of the ABC posteriors before and after adjustment are shown in the bottom-right plot. From the latter we see the adjusted ABC posterior has a smaller variance and has more posterior mass close to the true parameter value.

The vector $\hat{\boldsymbol{\alpha}}$ and the matrix $\hat{B}$ can be viewed as estimates of the vector $\boldsymbol{\alpha}$ and the matrix $B$ that minimises the expectation of

[TABLE]

where expectation is with respect to parameter, summary statistic pairs drawn from our ABC algorithm. Li and Fearnhead (2016) show that if we adjust our ABC output using this optimal $B$ then, for any $\epsilon_{n}=o(n^{-3/10})$, the adjusted ABC posterior has the same asymptotic limit as the true posterior given the summaries. Obviously the asymptotic distribution of the mean of this adjusted posterior will also have the same asymptotic distribution as the mean of the true posterior given the summaries.

The intuition behind this result is that, asymptotically, if we choose $\epsilon_{n}=o(n^{-3/10})$, then our accepted samples will concentrate around the true parameter value. As we focus on an increasingly small ball around the true parameter value, the binding function will be well approximated by the linear regression model we are fitting. Thus the regression correction step is able to correct for the biases we obtain from accepting summaries that are slightly different from the observed summary statistics. From this intuition we see that a key requirement of our model, implicit within the assumptions needed for the theoretical result, is that the binding function is differentiable at the true parameter value: as such a differentiability condition is needed for the linear regression model to be accurate.

In practice we use an estimate $\hat{B}$, and this will inflate the asymptotic variance of the adjusted posterior mean by a factor that is $1+O(1/N_{acc})$, a similar effect to that of using Monte Carlo draws to estimate the mean. Importantly we get these strong asymptotic results even when $\epsilon_{n}$ decays more slowly than $1/\sqrt{n}$. For such a choice, for example $\epsilon_{n}=O(n^{-1/3})$, and with a good importance sampling or MCMC implementation, the asymptotic acceptance rate of the algorithm will tend to 1 as $n$ increases.

6 Discussion

The theoretical results we have reviewed are positive for ABC. If initially we ignore using regression adjustment, then the results suggest that ABC with $\epsilon_{n}=O(1/\sqrt{n})$ and with an efficient adaptive importance sampling or MCMC algorithm will have performance that is close to that of using the true posterior given the summaries. Ignoring Monte Carlo error, the accuracy of using the ABC posterior mean will be the same as that of using the true posterior mean if either we have the same number of summaries as parameters, or we choose an appropriate Mahalanobis distance for measuring the discrepancy in summary statistics. However, for this scenario the ABC posterior will over-estimate the uncertainty in our point estimate. The impact of Monte Carlo error will only be to inflate the asymptotic variance of our estimator by a factor $1+O(1/N)$, where $N$ is the Monte Carlo sample size.

We suggest that this scaling of the bandwidth, $\epsilon_{n}=O(1/\sqrt{n})$, is optimal if we do not use regression adjustment. Choosing either a faster or slower rate will result in Monte Carlo error that will dominate. One way of achieving this scaling is by using an adaptive importance sampling algorithm and fixing the proportion of samples to accept. Thus the theory supports the common practice of choosing the bandwidth indirectly in this manner.

Also based on these results, we suggest choosing the number of summary statistics to be close to, or equal to, the number of parameters, and choosing a distance for measuring the discrepancy in summary statistics that is based on the variance of the summary statistics. In situations where there are many potentially informative summary statistics then one of the many dimension reduction approaches, that try to construct low dimensional summaries that are information about the parameters, should be used (e.g. Wegmann et al., 2009; Fearnhead and Prangle, 2012; Blum et al., 2013; Prangle et al., 2014).

The results for ABC with regression adjustment are stronger still. These show that the ABC posterior and its mean can have the same asymptotics as the true ABC posterior and mean given the summaries. Furthermore this is possible with $\epsilon_{n}$ decreasing more slowly than $1/\sqrt{n}$, in which case the acceptance rate of a good ABC algorithm will increase as $n$ increases. These strong results suggest that regression adjustment should be routinely applied. One word of caution is that the regression adjustment involves fitting a number of linear-models to predict the parameters from the summaries. If a large number of summaries are used then the errors in fitting these models can be large Fearnhead and Prangle (2012) and lead to under-estimation of uncertainty in the adjusted posterior Marin et al. (2016). This again suggests using a small number of summary statistics, close or equal to the number of parameters.

Whilst the choice of bandwidth is crucial to the performance of ABC, and the choice of distance can also have an important impact on the asymptotic accuracy, the actual choice of kernel asymptotically has little impact. It affects the form of the ABC posterior, but does not affect the asymptotic variance of the ABC posterior mean (at least under relatively mild conditions).

These asymptotic results ignore any “higher-order” effects of the kernel that become negligible as $n$ gets large; so there may be some small advantages of one kernel over another for finite $n$, but these are hard to quantify. Intuitively the uniform kernel seems the most sensible choice – as for a fixed acceptance proportion it accepts the summaries closest to the observed. Furthermore in situations where there is model error it is natural to conjecture that a kernel with bounded support, such as the uniform kernel, will be optimal. For such a case we want to only accept summaries that are $d_{0}+O(1/\sqrt{n})$, for some constant distance $d_{0}>0$, away from the observed summary (see Figure 4). This is only possible for a kernel with bounded support.

Acknowledgements This work was supported by EPSRC through the i-like programme grant. It also benefitted from discussions during the BIRS workshop on Validating and Expanding ABC Methods in February 2017.