Asymptotics of ABC
Paul Fearnhead1,†
1

Department of Mathematics and Statistics, Lancaster University

arXiv:1706.07712v1 [stat.ME] 23 Jun 2017

†

Correspondence: p.fearnhead@lancaster.ac.uk

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, θ , from a prior p(θθ ) (we use the common
convention of denoting vectors in bold, and we will assume these are column vectors). We assume we
1

have data of size n that is summarised through a d-dimensional summary statistic. The asymptotic
results we review consider the limit n → ∞, 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 ≥ 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 s obs , and the latter by 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, fn (ss; θ ), 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 S n,θθ . As is standard with ABC, we assume that we can simulate from the model but cannot
calculate fn (ss; θ ).
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: θi ∼ p(θθ ).
(RS2) Simulate a summary statistic from the model given θ i : si ∼ fn (ss|θθ i ).
(RS3) Accept θ i if kssobs − s i k < .
Here kssobs −ssi k is a suitably chosen distance between the observed and simulated summary statistics,
x|| is either Euclidean
and  is a suitably chosen bandwidth. In the following we will assume that ||x
x||2 = x T x , or a Mahalanobis distance, ||x
x||2 = x T Γx
x for some chosen positive-definite
distance, ||x
d × d matrix Γ.
x), to be 1 if kx
xk < 1 and 0 otherwise, then this rejection
If we define a (uniform) kernel function, K(x
sampler is drawing from the following distribution


Z
sobs − s
dss.
πABC (θθ ) ∝ p(θθ ) fn (ss|θθ )K

We call this the ABC posterior . If we are interested in estimating a function of the parameter h(θθ )
we can use the ABC posterior mean
Z
hABC =

h(θθ )πABC (θθ )dθθ .

2

In practice we cannot calculate this posterior mean analytically, but would have to estimate it based
on the sample mean of h(θθ i ) for parameter values θ 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 → ∞. 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 → ∞ 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
x), that is monotonically decreasing in kx
xk. Furthermore the specific
kernel by a different kernel, K(x
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
3

mass around the true parameter value as n → ∞. 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 θ 0 . If we define
Z
PrABC (kθθ − θ 0 k < δ) =
πABC (θθ )dθθ ,
θ :kθθ −θθ 0 k<δ

the ABC posterior probability that θ is within some distance δ of the true parameter value, then
for posterior concentration we want that for any δ > 0
PrABC (kθθ − θ 0 k < δ) → 1
as n → ∞. That is, for any strictly positive choice of distance, δ, regardless of how small it is, as
n → ∞ we need the ABC posterior to place all its probability on the event that θ is within δ 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 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, b (θθ ),
such that for any θ
S n,θθ → b (θθ )
in probability as n → ∞. If this holds, and the binding function satisfies an identifiability condition:
that b (θθ ) = b (θθ 0 ) implies θ = θ 0 , then we have posterior concentration providing the bandwidth
tends to zero, n → 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 θ0 and sobs , 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 θ value the summaries should tend to
b(θθ ) as n increases.
We also have that the observed summary statistic, sobs , converges towards b(θ0 ). Furthermore we
are decreasing the bandwidth as we increase n, which corresponds to narrower acceptance regions
4

7
6
5

7
6
5

7
6
5
0.5

1.0

1.5

2.0

2.5

3.0

s

3.5

2

3

4

s
4
3
2

2

3

4

s

●

0.5

1.0

1.5

2.0

2.5

3.0

3.5

0.5

2.0

2.5

3.0

3.5

2.5

3.0

3.5

7
6

7
6

5
s
2

3

4

s
4
3
1.5

2.0

●

●

2
1.0

1.5

θ

5

6
5
2

3

4

s

●

0.5

1.0

θ

7

θ

0.5

1.0

1.5

2.0

θ

2.5

3.0

3.5

0.5

θ

1.0

1.5

2.0

2.5

3.0

3.5

θ

Figure 1: Example binding function, b(θ) (top-left plot). Pairs of parameter and summary statistic
values proposed by a rejection sampler (top-middle). Output of rejection sampler (top-right): θ0 and
b(θ0 ) (blue dotted vertical and horizontal lines respectively); sobs (bold red circle, and red dashed
horizontal line) and acceptance region for proposed summaries (bold red dashed horizonal lines);
pairs of parameter and summary statistic values accepted (bold) and rejected (grey) by the rejection
sampler. Bottom-row plots are the same as top-right plot but for increasing n and decreasing n .
Here, and for all plots, our results are for a simple scenario where data is IID Gaussian with a mean
that is a function of the parameter, and the summary statistic is the sample mean. (In this case the
binding function is, by definition, equal to the mean function.)
for the summaries, which means that the accepted summary statistics converge towards b(θ0 ).
Asymptotically, only parameter values close to θ0 , which have values b(θ) which are close to b(θ0 ),
will simulate summaries close to b(θ0 ). Hence the only accepted parameter values will be close to,
and asymptotically will concentrate on, θ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 b (θθ 0 ) results in ABC posterior concentration around θ 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, λn → 0, such that
PrABC (kθθ − θ 0 k < λn ) → 1

5

7
6

7
6

7

1.0

1.5

2.0

2.5

3.0

s

3.5

2

3

4

s
4
3
2
0.5

●

5

●

5

6
5
2

3

4

s

●

0.5

1.0

1.5

2.0

2.5

3.0

3.5

0.5

1.0

1.5

2.5

3.0

3.5

7
6

●

5

6
5

6
5

2.0
θ

7

θ

7

θ

●

0.5

1.0

1.5

2.0

2.5

3.0

3.5

s
2

3

4

s
4
3
2

2

3

4

s

●

0.5

1.0

1.5

2.0

θ

2.5

3.0

3.5

0.5

1.0

1.5

θ

2.0

2.5

3.0

3.5

θ

Figure 2: Example of ABC concentration for differing rates of the noise in the summary statistics
and rates of n . Plots are as in Figure 1. Top-row: noise in summary statistics halving, or equivalently
√
sample size increasing by a factor of 4, while n decreasing by 1/ 2 as we move from left to right.
√
Bottom-row: noise in summary statistics decreasing by 1/ 2, or equivalently sample size doubling,
while n halving as we move from left to right.
as n → ∞. For parametric Bayesian inference with independent and identically distributed data
√
this rate would be 1/ n.
Assuming the binding function is continuous at θ 0 , then the rate of concentration will be determined
by the rate at which accepted summaries concentrate on b (θθ 0 ). As described above, this depends
on the variability (or ‘noise’) of the simulated summaries around the binding function and on the
bandwidth, 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 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 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 θ for which b (θθ ) is close to b (θθ 0 ) affects the ABC
posterior as it affects the range of θ 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 θ 0 then the value
of this gradient will directly impact the ABC posterior variance. This is shown in the top row
6

7
6

7
6

5
s
4

4

s

5

●

3
2

2

3

●

0.5

1.0

1.5

2.0

2.5

3.0

3.5

0.5

1.0

1.5

2.0

2.5

3.0

3.5

2.5

3.0

3.5

θ

s
4

s

0.5

5

1.0

6

1.5

θ

●

2

−1.0

−0.5

3

0.0

●

0.5

1.0

1.5

2.0

2.5

3.0

3.5

0.5

θ

1.0

1.5

2.0
θ

Figure 3: Example of the effect of the shape of binding function on the ABC posterior (plots are
as in Figure 1). Top row: gradient of binding function at b(θ0 ) affects the ABC posterior variance,
with larger gradient (left-hand plot) resulting in lower ABC posterior variance than smaller gradient
(right-hand plot). Bottom row: effect of non-identifiability on ABC posterior.
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(θ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(θ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.

7

5
3

s2

4

5
4
3

s2

●

2
1

1

2

●

−2

−1

0

1

2

−2

−1

1

2

5
4
1

1

2

s2

●

3

4
3

●

2

s2

0
s1

5

s1

−2

−1

0

1

2

−2

s1

−1

0

1

2

s1

Figure 4: Example of the effect of model error in ABC for the Gaussian model with incorrect
variance described in the text. The plots, from left to right and top to bottom, correspond to
increasing sample size. Each plot shows the 2-dimensional binding function as we vary θ (line); the
observed summary statistic (red circle) and accepted (black dots) and rejected (grey dots) summary
statistic values. (For this model the parameter value used to simulate the summary statistics will
be close to the first summary statistic, s1 .)

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 b (θθ 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
b (θ̃θ ) for some parameter value θ̃θ 6= θ 0 . In this case, by the arguments above, we can still expect
posterior concentration but to θ̃θ and not θ 0 .
The other possibility is that the observed summary statistic converges to a value that is not equal
to b (θθ ) for any θ . 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 θ and variance θ2 + 2, but our model assumes the mean and variance are θ and θ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 b0 , then the posterior concentrates on

8

parameter value, or values, whose binding function value is closest, according to the distance we
use for deciding whether to accept simulated summaries, to b0 .
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

√

S n,θθ − b (θθ )} → N {0, A(θθ )} ,
n {S

for some d × d positive definite matrix A(θθ ). In the above central limit theorem we have assumed
√
a 1/ 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 θ 0 . Thus asymptotically we need only consider
the behaviour of the model for θ close to θ 0 . Also asymptotically the noise in the summaries is
Gaussian. So if we make a linear approximation to b (θθ ) for θ close to θ 0 , our model will be well
approximated by
1
S n,θθ = b (θθ 0 ) + D0 (θθ − θ 0 ) + √ Z ,
n
where D0 is the d × p matrix of first derivatives of b (θθ ) with respect to θ , with these derivatives
evaluated at θ 0 ; and Z is a d-dimensional Gaussian random variable with covariance matrix A(θθ 0 ).
Furthermore, for θ close to θ 0 the prior will be well approximated by a uniform prior. For the
following we assume that D0 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
9

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(·), we use in ABC. So if V is a ddimensional random variable with density K(·), independent of Z , then our ABC posterior will
behave like the true posterior for the model
1
S n,θθ = b (θθ 0 ) + D0 (θθ − θ 0 ) + √ Z + nV .
n

(1)

From Section 2.2, we know that the rate of concentration is the slower of the rate of the noise in
√
the summaries, 1/ n under our central limit theorem, and the bandwidth n . This means that we
√
get different limiting results depending on whether n = O(1/ n) or not. This can be seen from
√
(1), as whether n = O(1/ n) or not will affect whether the nV noise term dominates or not.
√
If nn → ∞, so n is the slower rate, then to get convergence of the ABC posterior we need to
S n,θθ = {S
S n,θθ − b (θθ 0 )}/n then we
consider the re-scaled variable t = (θθ − θ 0 )/n . If we further define S̃
can re-write (1) as
S n,θθ = D0t + V +
S̃

1
√ Z → D0t + V .
n n

Thus the limiting form of the ABC posterior is equivalent to the true posterior for this model, given
observation s̃sobs = {ssobs − b (θθ 0 )}/n , with a uniform prior for 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 t to θ we see that the asymptotic variance for θ is
O(1/2n ) in this case.
√

nn → c for some positive, finite constant c. In this case we consider the
√
√
S n,θθ = n{S
S n,θθ − b (θθ 0 )}. The ABC
re-scaled variable t = n(θθ − θ 0 ), and re-scaled observation S̃

The other case is that

posterior will asymptotically be equivalent to the true posterior for t under a uniform prior, for a
model

√
S n,θθ = D0t + Z + n nV
V → D0t + Z + cV
V,
S̃

and given an observation s̃sobs =

√

n{ssobs − b (θθ 0 )}.

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



D0T A(θ0 )−1 D0 }−1 D0T A(θ0 )−1s̃sobs ,

(2)

and variance I −1 where I = D0T A(θθ 0 )−1 D0 . 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, D0 , are larger; in line with the intuition we presented in
Section 2.3.

10

Second if c 6= 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,
V , which is a convolution of the limiting Gaussian noise of the summaries and a random
Z + cV
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 6= 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 Σ, which corresponds
xk = x T Σ−1x , then
to measuring distances between summary statistics using the scaled distance kx
the ABC posterior for t will converge to a Gaussian with mean


D0T (A(θ0 ) + c2 Σ)−1 D0

−1

D0T {A(θ0 ) + c2 Σ}−1s̃sobs

(3)

and variance, I˜−1 , where
I˜ = D0T {A(θ0 ) + c2 Σ}−1 D0 .

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 n = O(1/ n), but note that results hold for the situation
where n decays more slowly; in fact Li and Fearnhead (2015) show that if n = o(n−3/10 ) then the
ABC posterior mean will have the same asymptotic distribution as for the case we consider, where
√
n = O(1/ n).
√
The results we stated for the ABC posterior in section 3.1 for the case n = O(1/ 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 s̃sobs has a Gaussian distribution with mean 0 and variance A(θθ 0 ) due to the central
limit theorem assumption, and the posterior mean for t is a linear transformation of s̃sobs . This
immediately gives that the asymptotic distribution of the ABC posterior mean of t is Gaussian
with mean 0 and variance I −1 . Equivalently, for large n, the ABC posterior mean for θ will be
approximately normally distributed with mean θ 0 and variance I −1 /n.
11

The case where

√
nn → 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 D0 is a square matrix. Assuming this matrix
is invertible, we see that the ABC posterior mean simplifies to D0−1s̃sobs . Alternatively if d > p but
Σ = γA(θθ 0 ) for some scalar γ > 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
D0T A(θ0 )−1 D0

−1

D0T A(θ0 )−1s̃sobs .

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 n = o(1/ 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

12

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(θθ ), and involves iterating the following
N times
(IS1) Simulate a parameter from the proposal distribution: θ i ∼ q(θθ ).
(IS2) Simulate a summary statistic from the model given θ i : s i ∼ fn (ss|θθ i ).
(IS3) If kssobs − s i k < n accept θ i and assign it a weight proportional to π(θθ i )/q(θθ i ).
The output is a set of, Nacc 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 θ k and their weights wk for k = 1, . . . , Nacc then we would estimate the posterior mean of
θ by
1
θ̂θ N = PNacc
k=1

N
acc
X

wk

wkθ k .

k=1

The use of this Monte Carlo estimator will inflate the error in our point estimate of the parameter
by Var(θ̂θ 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 n = O(1/ n), as for slower rates the ABC posterior variance will decay more
slowly than O(1/n).
√
Thus we will focus on n = O(1/ 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):
kssobs − s i k = k{ssobs − b (θθ 0 )} + {bb(θθ 0 ) − b (θθ i )} + {bb(θθ i ) − s i }k.

13

We need this distance to have a non-negligible probability of being less than n . Now the first and
√
third bracketed terms on the right-hand side will be Op (1/ n) under our assumption for the central
√
√
limit theorem for the summaries. Thus this distance is at best Op (1/ n), and if n = o(1/ n) the
probability of the distance being less than n should tend to 0 as n increases.
√
This suggests we need nn → c for some c > 0. For this choice, if we have a proposal which
√
has a reasonable probability of simulating θ values within O(1/ n) of θ 0 , then we could expect
the distance to have a non-zero probability of being less than 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 → ∞ 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 θ values will be a distance Op (1/ n) away from θ 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 θ 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 n decays slightly slower than 1/ 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 (θθ k , s k ) for k = 1 . . . , Nacc . These are
14

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
α, the intercepts in the p linear models, and a p × d matrix B̂, whose
is a p-dimensional vector α̂
ijth entry is the coefficient of the j summary statistic in the linear model for estimating the ith
component of θ .
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 θ which correspond to different values of the summaries. In our example,
the fit suggests that θ 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
k

correct for this bias. In particular we can adjust each of the accepted parameter values, to θ̃θ for
k = 1, . . . , Nacc where
k

θ̃θ = θ k − B̂(ssk − s obs ).
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.
α and the matrix B̂ can be viewed as estimates of the vector α and the matrix B that
The vector α̂
minimises the expectation of
p
X

θ i − αi −

i=1

d
X

!2
Bij S j

j=1

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 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 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.
15

9
5

s

6

7

8

9
8
7
6
s
5

4
3
2

2

3

4

●

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

θ

0.8
0.6
0.4
0.0

2

3

0.2

4

5

s

6

ABC Posterior

7

1.0

8

1.2

9

θ

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

1.0

θ

1.5

2.0

2.5

3.0

3.5

θ

Figure 5: Example of the regression correction procedure of Beaumont et al. (2002) for a single
parameter, single summary statistic. Output of an ABC algorithm (top-left) showing accepted pairs
of parameter and summary values (dots), the binding function for this model (solid black line), and
θ0 and sobs (red circle and also blue vertical and red horizonal lines respectively). Top-right: the fit
from a linear model predicting the parameter value from the summary (blue solid line). Bottomleft: the adjusted output (black dots; with original output in grey); we plot both old and adjusted
parameter values against original summary statistic values. Bottom-right: the ABC posterior based
on the original accepted parameter values (black solid line) and the adjusted values (red dashed
line).

16

In practice we use an estimate B̂, and this will inflate the asymptotic variance of the adjusted
posterior mean by a factor that is 1 + O(1/Nacc ), a similar effect to that of using Monte Carlo draws
to estimate the mean. Importantly we get these strong asymptotic results even when n decays
√
more slowly than 1/ n. For such a choice, for example 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 n = O(1/ 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, n = O(1/ 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 n decreasing more slowly than 1/ n, in which case

17

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
√
d0 + O(1/ n), for some constant distance d0 > 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.

References
Beaumont, M. A., Zhang, W. and Balding, D. J. (2002). Approximate Bayesian computation in
population genetics. Genetics 162, 2025–2035.
Beaumont, M. A., Cornuet, J.-M., Marin, J.-M. and Robert, C. P. (2009). Adaptive approximate
Bayesian computation. Biometrika 96(4), 983–990.
Bernton, E., Jacob, P. E., Gerber, M. and Robert, C. P. (2017). Inference in generative models
using the Wasserstein distance. arXiv:1701.05146.
Blum, M. G., Nunes, M. A., Prangle, D., Sisson, S. A. et al. (2013). A comparative review of

18

dimension reduction methods in approximate Bayesian computation. Statistical Science 28(2),
189–208.
Bortot, P., Coles, S. G. and Sisson, S. A. (2007). Inference for stereological extremes. Journal of the
American Statistical Association 102(477), 84–92.
Drovandi, C. C., Pettitt, A. N., Lee, A. et al. (2015). Bayesian indirect inference using a parametric
auxiliary model. Statistical Science 30(1), 72–95.
Fearnhead, P. and Prangle, D. (2012). Constructing summary statistics for approximate Bayesian
computation: semi-automatic approximate Bayesian computation. Journal of the Royal Statistical
Society: Series B (Statistical Methodology) 74(3), 419–474.
Frazier, D. T., Martin, G. M., Robert, C. P. and Rousseau, J. (2016). Asymptotic Properties of
Approximate Bayesian Computation. arXiv.1607.06903.
Frazier, D. T., Robert, C. P. and Rousseau, J. (2017). Model misspecification in ABC: Consequences
and diagnostics. In preparation.
Gourieroux, C., Monfort, A. and Renault, E. (1993). Indirect inference. Journal of Applied Econometrics 8(S1), S85–S118.
Heggland, K. and Frigessi, A. (2004). Estimating functions in indirect inference. Journal of the
Royal Statistical Society: Series B 66, 447–462.
Li, W. and Fearnhead, P. (2015). On the asymptotic efficiency of ABC estimators. arXiv:1506.03481.
Li, W. and Fearnhead, P. (2016). Improved convergence of regression adjusted Approximate
Bayesian Computation. arXiv:1609.07135.
Marin, J.-M., Pillai, N. S., Robert, C. P. and Rousseau, J. (2014). Relevant statistics for Bayesian
model choice. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 76(5),
833–859.
Marin, J.-M., Raynal, L., Pudlo, P., Ribatet, M. and Robert, C. P. (2016). ABC random forests for
Bayesian parameter inference. arXiv.1605.05537.
Marjoram, P., Molitor, J., Plagnol, V. and Tavare, S. (2003). Markov chain Monte Carlo without
likelihoods. Proceedings of the National Academy of Sciences 100, 15324–15328.
Martin, G. M., McCabe, B. P., Maneesoonthorn, W. and Robert, C. P. (2016). Approximate
Bayesian computation in state space models. arXiv:1409.8363.

19

Nott, D. J., Fan, Y., Marshall, L. and Sisson, S. (2014). Approximate Bayesian computation and
Bayes’ linear analysis: toward high-dimensional ABC. Journal of Computational and Graphical
Statistics 23(1), 65–86.
Peters, G. W., Fan, Y. and Sisson, S. A. (2012). On sequential Monte Carlo, partial rejection control
and approximate Bayesian computation. Statistics and Computing 22(6), 1209–1222.
Prangle, D., Fearnhead, P., Cox, M. P., Biggs, P. J. and French, N. P. (2014). Semi-automatic
selection of summary statistics for ABC model choice. Statistical Applications in Genetics and
Molecular Biology 13(1), 67–82.
Pritchard, J. K., Seielstad, M. T., Perez-Lezaun, A. and Feldman, M. W. (1999). Population growth
of human Y chromosomes: a study of Y chromosome microsatellites. Molecular Biology and Evolution 16, 1791–1798.
Roberts, G. O., Rosenthal, J. S. et al. (2001). Optimal scaling for various Metropolis-Hastings
algorithms. Statistical Science 16(4), 351–367.
Ruli, E., Sartori, N. and Ventura, L. (2016). Approximate Bayesian computation with composite
score functions. Statistics and Computing 26(3), 679–692.
Sisson, S. A., Fan, Y. and Tanaka, M. M. (2007). Sequential Monte Carlo without likelihoods.
Proceedings of the National Academy of Sciences 104(6), 1760–1765.
Wegmann, D., Leuenberger, C. and Excoffier, L. (2009). Efficient approximate Bayesian computation coupled with Markov chain Monte Carlo without likelihood. Genetics 182(4), 1207–1218.
Wilkinson, R. D. (2013). Approximate Bayesian computation (ABC) gives exact results under the
assumption of model error. Statistical Applications in Genetics and Molecular Biology 12(2),
129–141.

20

