arXiv:1609.07135v2 [math.ST] 28 Nov 2017

Convergence of Regression-Adjusted Approximate Bayesian
Computation
B Y W ENTAO L I
School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne NE1
7RU, U.K.
wentao.li@newcastle.ac.uk
PAUL F EARNHEAD
Department of Mathematics and Statistics, Lancaster University, Lancaster LA1 4YF, U.K.
p.fearnhead@lancaster.ac.uk
AND

S UMMARY
We present asymptotic results for the regression-adjusted version of approximate Bayesian
computation introduced by Beaumont et al. (2002). We show that for an appropriate choice of
the bandwidth, regression adjustment will lead to a posterior that, asymptotically, correctly quantifies uncertainty. Furthermore, for such a choice of bandwidth we can implement an importance
sampling algorithm to sample from the posterior whose acceptance probability tends to unity
as the data sample size increases. This compares favourably to results for standard approximate
Bayesian computation, where the only way to obtain a posterior that correctly quantifies uncertainty is to choose a much smaller bandwidth; one for which the acceptance probability tends to
zero and hence for which Monte Carlo error will dominate.
Keywords: Approximate Bayesian computation; Importance sampling; Local-linear regression; Partial information.

1. I NTRODUCTION
Modern statistical applications increasingly require the fitting of complex statistical models,
which are often intractable in the sense that it is impossible to evaluate the likelihood function.
This excludes standard implementation of likelihood-based methods, such as maximum likelihood estimation or Bayesian analysis. To overcome this there has been substantial interest in
likelihood-free or simulation-based methods, which replace calculating the likelihood by simulation of pseudo datasets from the model. Inference can then be performed by comparing these
pseudo datasets, simulated for a range of different parameter values, to the actual data.
Examples of such likelihood-free methods include simulated methods of moments (Duffie &
Singleton, 1993), indirect inference (Gouriéroux & Ronchetti, 1993; Heggland & Frigessi, 2004),
synthetic likelihood (Wood, 2010) and approximate Bayesian computation (Beaumont et al.,
2002). Of these, approximate Bayesian computation methods are arguably the most common
methods for Bayesian inference, and have been popular in population genetics (e.g., Beaumont
et al., 2002; Cornuet et al., 2008), ecology (e.g., Beaumont, 2010) and systems biology (e.g.,
Toni et al., 2009); more recently they have seen increased use in other application areas, such as
econometrics (Calvet & Czellar, 2015) and epidemiology (Drovandi & Pettitt, 2011).
1

2
The idea of approximate Bayesian computation is to first summarize the data using lowdimensional summary statistics, such as sample means or autocovariances or suitable quantiles
of the data. The posterior density given the summary statistics is then approximated as follows.
Assume the data are Yobs = (yobs,1 , . . . , yobs,n ) and modelled as a draw from a parametric model
with parameter θ ∈ Rp . Let K(x) be a positive kernel, where maxx K(x) = 1, and ε > 0 is the
bandwidth. For a given d-dimensional summary statistic s(Y ), our model will define a density
fn (s | θ). We then define a joint density, πε (θ, s | sobs ), for (θ, s) as
´

π(θ)fn (s | θ)K{ε−1 (s − sobs )}
,
−1
Rp ×Rd π(θ)fn (s | θ)K{ε (s − sobs )} dθds

where sobs = s(Yobs ). Our approximation to the posterior density is the marginal density
ˆ
πε (θ | sobs ) = πε (θ, s | sobs ) ds,

(1)

(2)

which we call the approximate Bayesian computation posterior density. For brevity we will often
shorten this to posterior density in the following. We will always call the actual posterior given
the summary the true posterior.
The idea of approximate Bayesian computation is that we can sample from πε (θ | sobs ) without needing to evaluate the likelihood function or fn (s | θ). The simplest approach is via rejection sampling (Beaumont et al., 2002), which proceeds by simulating a parameter value and
an associated summary statistic from π(θ)fn (s | θ). This pair is then accepted with probability
K{ε−1 (s − sobs )}. The accepted pairs will be drawn from (1), and the accepted parameter values will be drawn from the posterior (2). Implementing this rejection sampler requires only the
ability to simulate pseudo data sets from the model, and then to be able to calculate the summary
statistics for those data sets.
Alternative algorithms for simulating from the posterior include adaptive or sequential importance sampling (Beaumont et al., 2009; Bonassi & West, 2015; Lenormand et al., 2013; Filippi
et al., 2013) and Markov chain Monte Carlo approaches (Marjoram et al., 2003; Wegmann et al.,
2009). These aim to propose parameter values in areas of high posterior probability, and thus can
be substantially more efficient than rejection sampling. However, the computational efficiency of
all these methods is limited by the probability of acceptance for data simulated with a parameter
value that has high posterior probability.
This paper is concerned with the asymptotic properties of approximate Bayesian computation. It builds upon Li & Fearnhead (2018) and Frazier et al. (2016), who present results on
the asymptotic behaviour of the posterior distribution and the posterior mean of approximate
Bayesian computation as the amount of data, n, increases. Their results highlight the tension
in approximate Bayesian computation between choices of the summary statistics and bandwidth
that will lead to more accurate inferences, against choices that will reduce the computational cost
or Monte Carlo error of algorithms for sampling from the posterior.
An informal summary of some of these results is as follows. Assume a fixed dimensional
summary statistic and that the true posterior variance given this summary decreases like 1/n
as n increases. The theoretical results compare the posterior, or posterior mean, of approximate
Bayesian computation, to the true posterior, or true posterior mean, given the summary of the
data. The accuracy of using approximate Bayesian computation is governed by the choice of
bandwidth, and this choice should depend on n. Li & Fearnhead (2018) shows that the optimal
choice of this bandwidth will be O(n−1/2 ). With this choice, estimates based on the posterior
mean of approximate Bayesian computation can, asymptotically, be as accurate as estimates
based on the true posterior mean given the summary. Furthermore the Monte Carlo error of

3
an importance sampling algorithm with a good proposal distribution will only inflate the mean
square error of the estimator by a constant factor of the form 1 + O(1/N ), where N is the number
of pseudo data sets. These results are similar to the asymptotic results for indirect inference,
where error for a Monte Carlo sample of size N also inflates the overall mean square error of
estimators by a factor 1 + O(1/N ) (Gouriéroux & Ronchetti, 1993). By comparison choosing
a bandwidth which is o(n−1/2 ) will lead to an acceptance probability that tends to zero as n →
∞, and the Monte Carlo error of approximate Bayesian computation will blow up. Choosing a
bandwidth that decays more slowly than O(n−1/2 ) will also lead to a regime where the Monte
Carlo error dominates, and can lead to a non-negligible bias in the posterior mean that inflates
the error.
While the above results for a bandwidth that is O(n−1/2 ) are positive in terms of point estimates, they are negative in terms of the calibration of the posterior. With such a bandwidth
the posterior density of approximate Bayesian computation always over-inflates the parameter
uncertainty: see Proposition 1 below and Theorem 2 of Frazier et al. (2016).
The aim of this paper is to show that a variant of approximate Bayesian computation can
yield inference that is both accurate in terms of point estimation, with its posterior mean having the same frequentist asymptotic variance as the true posterior mean given the summaries,
and calibrated, in the sense that its posterior variance equals this asymptotic variance, when
the bandwidth converges to zero at a rate slower than O(n−1/2 ). This means that the acceptance
probability of a good approximate Bayesian computation algorithm will tend to unity as n → ∞.

2. N OTATION AND S ET- UP
We denote the data by Yobs = (yobs,1 , . . . , yobs,n ), where n is the sample size, and each observation, yobs,i , can be of arbitrary dimension. Assume the data are modelled as a draw from
a parametric density, fn (y | θ), and consider asymptotics as n → ∞. This density depends on
an unknown parameter θ ∈ Rp . Let B p be the Borel sigma-field on Rp . We will let θ0 denote
the true parameter value, and π(θ) the prior distribution for the parameter. Denote the support of
π(θ) by P. Assume that a fixed-dimensional summary statistic sn (Y ) is chosen and its density
under our model is fn (s | θ). The shorthand Sn is used to denote the random variable with density fn (s | θ). Often we will simplify notation and write s and S for sn and Sn respectively. Let
N (x; µ, Σ) be the normal density at x with mean µ and variance Σ. Let Ac be the complement
of a set A with respect to the whole space. For a series xn we write xn = Θ(an ) if there exist
constants m and M such that 0 < m < |xn /an | < M < ∞ as n → ∞. For a real function g(x),
denote its gradient function at x = x0 by Dx g(x0 ). To simplify notation, Dθ is written as D.
Hereafter ε is considered to depend on n, so the notation εn is used.
The conditions of the theoretical results are stated below.
Condition 1. There exists some δ0 > 0, such that P0 = {θ : |θ − θ0 | < δ0 } ⊂ P, π(θ) ∈
C 2 (P0 ) and π(θ0 ) > 0.
´
´ Ql
Condition 2. The kernel satisfies (i) vK(v) dv = 0; (ii)
k=1 vik K(v) dv < ∞ for any
coordinates (vi1 , . . . , vil ) of v and l ≤ p + 6; (iii) K(v) ∝ K(kvk2Λ ) where kvk2Λ = v T Λv
and Λ is a positive-definite matrix, and K(v) is a decreasing function of kvkΛ ; (iv) K(v) =
α
O(e−c1 kvk 1 ) for some α1 > 0 and c1 > 0 as kvk → ∞.
Condition 3. There exists a sequence an , satisfying an → ∞ as n → ∞, a d-dimensional vector s(θ) and a d × d matrix A(θ), such that for all θ ∈ P0 ,
an {Sn − s(θ)} → N {0, A(θ)}, n → ∞,

4
in distribution. We also assume that sobs → s(θ0 ) in probability. Furthermore, (i) s(θ) and
A(θ) ∈ C 1 (P0 ), and A(θ) is positive definite for any θ; (ii) for any δ > 0 there exists
δ 0 > 0 such that ks(θ) − s(θ0 )k > δ 0 for all θ satisfying kθ − θ0 k > δ; and (iii) I(θ) =
Ds(θ)T A−1 (θ)Ds(θ) has full rank at θ = θ0 .
Let fen (s | θ) = N {s; s(θ), A(θ)/a2n } be the density of the normal approximation to S and
introduce the standardized random variable Wn (s) = an A(θ)−1/2 {S − s(θ)}. We further let
fWn (w | θ) and f˜Wn (w | θ) be the densities for Wn under the true model for S and under our
normal approximation to the model for S.
2/5

Condition 4. There exists αn satisfying αn /an → ∞ and a density rmax (w)
satisfying Condition 2 (ii)–(iii) where K(v) is replaced with rmax (w), such that
supθ∈P0 αn fWn (w | θ) − feWn (w | θ) ≤ c3 rmax (w) for some positive constant c3 .
Condition 5. The following statements hold: (i) rmax (w) satisfies Condition 2 (iv); and (ii)
α
supθ∈Pc0 fWn (w | θ) = O(e−c2 kwk 2 ) as kwk → ∞ for some positive constants c2 and α2 , and
A(θ) is bounded in P.
Conditions 1–5 are from Li & Fearnhead (2018). Condition 2 is a requirement for the kernel function and is satisfied by all commonly used kernels, such as any kernel with compact
support or the Gaussian kernel. Condition 3 assumes a central limit theorem for the summary
statistic with rate an , and, roughly speaking, requires the summary statistic to accumulate information when n. This is a natural assumption, since many common summary statistics are sample
moments, proportions, quantiles and autocorrelations, for which a central limit theorem would
apply. It is also possible to verify the asymptotic normality of auxiliary model-based or composite likelihood-based summary statistics (Drovandi et al., 2015; Ruli et al., 2016) by referring to
the rich literature on asymptotic properties of quasi maximum-likelihood estimators (Varin et al.,
2011) or quasi-posterior estimators (Chernozhukov & Hong, 2003). This assumption does not
cover ancillary statistics, using the full data as a summary statistic, or statistics based on distances, such as an asymptotically chi-square distributed test statistic. Condition 4 assumes that,
in a neighborhood of θ0 , fn (s | θ) deviates from the leading term of its Edgeworth expansion
−2/5
by a rate an . This is weaker than the standard requirement, o(a−1
n ), for the remainder from
Edgeworth expansion. It also assumes that the deviation is uniform, which is not difficult to satisfy in a compact neighborhood. Condition 5 further assumes that fn (s | θ) has exponentially
decreasing tails with rate uniform in the support of π(θ). This implies that posterior moments
from approximate Bayesian computation are dominated by integrals in the neighborhood of θ0
and have leading terms with concise expressions. With Condition 5 weakened, the requirement
of εn for the proper convergence to hold might depend on the specific tail behavior of fn (s | θ).
Additionally, for the results regarding regression adjustment the following moments of the
summary statistic are required to exist.
´
´
Condition 6. The first two moments, Rd sfn (s | θ) ds and Rd ssT fn (s | θ) ds, exist.

3.

A SYMPTOTICS OF APPROXIMATE BAYESIAN COMPUTATION
3·1. Posterior
First we consider the convergence of the posterior distribution of approximate Bayesian computation, denoted by Πε (θ ∈ A | sobs ) for A ∈ B p , as n → ∞. The distribution function is a
random function with the randomness due to sobs . We present two convergence results. One is

5
the convergence of the posterior distribution function of a properly scaled and centered version
of θ, see Proposition 1. The other is the convergence of the posterior mean, a result which comes
from Li & Fearnhead (2018) but, for convenience, is repeated as Proposition 2.
The following proposition gives three different limiting forms for Πε (θ ∈ A | sobs ), corresponding to different rates for how the bandwidth decreases relative to the rate of the central limit
theorem in Condition 3. We summarize these competing rates by defining cε = limn→∞ an εn .
P ROPOSITION 1. Assume Conditions 1–5. Let θε denote the posterior mean of approximate
−3/5
Bayesian computation. As n → ∞, if εn = o(an ) then the following convergence holds, depending on the value of cε .
(i) If cε = 0 then

ˆ
ψ(t) dt → 0,

sup Πε {an (θ − θε ) ∈ A | sobs } −

A∈B p

A

in probability, where ψ(t) = N {t; 0, I(θ0 )−1 }.
(ii) If cε ∈ (0, ∞) then for any A ∈ B p ,

ˆ

Πε {an (θ − θε ) ∈ A | sobs } →

ψ(t) dt,
A

in distribution, where
ˆ
ψ(t) ∝
N [t; cε β0 {v − EG (v)}, I(θ0 )−1 ]G(v) dv, β0 = I(θ0 )−1 Ds(θ0 )T A(θ0 )−1 ,
Rd

and G(v) is a random density of v, with mean EG (v), which depends on cε and Z ∼ N (0, Id ).
(iii) If cε = ∞ then
ˆ
−1
sup Πε {εn (θ − θε ) ∈ A | sobs } −
ψ(t) dt → 0,
A∈B p

A

in probability, where ψ(t) ∝ K{Ds(θ0 )t}.
For a similar result, under different assumptions, see Theorem 2 of Frazier et al. (2016). See also
Soubeyrand & Haon-Lasportes (2015) for related convergence results for the true posterior given
the summaries for some specific choices of summary statistics.
The explicit form of G(v) is stated in the Supplementary Material. When we have the same
number of summary statistics and parameters, d = p, the limiting distribution simplifies to
ˆ
ψ(t) ∝
N {Ds(θ0 )t; cε v, A(θ0 )}K(v) dv.
Rd

The more complicated form in Proposition 1 (ii) above arises from the need to project the summary statistics onto the parameter space. The limiting distribution may depend on the value of the
summary statistic, sobs , in the space orthogonal to Ds(θ0 )T A(θ0 )−1/2 . Hence the limit depends
on a random quantity, Z, which can be interpreted as the noise in sobs .
The main difference between the three convergence results is the form of the limiting density
ψ(t) for the scaled random variable an,ε (θ − θε ), where an,ε = an 1cε <∞ + ε−1
n 1cε =∞ . For case
(i) the bandwidth is sufficiently small that the approximation in approximate Bayesian computation due to accepting summaries close to the observed summary is asymptotically negligible. The
asymptotic posterior distribution is Gaussian, and equals the limit of the true posterior for θ given
the summary. For case (iii) the bandwidth is sufficiently big that this approximation dominates

6
and the asymptotic posterior distribution of approximate Bayesian computation is determined by
the kernel. For case (ii) the approximation is of the same order as the uncertainty in θ, which
leads to an asymptotic posterior distribution that is a convolution of a Gaussian distribution and
the kernel. Since the limit distributions of cases (i) and (iii) are non-random in the space L1 (Rp ),
the weak convergence is strengthened to convergence in probability in L1 (Rp ). See the proof in
Appendix A.
P ROPOSITION 2. (Theorem 3.1 of Li & Fearnhead, 2018) Assume conditions of Proposition
−3/5
−1
1. As n → ∞, if εn = o(an ), an (θε − θ0 ) → N {0, IABC
(θ0 )} in distribution. If εn = o(a−1
n )
or d = p or the covariance matrix of the kernel is proportional to A(θ0 ) then IABC (θ0 ) = I(θ0 ).
For other cases, I(θ0 ) − IABC (θ0 ) is semi-positive definite.
Proposition 2 helps us to compare the frequentist variability in the posterior mean of approximate Bayesian computation with the asymptotic posterior distribution given in Proposition 1.
If εn = o(a−1
n ) then the posterior distribution is asymptotically normal with variance matrix
−2
−1
an I(θ0 ) , and the posterior mean is also asymptotically normal with the same variance matrix. These results are identical to those we would get for the true posterior and posterior mean
given the summary.
For an εn which is the same order as a−1
n , the uncertainty in approximate Bayesian compu−1
tation has rate an . However the limiting posterior distribution, which is a convolution of the
true limiting posterior given the summary with the kernel, will overestimate the uncertainty by a
constant factor. If εn decreases slower than a−1
n , the posterior contracts at a rate εn , and thus will
over-estimate the actual uncertainty by a factor that diverges as n → 0.
In summary, it is much easier to get approximate Bayesian computation to accurately estimate
−3/5
the posterior mean. This is possible with εn as large as o(an ) if the dimension of the summary
statistic equals that of the parameter. However, accurately estimating the posterior variance, or
getting the posterior to accurately reflect the uncertainty in the parameter, is much harder. As
commented in Section 1, this is only possible for values of εn for which the acceptance probability in a standard algorithm will go to zero as n increases. In this case the Monte Carlo sample size,
and hence the computational cost, of approximate Bayesian computation will have to increase
substantially with n.
As one application of our theoretical results, consider observations that are independent and
identically distributed from a parametric density f (· | θ). One approach to construct the summary
statistics is to use the score vector of some tractable approximating auxiliary model evaluated at
the maximum auxiliary likelihood estimator (Drovandi et al., 2015). Ruli et al. (2016) constructs
an auxiliary
Q model from a composite likelihood, so the auxiliary likelihood for a single observation is i∈I f (y ∈ Ai | θ) where {Ai : i ∈ I } is a set of marginal or conditional events for y.
Denote the auxiliary score vector for a single observation by clθ (· | θ) and the maximum auxiliary likelihood estimator for our data set by θbcl . Then the summary statistic, s, for any pseudo
P
data set {y1 , . . . , yn } is nj=1 clθ (yj | θbcl )/n.
For y ∼ f (· | θ), assume the first two moments of clθ (y | θ0 ) exist and clθ (y | θ) is differentiable at θ. Let H(θ) = Eθ {∂clθ (y | θ0 )/∂θ} and J(θ) = varθ {clθ (y | θ0 )}. Then if θbcl is consistent for θ0 , Condition 3 is satisfied with
n1/2 [S − Eθ {clθ (Y | θ0 )}] → N {0, J(θ)}, n → ∞,
in distribution, and with I(θ0 ) = H(θ0 )T J(θ0 )−1 H(θ0 ).
Our results show that the posterior mean of approximate Bayesian computation, using εn =
O(n−1/2 ), will have asymptotic variance I(θ0 )−1 /n. This is identical to the asymptotic variance

7
of the maximum composite likelihood estimator (Varin et al., 2011). Furthermore, the posterior
variance will overestimate this just by a constant factor. As we show below, using the regression
correction of Beaumont et al. (2002) will correct this overestimation and produce a posterior that
correctly quantifies the uncertainty in our estimates.
An alternative approach to construct an approximate posterior using composite likelihood is
to use the product of the prior and the composite likelihood. In general, this leads to a poorly
calibrated posterior density which substantially underestimates uncertainty (Ribatet et al., 2012).
Adjustment of the composite likelihood is needed to obtain calibration, but this involves estimation of the curvature and the variance of the composite score (Pauli et al., 2011). Empirical
evidence that approximate Bayesian computation more accurately quantifies uncertainty than
alternative composite-based posteriors is given in Ruli et al. (2016).
3·2. Regression Adjusted Approximate Bayesian Computation
The regression adjustment of Beaumont et al. (2002) involves post-processing the output of
approximate Bayesian computation to try to improve the resulting approximation to the true posterior. Below we will denote a sample from the posterior of approximate Bayesian computation
by {(θi , si )}i=1,...,N . Under the regression adjustment, we obtain a new posterior sample by using {θi − βbε (si − sobs )}i=1,...,N where βbε is the least square estimate of the coefficient matrix
in the linear model
θi = α + β(si − sobs ) + ei ,

i = 1, . . . , N,

where ei are independent identically distributed errors.
We can view the adjusted sample as follows. Define a constant, αε , and a vector βε as
(αε , βε ) = arg min Eε [kθ − α − β(s − sobs )k2 | sobs ],
α,β

where expectation is with respect the joint posterior distribution of (θ, s) given by approximate
Bayesian computation. Then the ideal adjusted posterior is the distribution of θ∗ = θ − βε (s −
sobs ) where (θ, s) ∼ πε (θ, s). The density of θ∗ is
ˆ
πε∗ (θ∗ | sobs ) =
πε {θ∗ + βε (s − sobs ), s | sobs } ds
Rd

and the sample we get from regression-adjusted approximate Bayesian computation is a draw
from πε∗ (θ∗ | sobs ) but with βε replaced by its estimator.
The variance of πε∗ (θ∗ | sobs ) is strictly smaller than that of πε (θ | sobs ) provided s is correlated with θ. The following results, which are analogous to Propositions 1 and 2, show that this
reduction in variation is by the correct amount to make the resulting adjusted posterior correctly
quantify the posterior uncertainty.
T HEOREM 1. Assume Conditions 1–6. Denote the mean of πε∗ (θ∗ | sobs ) by θε∗ . As n → ∞, if
−3/5
εn = o(an ),
ˆ
∗
∗
sup Πε {an (θ − θε ) ∈ A | sobs } −
N {t; 0, I(θ0 )−1 } dt → 0,
A∈B p

A

in probability, and an (θε∗ − θ0 ) → N {0, I(θ0 )−1 } in distribution. Moreover, if βε is replaced by
βeε satisfying an εn (βeε − βε ) = op (1), the above results still hold.
The limit of the regression adjusted posterior distribution is the true posterior given the summary
−3/5
provided εn is o(an ). This is a slower rate than that at which the posterior contracts, which,

8
as we will show in the next section, has important consequences in terms of the computational
efficiency of approximate Bayesian computation. The regression adjustment corrects both the
additional noise of the posterior mean when d > p and the overestimated uncertainty of the
posterior. This correction comes from the removal of the first order bias caused by ε. Blum (2010)
shows that the regression adjustment reduces the bias of approximate Bayesian computation
when E(θ | s) is linear and the residuals θ − E(θ | s) are homoscedastic. Our results do not
require these assumptions, and suggest that the regression adjustment should be applied routinely
with approximate Bayesian computation provided the coefficients βε can be estimated accurately.
With the simulated sample, βε is estimated by βbε . The accuracy of βbε can be seen by the
following decomposition,
βbε = covN (s, θ)varN (s)−1




s − sε
s − sε −1
1
∗
∗
= βε +
covN
, an (θ − θε ) varN
,
an εn
εn
εn
where covN and varN are the sample covariance and variance matrices, and s is the sample
mean. Since cov(s, θ∗ ) = 0 and the distributions of s − sε and θ∗ − θε∗ contract at rates εn
−1 −1/2 }
b
and a−1
n respectively, the error βε − βε can be shown to have the rate Op {(an εn ) N
as n → ∞ and N → ∞. We omit the proof, since it is tedious and similar to the proof of the
asymptotic expansion of βε in Lemma 4. Thus, if N increases to infinity with n, βbε − βε will be
op {(an εn )−1 } and the convergence of Theorem 1 will hold instead.
Alternatively we can get an idea of the additional error for large N from the following proposition.
P ROPOSITION 3. Assume Conditions 1–6. Consider θ∗ = θ − βbε (s − sobs ). As n → ∞, if
−3/5
εn = o(an ) and N is large enough, for any A ∈ B p ,
ˆ
∗
∗
Πε {an (θ − θε ) ∈ A | sobs } →
ψ(t) dt,
A

in distribution, where

ˆ

h
N t;

ψ(t) ∝
Rd

when cε < ∞,

ˆ
ψ(t) ∝

n
N t;

Rp

i
η
−1
{v
−
E
(v)},
I(θ
)
G(v) dv,
0
G
N 1/2

o
η
0
−1
Ds(θ
)t
,
I(θ
)
K{Ds(θ0 )t0 } dt0 ,
0
0
N 1/2

when cε = ∞, and η = Op (1).
The limiting distribution here can be viewed as the convolution of the limiting distribution obtained when the optimal coefficients are used and that of a random variable, which relates to the
error in our estimate of βε , and that is Op (N −1/2 ).
3·3. Acceptance Rates when ε is Negligible
Finally we present results for the acceptance probability of approximate Bayesian computation, the quantity that is central to the computational cost of importance sampling or Markov
chain Monte Carlo-based algorithms. We consider a set-up where we propose the parameter
value from a location-scale family. That is, we can write the proposal density as the density of a
random variable, σn X + µn , where X ∼ q(·), E(X) = 0 and σn and µn are constants that can

9
depend on n. The average acceptance probability pacc,q would then be
ˆ
qn (θ)fn (s | θ)K{ε−1
n (s − sobs )} dsdθ,
P×Rd

where qn (θ) is the density of σn X + µn . This covers the proposal distribution in fundamental
sampling algorithms, including the random-walk Metropolis algorithm and importance sampling
with unimodal proposal distribution, and serves as the building block for many advanced algorithms where the proposal distribution is a mixture of distributions from location-scale families,
such as iterative importance sampling-type algorithms.
We further assume that σn (µn − θ0 ) = Op (1), which means θ0 is in the coverage of qn (θ).
This is a natural requirement for any good proposal distribution. The prior distribution and θ0 as a
point mass are included in this proposal family. This condition would also apply to many Markov
chain Monte Carlo implementations of approximate Bayesian computation after convergence.
−1
As above, define an,ε = an 1cε <∞ + ε−1
n 1cε =∞ to be the smaller of an and εn . Asymptotic
−1
results for pacc,q when σn has the same rate as an,ε are given in Li & Fearnhead (2018). Here we
extend those results to other regimes.
−1/2

T HEOREM 2. Assume the conditions of Proposition 1. As n → ∞, if εn = o(an ): (i) if
−1
cε = 0 or σn /a−1
n,ε → ∞, then pacc,q → 0 in probability; (ii) if cε ∈ (0, ∞) and σn /an,ε →
−1
r1 ∈ [0, ∞), or cε = ∞ and σn /an,ε → r1 ∈ (0, ∞), then pacc,q = Θp (1); (iii) if cε = ∞ and
σn /a−1
n,ε → 0, then pacc,q → 1 in probability.
The proof of Theorem 2 can be found in the Supplementary Material. The underlying intuition
is as follows. For the summary statistic, s, sampled with parameter value θ, the acceptance probability depends on
s − sobs
1
= [{s − s(θ)} + {s(θ) − s(θ0 )} + {s(θ0 ) − sobs }],
εn
εn

(3)

where s(θ) is the limit of s in Condition 3. The distance between s and sobs is at least Op (a−1
n ),
−1
−1
since the first and third bracketed terms are Op (an ). If εn = o(an ) then, regardless of the value
of θ, (3) will blow up as n → ∞ and hence pacc,q goes to 0. If εn decreases with a rate slower
than a−1
n , (3) will go to zero providing we have a proposal which ensures that the middle term is
op (εn ), and hence pacc,q goes to unity.
Theorem 1 shows that, without the regression adjustment, approximate Bayesian computation
requires εn to be o(a−1
n ) if its posterior is to converge to the true posterior given the summary. In
this case Theorem 2 shows that the acceptance rate will degenerate to zero as n → ∞ regardless
of the choice of q(·). On the other hand, with the regression adjustment, we can choose εn =
−3/5
o(an ) and still have convergence to the true posterior given the summary. For such a choice,
if our proposal density satisfies σn = o(εn ), the acceptance rate will go to unity as n → ∞.

4. N UMERICAL E XAMPLE
Here we illustrate the gain of computational efficiency from using the regression adjustment on
the g-and-k distribution, a popular model for testing approximate Bayesian computation methods
(e.g., Fearnhead & Prangle, 2012; Marin et al., 2014). The data are independent and identically
distributed from a distribution defined by its quantile function,


1 − exp{−γz(x)}
−1
F (x; α, β, γ, κ) =α + β 1 + 0.8
{1 + z(x)2 }κ z(x), x ∈ [0, 1],
1 + exp{−γz(x)}

10
where α and β are location and scale parameters, γ and κ are related to the skewness and kurtosis
of the distribution, and z(x) is the corresponding quantile of a standard normal distribution. No
closed form is available for the density but simulating from the model is straightforward by
transforming realisations from the standard normal distribution.
In the following we assume the parameter vector (α, β, γ, κ) has a uniform prior in [0, 10]4
and multiple datasets are generated from the model with (α, β, γ, κ) = (3, 1, 2, 0.5). To illustrate
the asymptotic behaviour of approximate Bayesian computation, 50 data sets are generated for
each of a set of values of n ranging from 500 to 10, 000. Consider estimating the posterior
means, denoted by µ = (µ1 , . . . , µ4 ), and standard deviations, denoted by σ = (σ1 , . . . , σ4 ), of
the parameters. The summary statistic is a set of evenly spaced quantiles of dimension 19.
The bandwidth is chosen via fixing the proportion of the Monte Carlo sample to be accepted,
and the accepted proportions needed to achieve certain approximation accuracy for estimates
with and without the adjustment are compared. A higher proportion means more simulated parameter values can be kept for inference. The accuracy is measured by the average relative errors
of estimating µ or σ,
4

REµ =

µk − µ k |
1 X |b
,
4
µk
k=1

4

REσ =

1 X |b
σk − σk |
,
4
σk
k=1

for estimators µ
b = (b
µ1 , . . . , µ
b4 ) and σ
b = (b
σ1 , . . . , σ
b4 ). The proposal distribution is normal, with
the covariance matrix selected to inflate the posterior covariance matrix by a constant factor c2
and the mean vector selected to differ from the posterior mean by half of the posterior standard
deviation, which avoids the case that the posterior mean can be estimated trivially. We consider a
series of increasing c in order to investigate the impact of the proposal distribution getting worse.
Figure 1 shows that the required acceptance rate for the regression adjusted estimates is higher
than that for the unadjusted estimates in almost all cases. For estimating the posterior mean, the
improvement is small. For estimating the posterior standard deviations, the improvement is much
larger. To achieve each level of accuracy, the acceptance rates of the unadjusted estimates are
all close to zero. Those of the regression-adjusted estimates are higher by up to two orders of
magnitude, so the Monte Carlo sample size needed to achieve the same accuracy can be reduced
correspondingly.

5. D ISCUSSION
One way to implement approximate Bayesian computation so that the acceptance probability
tends to unity as n increases is to use importance sampling with a suitable proposal from a
location-scale family. The key difficulty with finding a suitable proposal is to ensure that the
location parameter is close to the true parameter, where close means the distance is O(εn ). This
can be achieved by having a preliminary analysis of the data, and using the point estimate of the
parameter from this preliminary analysis as the location parameter (Beaumont et al., 2009; Li &
Fearnhead, 2018).

ACKNOWLEDGMENT
This work was funded by the Engineering and Physical Sciences Research Council, under the
i-like programme grant.

11

0.8
0.6
0.4
0.0

0.2

Acceptance rate

0.8
0.6
0.4
0.2
0.0

Acceptance rate

1.0

(b)

1.0

(a)

2.0

2.5

3.0

3.5

4.0

4.5

2.0

2.5

3.0

c

4.0

4.5

3.5

4.0

4.5

c

0.8
0.6
0.4
0.2
0.0

0.0

0.2

0.4

0.6

Acceptance rate

0.8

1.0

(d)

1.0

(c)

Acceptance rate

3.5

2.0

2.5

3.0

3.5
c

4.0

4.5

2.0

2.5

3.0
c

Fig. 1: Acceptance rates required for different degrees of accuracy of approximate Bayesian
computation and different variances of the proposal distribution (which are proportional to c). In
each plot we show results for standard (grey-line) and regression adjusted (black-line) approximate Bayesian computation and for different values of n: n = 500 (dotted), n = 3, 000 (dashed)
and n = 10, 000 (solid). The averages over 50 data sets (thick) and their 95% confidence intervals (thin) are reported. Results are for a relative error of 0.08 and 0.05 in the posterior mean, in
(a) and (b) respectively, and for a relative error of 0.2 and 0.1 in the posterior standard deviation,
in (c) and (d) respectively.

S UPPLEMENTARY M ATERIAL
Proofs of lemmas and Theorem 2 are included in the online supplementary material.

12

A PPENDIX
Proof of Result from Section 3·1
Throughout the data are considered to be random. For any integer l > 0 and a set A ⊂ Rl ,
we use the convention that cA + x denotes the set {ct + x : t ∈ A} for c ∈ R and x ´∈ Rl . For a
non-negative function h(x), integrable in Rl , denote the normalised function h(x)/ Rl h(x) dx
by h(x)(norm) . For a vector x, denote a general polynomial of elements of x with degree up to l by
Pl (x). For any fixed δ < δ0 , let Bδ be the neighborhood {θ : kθ − θ0 k < δ}. Let πδ (θ) be π(θ)
truncated in Bδ . Recall that fen (s | θ) denotes the normal density with mean s(θ) and covariance
matrix A(θ)/a2n . Let t(θ) = an,ε (θ − θ0 ) and v(s) = ε−1
n (s − sobs ), rescaled versions θ and s.
For any A ∈ B p , let t(A) be the set {φ : φ = t(θ) for some θ ∈ A}.
e ε (θ ∈ A | sobs ) to be the normal counterpart of Πε (θ ∈ A | sobs ) with truncated prior,
Define Π
obtained by replacing π(θ) and fn (s | θ) in Πε by πδ (θ) and fen (s | θ). So let
π
eε (θ, s | sobs ) = ´

´

Bδ

πδ (θ)fen (s | θ)K{ε−1
n (s − sobs )}
,
−1
e
d πδ (θ)fn (s | θ)K{εn (s − sobs )} dθds

R

´

e ε (θ ∈ A | sobs ) be the distribution function with denπ
eε (θ | sobs ) = Rd π
eε (θ, s | sobs ) ds and Π
e ε by θeε . Let Wobs = an A(θ0 )−1/2 {sobs − s(θ0 )} and
sity π
eε (θ | sobs ). Denote the mean of Π
−1
T
−1
β0 = I(θ0 ) Ds(θ0 ) A(θ0 ) . By Condition 3, Wobs → Z in distribution as n → ∞, where
Z ∼ N (0, Id ).
e ε is an incorrect model for
Since the approximate Bayesian computation likelihood within Π
sobs , standard posterior convergence results do not apply. However, if we condition on the value
of the summary, s, then the distribution of θ is just the true posterior given s. Thus we can express
the posterior from approximate Bayesian computation as a continuous mixture of these true
−1
−1 e
posteriors. Let π
eε,tv (t, v) = a−d
n,ε πδ (θ0 + an,ε t)fn (sobs + εn v | θ0 + an,ε t)K(v). For any A ∈
e ε as,
B p , we rewrite Π
ˆ ˆ
e ε (θ ∈ A | sobs ) =
e ∈ A | sobs + εn v)e
Π
Π(θ
πε,tv (t, v)(norm) dtdv,
(4)
Rd

t(Bδ )

e ∈ A | s) is the posterior distribution with prior πδ (θ) and likelihood fen (s | θ).
where Π(θ
e ∈ A | sobs + εn v)
Using results from Kleijn & van der Vaart (2012), the leading term of Π(θ
can be obtained and is stated in the following lemma.
d
L EMMA 1. Assume Conditions 3 and 4. If εn = O(a−1
n ), for any fixed v ∈ R and small
enough δ,
ˆ
e
N [t; β0 {A(θ0 )1/2 Wobs + cε v}, I(θ0 )−1 ] dt → 0,
sup Π{an (θ − θ0 ) ∈ A | sobs + εn v} −
A∈B p

A

in probability as n → ∞.
This leading term is Gaussian, but with a mean that depends on v. Thus asymptotically, the
posterior of approximate Bayesian computation is the distribution
´ of the sum of a Gaussian
random variable and β0 cε V , where V has density proportional to π
eε,tv (t, v) dt.
To make this argument rigorous, and to find the distribution of this sum of random variables we
need to introduce several functions that relate to the limit of π
eε,tv (t0 , v). For a rank-p d × p matrix
A, a rank-d d × d matrix B and a d-dimensional vector c, define g(v; A, B, c) = exp[−(c +

13
Bv)T {I − A(AT A)−1 AT }(c + Bv)/2]/(2π)(d−p)/2 . Let
( 
N nDs(θ0 )t; an εn v + A(θ0 )1/2 Wobs , A(θ0 ) K(v),
o
gn (t, v) =
N Ds(θ0 )t; v + an1εn A(θ0 )1/2 Wobs , a21ε2 A(θ0 ) K(v),
n n

cε < ∞,
cε = ∞,

Gn (v) be g{v; A(θ0 )−1/2 Ds(θ0 ), an εn A(θ0 )−1/2 , Wobs }K(v), and EGn (·) be´ the expectation
under the density Gn (v)(norm) . In both cases it is straightforward to show that Rp gn (t, v) dt =
|A(θ0 )|−1/2 Gn (v). Additionally, for the case cε = ∞, define v 0 (v, t) = A(θ0 )1/2 Wobs +
an εn v − an εn Ds(θ0 )t and


1 0
1
gn0 (t, v 0 ) = N {v 0 ; 0, A(θ0 )}K Ds(θ0 )t +
v −
A(θ0 )1/2 Wobs .
an εn
an εn
Then with the transformation v 0 = v 0 (v, t), gn (t, v)dv = gn0 (t, v 0 )dv 0 . Let
(
N {Ds(θ0 )t; cε v + A(θ0 )1/2 Z, A(θ0 )}K(v), cε < ∞,
g(t, v) =
K{Ds(θ0 )t}N {v; 0, A(θ0 )},
cε = ∞,
−1/2 , Z}K(v) and E (·) be the expectation under the
G(v) be g{v; A(θ0 )−1/2 Ds(θ0 ), cε A(θ
G
´ 0)
(norm)
−1/2 G(v).
density G(v)
.
When
c
<
∞,
g(t,
v)
dt
=
|A(θ
)|
ε
0
Rp
´ ´
Expansions of Rd t(Bδ ) π
eε,tv (t, v)dtdv are given in the following lemma.
−1/2

Conditions
1–3.
If
εn´ = ´o(an ),
then
´ L´EMMA 2. Assume
πε,tv (t, v) − π(θ0 )gn (t, v)| dtdv → 0 in probability and Rd t(Bδ ) gn (t, v) dtdv =
Rd t(Bδ ) |e
´ ´
Θp (1), as n → ∞. Furthermore, for l ≤ 6, Rd t(Bδ ) Pl (v)gn (t, v) dtdv converges
´
−1/2
to
in distribution when cε < ∞ and converges to
Rd Pl (v)G(v) dv
´ |A(θ0 )|
P
{Ds(θ
)t}K{Ds(θ
)t}
dt
in
probability
when cε = ∞, as n → ∞.
p
0
0
l
R
e ε are asymptotically the same and gives an expansion
The following lemma states that Πε and Π
e
of θε .
−1/2

L EMMA 3. Assume Conditions 1–4. If εn = o(an ), then
e ε (θ ∈ B c | sobs ) are op (1);
(a) for any δ < δ0 , Πε (θ ∈ Bδc | sobs ) and Π
δ
e ε (θ ∈ A ∩ Bδ | sobs ) =
(b) there exists a δ < δ0 such that supA∈Bp Πε (θ ∈ A ∩ Bδ | sobs ) − Π
op (1);
(c) if, in addition, Condition 5 holds, then an,ε (θε − θeε ) = op (1), and θeε = θ0 +
1/2 W
−1
a−1
obs + εn β0 EGn (v) + rn,1 where the remainder rn,1 = op (an ).
n β0 A(θ0 )

Proof of Proposition 1. Lemma 10 in the Supplementary Material shows that Πε {an,ε (θ −
e ε {an,ε (θ − θeε ) ∈ A | sobs } have the same limit, in distribution when cε ∈
θε ) ∈ A | sobs } and Π
(0, ∞) and in total variation form when cε = 0 or ∞. Therefore it is sufficient to only consider
e ε of the properly scaled and centered θ.
the convergence of Π
e ε {an (θ − θeε ) ∈ A | sobs } equals
When an εn → cε < ∞, according to (4), Π
ˆ ˆ
e n (θ − θ0 ) ∈ A + an (θeε − θ0 ) | sobs + εn v}e
Π{a
πε,tv (t0 , v)(norm) dt0 dv.
Rd

t(Bδ )

14
By Lemma 1 and Lemma 3(c), we have

ˆ
N {t; µn (v), I(θ0 )−1 } dt = op (1),

e n (θ − θ0 ) ∈ A + an (θeε − θ0 ) | sobs + εn v} −
sup Π{a

A∈B p

A

where µn (v) = β0 {cε v − an εn EGn (v)} − an rn,1 . Then with Lemma 2, the leading term of
e ε {an (θ − θeε ) ∈ A | sobs } equals
Π
ˆ ˆ
ˆ
e
e
sup Πε {an (θ − θε ) ∈ A | sobs } −
N {t; µn (v), I(θ0 )−1 }gn (t0 , v)(norm) dtdt0 dv = op (1).
A∈B p

Rd

t(Bδ )

A

(5)
The numerator of the leading term of (5) is in the form
ˆ ˆ
ˆ
N {t; cε β0 v + x3 , I(θ0 )−1 }N {Ds(θ0 )t0 ; x1 v + x2 , A(θ0 )}K(v) dtdt0 dv,
Rd

t(Bδ )

A

where x1 ∈ R, x2 ∈ Rd and x3 ∈ Rp . This is continuous by Lemma 9 in the Supplementary
Material. Then since EGn (v) → EG (v) in distribution as n → ∞ by Lemma 2, we have
ˆ ˆ
ˆ
N {t; µn (v), I(θ0 )−1 }gn (t0 , v) dtdt0 dv
Rd t(Bδ ) A
ˆ ˆ ˆ
→
N [t; cε β0 {v − EG (v)}, I(θ0 )−1 ]g(t0 , v) dtdt0 dv,
Rd

Rp

A

in distribution as n → ∞. Putting the above results together, it holds that
ˆ ˆ
e ε {an (θ − θeε ) ∈ A | sobs } →
Π
N [t; cε β0 {v − EG (v)}, I(θ0 )−1 ]G(v)(norm) dvdt,
A

Rp

in distribution, and statement (ii) of the proposition holds.
When cε = 0, since µn (v) does not depend on v, (5) becomes
ˆ
e
e
N {t; −an εn β0 EGn (v) − an rn,1 , I(θ0 )−1 } dt = op (1),
sup Πε {an (θ − θε ) ∈ A | sobs } −

A∈B p

A

and by the continuous mapping theorem (van der Vaart, 2000),
ˆ
N {t; −an εn β0 EGn (v) − an rn,1 , I(θ0 )−1 } − N {t; 0, I(θ0 )−1 } dt = op (1).
Rp

´
e ε {an (θ − θeε ) ∈ A | sobs } − N {t; 0, I(θ0 )−1 } dt = op (1), and stateTherefore supA∈Bp Π
A
ment (i) of the proposition holds.
When cε = ∞, Lemma 1 cannot be applied to the posterior distribution within (4) directly.
e ε {θ ∈ A | sobs } equals,
With transformation v 00 = an εn v, Π
ˆ ˆ
−1 00 (norm) 0 00
e ∈ A | sobs + a−1 v 00 )e
Π(θ
πε,tv (t0 , a−1
dt dv ,
n
n εn v )
Rd

t(Bδ )

e ε {ε−1 (θ − θeε ) ∈ A | sobs } equals
which implies that Π
n
ˆ ˆ
−1 00 (norm)
e n (θ − θ0 ) ∈ an εn A + an (θeε − θ0 ) | sobs + a−1 v 00 }e
Π{a
πε,tv (t0 , a−1
dt0 dv 00 .
n
n εn v )
Rd

t(Bδ )

15
Then Lemma 1 can be applied. Using Lemma 2 and transforming v 00 back to v we have
e ε {ε−1 (θ − θeε ) ∈ A | sobs }
sup Π
n

A∈B p

ˆ

ˆ

ˆ
(an εn )p N {an εn t; µ0n (v), I(θ0 )−1 }gn (t0 , v)(norm) dtdt0 dv = op (1),

−
Rd

t(Bδ )

(6)

A

where µ0n (v) = β0 {A(θ0 )1/2 Wobs + an εn v} − an (θeε − θ0 ).
Let t00 (t, t0 ) = an εn (t − t0 ) + an (θeε − θ0 ) and t00 (A, Bδ ) be the set {t00 (t, t0 ) : t ∈ A, t0 ∈
t(Bδ )}. With transformations v 0 = v 0 (v, t0 ) and t00 = t00 (t, t0 ), since β0 Ds(θ0 ) = Ip , we have
ˆ ˆ
ˆ
(an εn )p N {an εn t; µ0n (v), I(θ0 )−1 }gn (t0 , v) dtdt0 dv
d
R
t(B ) A
ˆ ˆδ
ˆ
(an εn )p N {an εn (t − t0 ); β0 v 0 − an (θeε − θ0 ), I(θ0 )−1 }gn0 (t0 , v 0 ) dtdt0 dv 0
=
Rd t(Bδ ) A


ˆ
ˆ ˆ
1
00
0
−1 0
00
0
e
e
=
N {t ; β0 v − an (θε − θ0 ), I(θ0 ) }gn t −
{t − an (θε − θ0 )}, v dtdt00 dv 0 .
a n εn
Rd t00 (A,Bδ ) A
The idea now is that as n → ∞, an εn → ∞, so the gn0 term in the integral will tend to gn0 (t, v 0 ).
Then by integrating first with respect to t00 and then with respect to v, we get the required result.
To make this argument rigorous, consider the following function,
ˆ ˆ
ˆ
N {t00 ; β0 v 0 , I(θ0 )−1 }N {v 0 ; 0, A(θ0 )}
Rd

t00 (A,Bδ )

A

× K{Ds(θ0 )t + x1 v 0 − x2 t00 + x3 } − K{Ds(θ0 )t} dtdt00 dv 0 ,



where x1 ∈ R, x2 ∈ R and x3 ∈ Rd . This is continuous by Lemma 9 in the Supplementary Material, so by the continuous mapping theorem,


ˆ
ˆ ˆ
ˆ
1 00 0
00
0
−1 0
00 0
N {t ; β0 v , I(θ0 ) }gn t −
sup
t , v dtdt dv −
K{Ds(θ0 )t} dt = op (1).
an εn
A∈B p
A
Rd t00 (A,Bδ ) A
Then using Lemma 2,

ˆ

e ε {ε−1 (θ − θeε ) ∈ A | sobs } −
sup Π
n

A∈B p

ˆ
K{Ds(θ0 )t} dt/

A

K{Ds(θ0 )t} dt = op (1).
Rp

Therefore statement (iii) of the proposition holds.
Proof of Result from Section 3·2
The intuition behind this result is that, as shown above, the joint posterior of (θ, s) under
approximate Bayesian computation can be viewed as a marginal distribution for s times a conditional for θ given s. The latter is just the true posterior for θ given s, and this posterior converges
to a Gaussian limit that depends on s only through its mean. Regression adjustment works because it corrects for the dependence of this mean on s, so if we work with the regression adjusted
parameter θ∗ then the conditional distribution for θ∗ given s will tend to a Gaussian limit whose
mean is the same for all s.
The following lemma gives an expansion of the optimal linear coefficient matrix βε . The order
of the remainder leads to successful removal of the first order bias in the posterior distribution of
approximate Bayesian computation.

16
−3/5

L EMMA 4. Assume Conditions 1–6. Then if εn = o(an

), an εn (βε − β0 ) = op (1).

The following lemma, similar to Lemma 3, says that the approximate Bayesian computation
e ε . Recall that θ∗ is the mean of π ∗ (θ∗ |
posterior distribution of θ∗ is asymptotically the same as Π
ε
ε
´
∗
∗
∗
sobs ). Let π
eε (θ | sobs ) = Rd π
eε {θ + βε (s − sobs ), s | sobs } ds and θeε∗ be the mean of π
eε∗ (θ∗ |
sobs ).
−3/5

L EMMA 5. Assume Conditions 1–6. If εn = o(an ), then
e ε (θ∗ ∈ B c | sobs ) are op (1);
(a) for any δ < δ0 , Πε (θ∗ ∈ Bδc | sobs ) and Π
δ
∗
e ε (θ∗ ∈ A ∩ Bδ | sobs ) =
(b) there exists δ < δ0 such that supA∈Bp Πε (θ ∈ A ∩ Bδ | sobs ) − Π
op (1);
1/2 W
(c) an (θε∗ − θeε∗ ) = op (1), and θeε∗ = θ0 + a−1
obs + εn (β0 − βε )EGn (v) + rn,2
n β0 A(θ0 )
−1
where the remainder rn,2 = op (an ).
Proof of Theorem 1. Similar to the proof of Proposition 1, it is sufficient to only consider the
e ε of the properly scaled and centered θ∗ . Similar to (4),
convergence of Π
(´ ´
e ∈ A + εn βε v | sobs + εn v)e
Π(θ
πε,tv (t0 , v)(norm) dt0 dv,
e ε (θ∗ ∈ A | sobs ) = ´Rd ´t(Bδ )
Π
−1 π
0 −1 −1 (norm) dt0 dv,
e
ε,tv (t , an εn v)
Rd t(Bδ ) Π(θ ∈ A + εn βε v | sobs + an v)e

cε < ∞,
cε = ∞.

Similar to (5), by Lemma 1 and Lemma 5(c), we have
ˆ
e ε {an (θ∗ − θeε∗ ) ∈ A | sobs } −
sup Π

A∈B p

Rd

ˆ

ˆ
N {t; µ∗n (v), I(θ0 )−1 }gn (t0 , v)(norm) dtdt0 dv = op (1),

t(Bδ )

A

(7)
where
µ∗n (v) = {an εn (β0 − βε ) + (cε − an εn )β0 1cε <∞ }{v − EGn (v)} + an rn,2 .
Since
an εn (βε − β0 ) = op (1), by Lemma 9 in the Supplementary Material and the continuous
mapping theorem,
ˆ ˆ
ˆ
N {t; µ∗n (v), I(θ0 )−1 } − N {t; 0, I(θ0 )−1 } gn (t0 , v) dtdt0 dv = op (1).
Rd

t(Bδ )

Rp

Then we have

ˆ

e ε {an (θ∗ − θe∗ ) ∈ A | sobs } −
sup Π
ε

A∈B p

N {t; 0, I(θ0 )−1 } dt = op (1),
A

and the first convergence in the theorem holds.
The second convergence in the theorem holds by Lemma 5(c).
0
Since the only requirement for βε in the above is an εn (βε − βε ) = op (1), the above arguments
will hold if βε is replaced by a p × d matrix βbε satisfying an εn (βbε − βε ) = op (1).

Proof of Proposition 3. Consider θ∗ where βε is replaced by βbε . Let ηn = N 1/2 an εn (βbε −
βε ). Since βbε − βε = Op {(an εn )−1 N −1/2 } as n → ∞, ηn = Op (1) as n → ∞ and let its limit
be η. In this case, if we replace µ∗n (v) in (7) with µ∗n (v) + N −1/2 ηn {v − EGn (v)}, denoted by
µ
b∗n (v), the equation still holds. Denote this equation by (70 ). Limits of the leading term in (70 )
can be obtained by arguments similar as those for (5).

17
When cε < ∞, since for fixed v, µ
b∗n (v) converges to N −1/2 η{v − EG (v)} in distribution, by
following the same line we have
ˆ ˆ
∗
∗
e
e
N [t; N −1/2 η{v − EG (v)}, I(θ0 )−1 ]G(v)(norm) dvdt,
Πε {an (θ − θε ) ∈ A | sobs } →
A

Rp

in distribution as n → ∞.
cε = ∞, by ´Lemma 2 we have EGn (v) → 0 in probability, and
´ When
´
0
0
g
Rd t(Bδ ) n (t , v) dt dv → Rp K{Ds(θ0 )t} dt in probability. Then with transformation
v 0 = v 0 (v, t0 ), for fixed v,


1
1 0
1/2
∗
−1/2
0
v −
A(θ0 ) Wobs − EGn (v)
µ
bn (v) = {an εn (β0 − βε ) + N
ηn } Ds(θ0 )t +
an εn
an εn
→ N −1/2 ηDs(θ0 )t0 ,
in distribution as n → ∞. Recall that gn (t0 , v) dv = gn0 (t0 , v 0 ) dv 0 . Then by Lemma 9 in the Supplementary Material and the continuous mapping theorem,
ˆ ˆ
ˆ
N {t; µ
b∗n (v), I(θ0 )−1 }gn0 (t0 , v 0 ) dtdt0 dv
d
R
t(B ) A
ˆ ˆ δ ˆ
N {t; N −1/2 ηDs(θ0 )t0 , I(θ0 )−1 }K{Ds(θ0 )t0 } dtdt0 dv,
→
Rd

t(Bδ )

A

in distribution as n → ∞. Therefore by (70 ) and the above convergence results,
ˆ ˆ
e ε {an (θ∗ − θeε∗ ) ∈ A | sobs } →
Π
N {t; N −1/2 ηDs(θ0 )t0 , I(θ0 )−1 }K{Ds(θ0 )t0 }(norm) dt0 dt,
A

in distribution as n → ∞.

Rp



S UPPLEMENTARY M ATERIAL
Notations and Set-up
First
and conventions´ are given.
For two sets A and B, the sum of inte´ some limit notations
´
´
grals A f (x) dx + B f (x) dx is written as ( A + B )f (x) dx. For a constant d × p matrix A, let
the minimum and maximum eigenvalues of AT A be λ2min (A) and λ2max (A) where λmin (A) and
λmax (A) are non-negative. Obviously, for any p-dimension vector x, λmin (A)kxk ≤ kAxk ≤
λmax (A)kxk. For two matrices A and B, we say A is bounded by B or A ≤ B if λmax (A) ≤
λmin (B). For a set of matrices {Ai : i ∈ I} for some index set I, we say it is bounded if
λmax (Ai ) are uniformly bounded in i. Denote the identity matrix with dimension d by Id . Notations from the main text will also be used.
The following basic asymptotic results (Serfling, 2009) will be used throughout.
L EMMA 6. (i) For a series of random variables Zn , if Zn → Z in distribution as n → ∞,
Zn = Op (1). (ii) (Continuous mapping) For a series of continuous function gn (x), if gn (x) =
O(1) almost everywhere, then gn (Zn ) = Op (1), and this also holds if O(1) and Op (1) are replaced by Θ(1) and Θp (1).
Some notations regarding the posterior distribution of approximate Bayesian computation are
given. For A ⊂ Rp and a scalar function h(θ, s), let
ˆ ˆ
−d
πA (h) =
h(θ, s)π(θ)fn (s | θ)K{ε−1
n (s − sobs )}εn dsdθ,
A

Rd

18
and

ˆ ˆ
π
eA (h) =

A

Rd

−d
h(θ, s)πδ (θ)fen (s | θ)K{ε−1
n (s − sobs )}εn dsdθ.

e ε (θ ∈ A | sobs ) =
Then Πε (θ ∈ A | sobs ) = πA (1)/πP (1) and its normal counterpart Π
π
eA (1)/e
πP (1).
The following results from Li & Fearnhead (2018) will be used throughout.
L EMMA 7. Assume Conditions 1–4. Then as n → ∞,
(i) if Condition 5 also holds then, for any δ < δ0 , πBδc (1) and π
eBδc (1) are op (1), and
αδ

Op (e−an,ε cδ ) for some positive constants cδ and αδ depending on δ;
eA (1)| /e
πBδ (1) = Op (αn−1 );
eBδ (1){1 + Op (αn−1 )} and supA⊂Bδ |πA (1) − π
(ii) πBδ (1) = π
−1/2

eP (1) and πP (1) are
(iii) if εn = o(an ), π
eBδ (1) and πBδ (1) are Θp (ad−p
n,ε ), and thus π
d−p
Θp (an,ε );
−1/2
−3/5
(iv) if εn = o(an ) and Condition 5 holds, θε = θeε + op (a−1
), θε = θeε +
n,ε ). If εn = o(an
op (a−1
n ).
Proof. (i) is from Li & Fearnhead (2018, Lemma 3) and a trivial modification of its proof
when Condition 5 does no hold; (ii) is from Li & Fearnhead (2018, equation 13 of supplements);
(iii) is from Li & Fearnhead (2018, Lemma 5 and equation 13 of supplements); and (iv) is from
Li & Fearnhead (2018, Lemma 3 and Lemma 6).


Proof for Results in Section 3·1
e ∈ A | sobs + εn v) is the posterior
Proof of Lemma 1. For any fixed v ∈ Rd , recall that Π(θ
distribution given sobs + εn v with prior πδ (θ) and the misspecified model fen (· | θ). By Kleijn &
van der Vaart (2012), if there exist ∆n,θ0 and Vθ0 such that,
(KV1) for any compact set K ⊂ t(Bδ ),
sup log
t∈K

fen (sobs + εn v | θ0 + a−1
1
n t)
− tT Vθ0 ∆n,θ0 + tT Vθ0 t → 0,
2
fen (sobs + εn v | θ0 )

in probability as n → ∞, and
e n kθ − θ0 k > Mn | sobs + εn v)} → 0 as n → ∞ for any sequence of constants
(KV2) E{Π(a
Mn → ∞,
then

ˆ
e n (θ − θ0 ) ∈ A | sobs + εn v} −
sup Π{a

A∈B p

A

N (t; ∆n,θ0 , Vθ−1
) dt → 0,
0

in probability as n → ∞.
For (KV1), by the definition of fen (s | θ),
log

−2
−1
fen (sobs + εn v | θ0 + a−1
N {sobs + εn v; s(θ0 + a−1
n t)
n t), an A(θ0 + an t)}
= log
.
N {sobs + εn v; s(θ0 ), a−2
fen (sobs + εn v | θ0 )
n A(θ0 )}

As xT Ax − y T By = xT (A − B)x + (x − y)T B(x + y), for vectors x and y and matrices A
and B, by applying a Taylor expansion on s(θ0 + xt) and A(θ0 + xt) around x = 0, the right

19
hand side of above equation equals
a−1
n
T
−1
{Ds(θ0 + e(1)
t)t}
A(θ
)
ζ
(v,
t)
−
ζn (v, t)T
0
n
n
2
oT
a−1 n
+ n
D log A(θ0 + e(3)
t,
n t)
2

( p
X

)
Dθi A−1 (θ0 + e(2)
n t)ti

ζn (v, t)

i=1

(1)

(j)

where ζn (v, t) = A(θ0 )1/2 Wobs + an εn v − 12 Ds(θ0 + en t)t and for j = 1, 2, 3, en is a func(j)

tion of t satisfying |en | ≤ a−1
n which is from the remainder of the Taylor expansions. Since
Ds(θ), DA−1 (θ) and D log |A(θ)| are bounded in Bδ when δ is small enough,
sup log
t∈K

fen (sobs + εn v | θ0 + a−1
1
n t)
− tT I(θ0 )β0 {A(θ0 )1/2 Wobs + cε v} + tT I(θ0 )t → 0,
2
fen (sobs + εn v | θ0 )

in probability as n → ∞, for any compact set K. Therefore (KV1) holds with ∆n,θ0 =
β0 {A(θ0 )1/2 Wobs + cε v} and Vθ0 = I(θ0 ).
For (KV2), let rn (s | θ0 ) = αn {fn (s | θ0 ) − fen (s | θ0 )}. Since rn (s | θ0 ) is bounded by a
function integrable in Rd by Condition 4,
ˆ
e
e n kθ − θ0 k > Mn | s + εn v)fen (s | θ0 ) ds
E{Π(an kθ − θ0 k > Mn | sobs + εn v)} −
Π(a
d
R
ˆ
−1
≤αn
|rn (s | θ0 )| ds = o(1).
Rd

Then it is sufficient for the expectation under fen (s | θ0 ) to be o(1). For any constant M > 0,
with the transformation v̄ = an {s − s(θ0 )},
ˆ
e n kθ − θ0 k > Mn | s + εn v)fen (s | θ0 ) ds
Π(a
d
R
´
ˆ
ˆ
e(t, v̄ | v) dt
ktk>Mn π
´
N {v̄; 0, A(θ0 )} dv̄ +
N {v̄; 0, A(θ0 )} dv̄,
≤
e(t, v̄ | v) dt
kv̄k>M
kv̄k≤M
t(Bδ ) π
−1
−1
e
where π
e(t, v̄ | v) = πδ (θ0 + a−1
n t)fn {s(θ0 ) + an v̄ + εn v | θ0 + an t}. For the first term in the
above upper bound, it is bounded ´by a series which does not depend on M and is o(1) as Mn →
∞, as shown below. Obviously t(Bδ ) π
e(t, v̄ | v) dt can be lower bounded for some constant

mδ > 0. Choose δ small enough such that Ds(θ) and A(θ)1/2 are bounded for θ ∈ Bδ . Let λmin
and λmax be their common bounds. When kv̄k < M and Mn is large enough,


supθ∈Bδ kDs(θ)tk
{t : ktk > Mn } ⊂ t :
≥ kan εn v + v̄k .
(8)
2
Then since for any v̄ satisfying kv̄k < M , by a Taylor expansion,
−1
d
(1)
−1
fen {s(θ0 ) + a−1
n v̄ + εn v | θ0 + an t} = an N {Ds(θ0 + en t)t; v̄ + an εn v, A(θ0 + an t)},

π
e(t, v̄ | v) ≤ cN (λ−1
max λmin ktk/2; 0, 1), where c is some positive constant, for t in the right hand
side of (8). Then
´
ˆ
ˆ
e(t, v̄ | v) dt
ktk>Mn π
−1
´
N (λ−1
N {v̄; 0, A(θ0 )} dv̄ ≤ mδ c
max λmin ktk/2; 0, 1) dt,
π
e
(t,
v̄
|
v)
dt
kv̄k≤M
ktk>Mn
t(Bδ )

20
the right hand side of which is o(1) when Mn → ∞. Meanwhile by letting M → ∞, it can be
seen that the expectation under fen (s | θ0 ) is o(1). Therefore (KV2) holds and the lemma holds.
´
The following lemma is used for equations Rp gn (t, v) dt = |A(θ0 )|−1/2 Gn (v) and
´
−1/2 G(v).
Rp g(t, v) dt = |A(θ0 )|
L EMMA 8. For a rank-p d × p matrix A, a rank-d d × d matrix B and a d-dimension vector
c,

N (At; Bv + c, Id ) = N t; (AT A)−1 AT (c + Bv), (AT A)−1 g(v; A, B, c),

(9)

where P = AT A, and
g(v; A, B, c) =

1
(2π)(d−p)/2




1
T
T
−1 T
exp − (c + Bv) (I − A(A A) A )(c + Bv) .
2

Proof. This can be verified easily by matrix algebra.



The following lemma regarding the continuity of a certain form of integral will be helpful
when applying the continuous mapping theorem.
L EMMA 9. Let l1 , l10 , l2 , l20 and l3 be positive integers satisfying l10 ≤ l1 and l20 ≤ l2 . Let A and
B be l1 × l10 and l2 × l20 matrices, respectively, satisfying that AT A and B T B are positive definite. Let g1 (·), g2 (·) and g3 (·) be functions in Rl1 , Rl2 and Rl3 , respectively, that are integrable
and continuous almost everywhere. Assume:
(i) gj (·) is bounded in Rlj for j = 1, 2;
(ii) gj (w) depends on w only through kwk and is a decreasing function of kwk, for j = 1, 2;
and
´ Ql10 +l20 +l
(iii) there exists a non-negative integer l such that Rl3 k=1
wik g3 (w) dw < ∞ for any
coordinates (wi1 , . . . , wil0 +l0 +l ) of w.
1 2
Then the function,
ˆ ˆ ˆ
Pl (w1 , w2 , w3 ) |g1 (Aw1 + x1 w2 + x2 w3 + x3 ) − g1 (Aw1 )| g2 (Bw2 + x4 w3 + x5 )g3 (w3 ) dw3 dw2 dw1 ,
0

where x1 ∈ Rl1 ×l2 , x2 ∈ Rl1 ×l3 , x4 ∈ Rl2 ×l3 , x3 ∈ Rl1 and x5 ∈ Rl2 , is continuous almost everywhere.
Proof. Let mA and mB be the lower bound of A and B respectively. For any
0
(x01, . . . , x05 ) ∈ Rl1 ×l2 × Rl1 ×l3 × Rl2 ×l3 × Rl1 × Rl2 such that the integrand in the target
integral is continuous, consider any sequence (xn1 , . . . , xn5 ) converging to (x01, . . . , x05 ).
It is sufficient to show the convergence of the target function at (xn1 , . . . , xn5 ). Let
VA = {w1 : kAw1 k/2 ≥ sup(xn1 ,xn2 ,xn3 ) kxn1 w2 + xn2 w3 + xn3 k}, VB = {w2 : kBw2 k/2 ≥
sup(xn4 ,xn5 ) kxn4 w3 + xn5 k}, UA = {w1 : kw1 k ≤ 4m−1
A (kx01 w2 k + kx02 w3 k + kx03 k)} and
−1
UB = {w2 : kw2 k ≤ 4mB (kx04 w3 k + kx05 k)}. We have VAc ⊂ UA and VBc ⊂ UB . Then according to the following upper bounds and condition (iii),
|g1 (Aw1 + xn1 w2 + xn2 w3 + xn3 ) − g1 (Aw1 )| ≤ g1 (Aw1 + xn1 w2 + xn2 w3 + xn3 ) + g1 (Aw1 ),
g1 (Aw1 + xn1 w2 + xn2 w3 + xn3 ) ≤ ḡ1 (mA kw1 k/2)1{w1 ∈VA } + sup g1 (w)1{w1 ∈UA } ,
w∈Rl1

g2 (Bw2 + x4 w3 + x5 ) ≤ ḡ2 (mB kw2 k/2)1{w2 ∈VB } + sup g2 (w)1{w2 ∈UB } ,
w∈Rl2

21
where g1 (w) = ḡ1 (kwk) and g2 (w) = ḡ2 (kwk), by applying the dominated convergence theorem, the target function at (xn1, . . . , xn5 ) converges to its value at (x01, . . . , x05 ).

Proof of Lemma 2. The first part holds according to Lemma 5 of Li & Fearnhead (2018). For
the second part, when cε = ∞, by the transformation v 0 = v 0 (v, t),
ˆ

ˆ

ˆ

ˆ


Pl Ds(θ0 )t +

Pl (v)gn (t, v) dtdv =
Rd

Rd

t(Bδ )

t(Bδ )


1 0
1
v −
A(θ0 )1/2 Wobs gn0 (t, v 0 ) dtdv 0 .
a n εn
a n εn

By applying Lemma 9 and the continuous
´ ´mapping theorem in Lemma 6 to the right hand
side of the above when cε = ∞, and to Rd t(Bδ ) Pl (v)gn (t, v) dtdv when cε < ∞, and using
´
−1/2 G(v), the lemma holds.

Rp g(t, v) dt = |A(θ0 )|
Proof of Lemma 3. (a), (b) and the first part of (c) hold immediately by Lemma 7. The second
part of (c) is stated in the proof of Theorem 1 of Li & Fearnhead (2018).

L EMMA 10. Assume conditions 1–5.
e ε {an (θ − θeε ) ∈ A | sobs } have the
(i) If cε ∈ (0, ∞) then Πε {an (θ − θε ) ∈ A | sobs } and Π
same limit in distribution.
(ii) If cε = 0 or cε = 0∞ then
e ε {an,ε (θ − θeε ) ∈ A | sobs } = op (1).
supA∈Bp Πε {an,ε (θ − θε ) ∈ A | sobs } − Π
(iii) If Condition 6 holds then
e ε {an (θ∗ − θeε∗ ) ∈ A | sobs } = op (1).
supA∈Bp Πε {an (θ∗ − θε∗ ) ∈ A | sobs } − Π
Proof. Let λn = an,ε (θε − θeε ), and by Lemma 3(c), λn = op (1). When cε ∈ (0, ∞), for any
A ∈ B p , decompose Πε {an (θ − θε ) ∈ A | sobs } into the following three terms,
h

i
e ε {an (θ − θε ) ∈ A | sobs }
Πε {an (θ − θε ) ∈ A | sobs } − Π

h
i
e ε {an (θ − θeε ) ∈ A + λn | sobs } − Π
e ε {an (θ − θeε ) ∈ A | sobs }
+ Π
e ε {an (θ − θeε ) ∈ A | sobs }.
+Π
For (i) to hold, it is sufficient that the first two terms in the above are op (1). The first term is
e ε {an (θ − θeε ) ∈
op (1) by Lemma 3. For the second term to be op (1), given the leading term of Π
A | sobs } stated in the proof of Proposition 1 in the main text, it is sufficient that
ˆ
v∈Rd

ˆ 
−

sup
A+λn

N {t; µn (v), I(θ0 )−1 } dt = op (1).

A

´
´
This holds by noting that the left hand side of the above is bounded by ( A+λn − A )c dt for
some constant c and this upper bound is op (1) since λn = op (1). Therefore (i) holds.

22
e ε {an,ε (θ − θeε ) ∈ A | sobs }
When cε = 0 or ∞, supA∈Bp Πε {an,ε (θ − θε ) ∈ A | sobs } − Π
is bounded by
e ε {an,ε (θ − θε ) ∈ A | sobs }
supA∈Bp Πε {an,ε (θ − θε ) ∈ A | sobs } − Π
ˆ
e
e
ψ(t) dt
+supA∈Bp Πε {an,ε (θ − θε ) ∈ A + λn | sobs } −
A+λn
ˆ
e ε {an,ε (θ − θeε ) ∈ A | sobs } −
+supA∈Bp Π
ψ(t) dt
A
ˆ
ˆ
ψ(t) dt .
+supA∈Bp
ψ(t) dt −

(10)

A

A+λn

With similar arguments as before, the first three
´ terms are op (1). For the fourth term, by transforming t to t + λn , it is upper bounded by Rp |ψ(t − λn ) − ψ(t)| dt which is op (1) by the
continuous mapping theorem. Therefore (ii) holds.
For (iii), the left hand side of the equation has the decomposed upper bound similar to (10),
with θ, θε , θeε and ψ(t) replaced by θ∗ , θε∗ , θeε∗ and N {t; 0, I(θ0 )−1 }. Then by Lemma 5, using
the leading term of Πε {an (θ∗ − θε∗ ) ∈ A | sobs } stated in the proof of Theorem 1, and similar
arguments to those used for the fourth term of (10), it can be seen that this upper bound is op (1).
Therefore (iii) holds.

Proof for Results in Section 3·2
To prove Lemmas 4 and 5, some notation regarding the regression adjusted approximate
Bayesian computation posterior, similar to those defined previously, are needed. Consider transp
p
d
formations t´= t(θ)
´ and v = v(s). For A ⊂ R and the scalar function h(t, v) in R × R , let
π
eA,tv (h) = t(A) Rd h(t, v)e
πε,tv (t, v) dvdt.
Proof of Lemma 4. Since βε = covε (θ, s)varε (s)−1 , to evaluate the covariance matrices, we
need to evaluate πRp {(θ − θ0 )k1 (s − sobs )k2 }/πRp (1) for (k1 , k2 ) = (0, 0), (1, 0), (1, 1), (0, 1)
and (0, 2).
First of all, we show that πBδc {(θ − θ0 )k1 (s − sobs )k2 } is ignorable for any δ < δ0 by showing
αδ

that it is Op (e−an,ε cδ ) for some positive constants cδ and αδ . By dividing Rd into {v : kεn vk ≤
δ 0 /3} and its complement,
ˆ
(s − sobs ) fn (s | θ)K

sup
θ∈Bδc



k2

Rd

θ∈Bδc

+



ε−d
n ds

ˆ

(
≤ sup

s − sobs
εn

sup
ks−sobs k≤δ 0 /3

(s − sobs )k2 K

fn (s | θ)

0
−d
K{λmin (Λ)ε−1
n δ /3}εn

ˆ

Rd



s − sobs
εn

(s − sobs )k2 fn (s | θ) ds.



)
ε−d
n ds
(11)

Rd

By Condition 2(ii), Condition 6 and following the arguments in the proof of Lemma 3 of Li &
αδ
Fearnhead (2018), the right hand side of (11) is Op (e−an,ε cδ ), which is sufficient for πBδc {(θ −
αδ

θ0 )k1 (s − sobs )k2 } to be Op (e−an,ε cδ ).

23
For the integration over Bδ , by Lemma 7 (ii),

π
eBδ ,tv (tk1 v k2 )
πBδ {(θ − θ0 )k1 (s − sobs )k2 }
−k1 k2
= an,ε εn
+
πBδ (1)
π
eBδ ,tv (1)
´
´ k k
)
−1 t)K(v) dvdt
1 v 2 π(θ + a−1 t)r (s
t
+
ε
v
|
θ
+
a
0
n
n
0
obs
n,ε
n,ε
t(B
)
δ
{1 + Op (αn−1 )}
αn−1
π
eBδ ,tv (1)
where rn (s | θ) is the scaled remainder αn {fn (s | θ) − fen (s | θ)}. In the above, the second term
in the first brackets is Op (αn−1 ) by the proof of Lemma 6 of Li & Fearnhead (2018). Then


π
eBδ ,tv (tk1 v k2 )
πBδ {(θ − θ0 )k1 (s − sobs )k2 }
−k1 k2
−1
= an,ε εn
+ Op (αn ) ,
πBδ (1)
π
eBδ ,tv (1)
and the moments π
eBδ ,tv (tk1 v k2 )/e
πBδ ,tv (1) need to be evaluated. Theorem 1 of Li & Fearnhead
(2018) gives the value of π
eBδ ,tv (t)/e
πBδ ,tv (1), and this is obtained by substituting the leading
term of π
eε,tv (t, v), that is π(θ0 )gn (t, v) as stated in Lemma 2, into the integrands. The other
moments can be evaluated similarly, and give
 −1
bn β0 {A(θ0 )1/2 Wobs + an εn EGn (v)},
(k1 , k2 ) = (1, 0),




k
k
1/2 W
T )}, (k , k ) = (1, 1),
π
eBδ ,tv (t 1 v 2 )
b−1
β
{A(θ
)
E
(v]
+
a
ε
E
(vv
0
0
n n Gn
1 2
obs Gn
n
=

π
eBδ ,tv (1)
E
(v),
(k
1 , k2 ) = (0, 1),
Gn



T
EGn (vv ),
(k1 , k2 ) = (0, 2),
2 4
+ Op (a−1
n,ε ) + Op (an εn ),

(12)

where bn = 1 when cε < ∞, and an εn when cε = ∞. By Lemma 2, EGn (vv T ) = Θp (1). Since
−2/5
−2/5
αn−1 = o(an ), covε (θ, s) =ε2n β0 varGn (v) + op (an ε2n ) and varε (s) = ε2n varGn (v){1 +
−2/5
op (an )}. Thus
βε = β0 + op (a−2/5
),
n
and the lemma holds.

´ ´

(13)


−d
For A ⊂ Rp and B ⊂ Rd , let π(A, B) = A B π(θ)fn (s | θ)K{ε−1
n (s − sobs )}εn dsdθ and
´ ´
−d
π
e(A, B) = A B π(θ)fen (s | θ)K{ε−1
n (s − sobs )}εn dsdθ. Denote the marginal mean values
of s for πε (θ, s | sobs ) and π
eε (θ, s | sobs ) by sε and seε respectively.

Proof of Lemma 5. For
(a),
write
Πε (θ∗ ∈ Bδc | sobs )
as
π[Rp , {s : θ∗ (θ, s) ∈
d−p
p
d
By Lemma 7, π(R , R ) = πP (1) = Θp (an,ε ). By the triangle inequality,

Bδc }]/π(Rp , Rd ).

c
π[Rp , {s : θ∗ (θ, s) ∈ Bδc }] ≤ π(Bδ/2
, Rd ) + π[Bδ/2 , {s : kβε (s − sobs )k ≥ δ/2}],

(14)

and it is sufficient that the right hand side of the above inequality is op (1). Since its first term is
c (1), by Lemma 7 the first term is op (1).
πBδ/2
−7/5

−7/5

When εn = Ω(an ) or Θ(an ), by (13), βε − β0 = op (1) and so βε is bounded in probability. For any constant βsup > 0 and β ∈ Rp×d satisfying β ≤ βsup ,


δ
π[Bδ/2 , {s : kβ(s − sobs )k ≥ δ/2}] ≤ K ε−1
ε−d
n ,
2βsup
and by Condition 2(iv), the second term in (14) is op (1).

24
−7/5

When εn = o(an ), βε is unbounded and the above argument does not apply. Let δ1 be a
constant less than δ0 such that inf θ∈Bδ1 /2 λmin {A(θ)−1/2 } ≥ m and inf θ∈Bδ1 /2 λmin {Ds(θ)} ≥
m for some positive constant m. In this case, it is sufficient to consider δ < δ1 . By Condition 4,
rn (s | θ) ≤ adn |A(θ)|1/2 rmax [an A(θ)−1/2 {s − s(θ)}].
Using the transformation t = t(θ) and v = v(s), fn (s | θ) = fen (s | θ) + αn−1 rn (s | θ) and applying the Taylor expansion of s(θ0 + xt) around x = 0,
π[Bδ/2 , {s : kβε (s − sobs )k ≥ δ/2}] ≤
ˆ
ˆ
−1/2
1/2
c
N [A(θ0 + a−1
{Ds(θ0 + e(1)
Wobs − an εn v}; 0, Id ]K(v) dvdt
n t)
n t)t − A(θ0 )
t(Bδ/2 )

ˆ

kβε εn vk≥δ/2

ˆ

+c
t(Bδ/2 )

kβε εn vk≥δ/2

−1/2
1/2
rmax [A(θ0 + a−1
{Ds(θ0 + e(1)
Wobs − an εn v}]K(v) dvdt,
n t)
n t)t − A(θ0 )

for some positive constant c. To show that the right hand side of the above inequality is op (1),
consider a function g4 (·) in Rd satisfying that g4 (v) can be written as g 4 (kvk) and g 4 (·) is
−1/2 , C (t) = Ds(θ + ξ ) and c = A(θ )1/2 W
decreasing. Let An (t) = A(θ0 + a−1
n
0
1
0
obs . For
n t)
p
c
each n divide R into Vn = {t : kCn (t)tk/2 ≥ kc + an εn vk} and Vn . In Vn , kAn (t){Cn (t)t −
c − an εn v}k ≥ m2 ktk/2 and in V c , ktk ≤ 2m−1 kc + an εn vk. Then
ˆ
ˆ
g4 [An (t){Cn (t)t − c − an εn v}]K(v) dvdt
t(Bδ/2 )

ˆ

kβε εn vk≥δ/2

(ˆ

≤
kβε εn vk≥δ/2

where

´
c
Vn

ˆ
2

Rp

g 4 (m ktk/2) dt + sup g4 (v)
v∈Rp

)
1 dt K(v) dv,

c
Vn

1 dt is the volume of Vnc in Rp . Then since βε εn = op (1), an εn = op (1) and

´
Vnc

1 dt

vkp , the right hand side of the above inequality is o

is proportional to kc + an εn
p (1). This implies
π(Bδ/2 , {s : kβε (s − sobs )k ≥ δ/2}) = op (1).
e ε (θ∗ ∈ B c | sobs ), since the supTherefore in both cases Πε (θ∗ ∈ Bδc | sobs ) = op (1). For Π
δ
e ε (θ∗ ∈ B c | sobs ) = 0.
port of its prior is Bδ , there is no probability mass outside Bδ , i.e. Π
δ
Therefore (a) holds.
For (b),
e ε (θ∗ ∈ Aθ ∩ Bδ | sobs )
supA∈Bp Πε (θ∗ ∈ Aθ ∩ Bδ | sobs ) − Π
supA∈Bp |π(Rp , {s : θ∗ (θ, s) ∈ Aθ ∩ Bδ }) − π
e(Rp , {s : θ∗ (θ, s) ∈ Aθ ∩ Bδ })|
+ op (1)
π
eBδ (1)
´ ´
−d
π(θ)|rn (s | θ)|K{ε−1
n (s − sobs )}εn dsdθ
−1 Bδ Rd
≤αn
+ op (1).
π
eBδ (1)
=

Then by the proof of Lemma 6 of Li & Fearnhead (2018), (b) holds.
For (c), to begin with, an (θε∗ − θeε∗ ) = an (θε − θeε ) − an βε (sε − seε ). By Lemma 7, an (θε −
θeε ) = op (1). For an βε (sε − seε ), similar to the arguments of the proof of Lemma 4,


π
eBδ,tv (v)
π
eB (v)
−1
sε − sobs = εn
+ Op (αn ) {1 + Op (αn−1 )}, seε − sobs = εn δ,tv
{1 + Op (αn−1 )}.
π
eBδ,tv (1)
π
eBδ,tv (1)

25
−3/5

Then an βε (sε − seε ) = Op (αn−1 an εn ) which is op (1) if εn = o(an ). Therefore the first part
of (c) holds. Since θeε∗ = θeε − βε (e
sε − sobs ), by the expansion of θeε in Lemma 3(c), the above
expansion of seε − sobs and (12), the second part of (c) holds.

Proof for Results in Section 3·3
Proof of Theorem 2. The integrand of pacc,q is similar to that of πRp (1). The expansion of
πRp (1) is given in Lemma 7(ii), and following the same reasoning, pacc,q can be expanded as
´ ´
εdn Bδ Rd qn (θ)fe(sobs + εn v | θ)K(v) dvdθ{1 + op (1)}. With transformation t = t(θ), plugging the expression of qn (θ) and π
eε,tv (t, v) gives that
ˆ
π
eε,tv (t, v)
−1
(rn,ε )−p q(rn,ε
t − cµ )
pacc,q = (an,ε εn )d
dvdt{1 + op (1)},
πδ (θ0 + a−1
t(Bδ )
n,ε t)
where rn,ε = σn /a−1
n,ε and cµ,n = σn (µn − θ0 ). By the assumption of µn , denote the limit of
cµ,n by cµ . Then by Lemma 2, pacc,q can be expanded as
ˆ
−1
pacc,q = (an,ε εn )d
(rn,ε )−p q(rn,ε
t − cµ,n )gn (t, v) dvdt{1 + op (1)}.
(15)
t(Bδ )×Rd

Denote the leading term of the above by Qn,ε .
For (1), when cε = 0, since supt∈Rp gn (t, v) ≤ c1 K(v) for some positive constant c1 , Qn,ε
is upper bounded by (an εn )d c1 almost surely. Therefore pacc,q → 0 almost surely as n → ∞.
When rn,ε → ∞, since q(·) is bounded
in Rp by some positive constant c2 , Qn,ε is upper
´
−p
d
bounded by (rn,ε ) c2 (an,ε εn ) Rp ×Rd gn (t, v) dvdt. Therefore pacc,q → 0 in probability as
´
n → ∞ since Rp ×Rd gn (t, v) dvdt = Θp (1) by Lemma 2.
−1 t(θ) − c
For (2), let t̃(θ) = rn,ε
µ,n and t̃(A) be the set {φ : φ = t̃(θ) for some θ ∈ A}. Since
−1
t̃ = σn (θ − θ0 ) − cµ,n and σn−1 → ∞, t̃(Bδ ) converges to Rp in probability as n → ∞. With
the transformation t̃ = t̃(θ),
(
´
(an εn )d t̃(Bδ )×Rd q(t̃)gn {rn,ε (t̃ + cµ,n ), v} dt̃dv, cε < ∞,
Qn,ε = ´
0
0
0
cε = ∞.
t̃(Bδ )×Rd q(t̃)gn {rn,ε (t̃ + cµ,n ), v } dt̃dv ,
By Lemma 9 and the continuous mapping theorem,
( ´
cd p d q(t̃)g{r1 (t̃ + cµ ), v} dt̃dv,
Qn,ε → ´ε R ×R
Rp ×Rd q(t̃)g{r1 (t̃ + cµ ), v} dt̃dv,

cε < ∞,
cε = ∞,

in distribution as n → ∞. Since the limits above are Θp (1), pacc,q = Θp (1).
´ For (3), when cε = ∞ and r1 = 0, in the above, the limit of Qn,ε in distribution is

Rp ×Rd q(t̃)g(0, v) dt̃dv = 1. Therefore pacc,q converges to 1 in probability as n → ∞.

R EFERENCES
B EAUMONT, M. A. (2010). Approximate Bayesian computation in evolution and ecology. Annual Review of Ecology,
Evolution, and Systematics 41, 379–406.
B EAUMONT, M. A., C ORNUET, J.-M., M ARIN , J.-M. & ROBERT, C. P. (2009). Adaptive approximate Bayesian
computation. Biometrika 96, 983–990.
B EAUMONT, M. A., Z HANG , W. & BALDING , D. J. (2002). Approximate Bayesian computation in population
genetics. Genetics 162, 2025–2035.
B LUM , M. G. (2010). Approximate Bayesian computation: a nonparametric perspective. Journal of the American
Statistical Association 105, 1178–1187.

26
B ONASSI , F. V. & W EST, M. (2015). Sequential Monte Carlo with adaptive weights for approximate Bayesian
computation. Bayesian Analysis 10, 171–187.
C ALVET, L. E. & C ZELLAR , V. (2015). Accurate methods for approximate Bayesian computation filtering. Journal
of Financial Econometrics 13, 798–838.
C HERNOZHUKOV, V. & H ONG , H. (2003). An MCMC approach to classical estimation. Journal of Econometrics
115, 293–346.
C ORNUET, J.-M., S ANTOS , F., B EAUMONT, M. A., ROBERT, C. P., M ARIN , J.-M., BALDING , D. J., G UILLE MAUD , T. & E STOUP, A. (2008). Inferring population history with DIY ABC: a user-friendly approach to approximate Bayesian computation. Bioinformatics 24, 2713–2719.
D ROVANDI , C. C. & P ETTITT, A. N. (2011). Estimation of parameters for macroparasite population evolution using
approximate Bayesian computation. Biometrics 67, 225–233.
D ROVANDI , C. C., P ETTITT, A. N. & L EE , A. (2015). Bayesian indirect inference using a parametric auxiliary
model. Statistical Science 30, 72–95.
D UFFIE , D. & S INGLETON , K. J. (1993). Simulated moments estimation of Markov models of asset prices. Econometrica 61, 929–952.
F EARNHEAD , P. & P RANGLE , D. (2012). Constructing summary statistics for approximate Bayesian computation:
semi-automatic approximate Bayesian computation (with Discussion). Journal of the Royal Statistical Society:
Series B (Statistical Methodology) 74, 419–474.
F ILIPPI , S., BARNES , C. P., C ORNEBISE , J. & S TUMPF, M. P. (2013). On optimality of kernels for approximate
Bayesian computation using sequential Monte Carlo. Statistical Applications in Genetics and Molecular Biology
12, 87–107.
F RAZIER , D. T., M ARTIN , G. M., ROBERT , C. P. & ROUSSEAU , J. (2016). Asymptotic Properties of Approximate
Bayesian Computation. arXiv.1607.06903 .
G OURI ÉROUX , C. & RONCHETTI , E. (1993). Indirect inference. Journal of Applied Econometrics 8, s85–s118.
H EGGLAND , K. & F RIGESSI , A. (2004). Estimating functions in indirect inference. Journal of the Royal Statistical
Society, Series B (Statistical Methodology) 66, 447–462.
K LEIJN , B. & VAN DER VAART, A. (2012). The Bernstein-von-Mises theorem under misspecification. Electronic
Journal of Statistics 6, 354–381.
L ENORMAND , M., JABOT, F. & D EFFUANT, G. (2013). Adaptive approximate Bayesian computation for complex
models. Computational Statistics 28, 2777–2796.
L I , W. & F EARNHEAD , P. (2018). On the asymptotic efficiency of approximate Bayesian computation estimators.
Biometrika , to appear.
M ARIN , J.-M., P ILLAI , N. S., ROBERT, C. P. & ROUSSEAU , J. (2014). Relevant statistics for Bayesian model
choice. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 76, 833–859.
M ARJORAM , P., M OLITOR , J., P LAGNOL , V. & TAVARE , S. (2003). Markov chain Monte Carlo without likelihoods.
Proceedings of the National Academy of Sciences 100, 15324–15328.
PAULI , F., R ACUGNO , W. & V ENTURA , L. (2011). Bayesian composite marginal likelihoods. Statistica Sinica 21,
149–164.
R IBATET, M., C OOLEY, D. & DAVISON , A. C. (2012). Bayesian inference from composite likelihoods, with an
application to spatial extremes. Statistica Sinica 22, 813–845.
RULI , E., S ARTORI , N. & V ENTURA , L. (2016). Approximate Bayesian computation with composite score functions. Statistics and Computing 26, 679–692.
S ERFLING , R. J. (2009). Approximation theorems of mathematical statistics, vol. 162. John Wiley & Sons.
S OUBEYRAND , S. & H AON -L ASPORTES , E. (2015). Weak convergence of posteriors conditional on maximum
pseudo-likelihood estimates and implications in ABC. Statistics & Probability Letters 107, 84–92.
T ONI , T., W ELCH , D., S TRELKOWA , N., I PSEN , A. & S TUMPF, M. P. (2009). Approximate Bayesian computation
scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society Interface
6, 187–202.
VAN DER VAART, A. W. (2000). Asymptotic Statistics. Cambridge University Press, Cambridge.
VARIN , C., R EID , N. & F IRTH , D. (2011). An overview of composite likelihood methods. Statistica Sinica 21, 5–42.
W EGMANN , D., L EUENBERGER , C. & E XCOFFIER , L. (2009). Efficient approximate Bayesian computation coupled
with Markov chain Monte Carlo without likelihood. Genetics 182, 1207–1218.
W OOD , S. N. (2010). Statistical inference for noisy nonlinear ecological dynamic systems. Nature 466, 1102–1104.

