Overview of Approximate Bayesian Computation
S. A. Sisson Y. Fan M. A. Beaumont


Discovering Key Nodes in a Temporal Social Network
    Jinshuo LIU1 
Chenghao MOU2
Donghong JI1 ()
    December 30, 2023
==================================================







§ INTRODUCTION






In Bayesian inference, complete knowledge about a vector of model parameters, θ∈Θ, obtained by fitting a model ℳ, is contained in the posterior distribution.
Here, prior beliefs about the model parameters as expressed through the prior distribution, π(θ), are updated by observing data y_obs∈𝒴 through the likelihood function π(y_obs|θ) of the model. Using Bayes' Theorem, the resulting posterior distribution

    π(θ|y_obs)=p(y_obs|θ)π(θ)/∫_Θ p(y_obs|θ)π(θ)dθ,

contains all necessary information required for analysis of the model, including model checking and validation, predictive inference and decision making. Typically, the complexity of the model and/or prior means that the posterior distribution, π(θ|y_obs), is not available in closed form, and so numerical methods are needed to proceed with the inference. A common approach makes use of Monte Carlo integration to enumerate the necessary integrals. This relies on the ability to draw samples θ^(1),θ^(2),…,θ^(N)∼π(θ|y_obs) from the posterior distribution so that a finite sample approximation to the posterior is given by the empirical measure

    π(θ|y_obs)≈1/N∑_i=1^Nδ_θ^(i)(θ),

where δ_Z(z) denotes the Dirac measure, defined as δ_Z(z)=1 if z∈ Z and δ_Z(z)=0 otherwise.
As the size of the sample from the posterior gets large, then the finite sample approximation better approximates the true posterior so that lim_N→∞1/N∑_i=1^Nδ_θ^(i)(θ) →π(θ|y_obs), by the law of large numbers.
As a result, the expectation of a function a(θ) under π(θ|y_obs) can be estimated as

    𝔼_π[a(θ)]    =   ∫_Θ a(θ)π(θ|y_obs)dθ
       ≈   ∫_Θ a(θ)1/N∑_i=1^Nδ_θ^(i)(θ)dθ
    	= 1/N∑_i=1^Na(θ^(i)).




There are a number of popular algorithms available for generating samples from posterior distributions, such as importance sampling, Markov chain Monte Carlo (MCMC) and sequential Monte Carlo (SMC) brooksgjm11,chen+si00,doucet+dg01,delmoral+dj06.


Inherent in such Monte Carlo algorithms is the need to numerically evaluate the posterior distribution, π(θ|y_obs), up to a normalisation constant, commonly many thousands or millions of times. For example, in the Metropolis-Hastings algorithm, an MCMC algorithm, this arises through computing the probability that the Markov chain accepts the proposed move from a current point θ to a proposed point θ'∼ q(θ,θ') where q is some proposal density, given by  α(θ,θ')=min{1,π(θ'|y_obs)q(θ',θ)/π(θ|y_obs)q(θ,θ')}. Similarly in SMC algorithms, the incremental particle weight is given by w_t(θ_t)=π_t(θ_t|y_obs)L_t-1(θ_t,θ_t-1)/π_t-1(θ_t-1|y_obs)M_t(θ_t-1,θ_t), where M_t and L_t-1 are transition kernels, and π_t denotes a function strongly related to the posterior distribution, such as π_t(θ_t|y_obs)=[π(θ_t|y_obs)]^t/Tπ(θ_t)^1-t/T. Evaluating acceptance probabilities or particle weights clearly requires evaluation of the likelihood function.


However, for an increasing range of scientific problems – see Section <ref> for a selection – numerical evaluation of the likelihood function, π(y_obs|θ), is either computationally prohibitive, or simply not possible. Examples of the former can occur where the size of the observed dataset, y_obs, is sufficiently large that, in the absence of low dimensional sufficient statistics, evaluating the likelihood function even once is impracticable.  This can easily occur in the era of Big Data, for example, through large genomic datsets.

Partial likelihood intractability can arise, for instance, in models for Markov random fields. Here, the likelihood function can be written as
p(y_obs|θ) = 1/Z_θp̃(y_obs|θ)
where p̃(y_obs|θ) is a function that can be evaluated, and where the normalisation constant, Z_θ=∑_𝒴p̃(y|θ), depends on the parameter vector θ.
Except for trivial datasets, the number of possible data configurations in the set 𝒴 means that brute-force enumeration of Z_θ is typically infeasible grelaud+rmrt09,moller+prb06.

While there are algorithmic techniques available that arrange for the intractable normalising constants to cancel out within e.g. Metropolis-Hastings acceptance probabilities moller+prb06, or that numerically approximate Z_θ through e.g. path sampling or thermodynamic integration, these are not viable when p̃(y|θ) itself is also computationally intractable.

Instances when the complete likelihood function is unavailable can also occur when the model density function is only implicitly defined, for example, through quantile or characteristic functions drovandi+p11,peters+sf12. Similarly, the likelihood function may only be implicitly defined as a data generation process.



In these scenarios, if the preferred model is computationally intractable, the need to repeatedly evaluate the posterior distribution to draw samples from the posterior makes the implementation of standard Bayesian simulation techniques impractical.
Faced with this challenge,
one option is simply to fit a different model that is more amenable to statistical computations. The disadvantage of this approach is that the model could then be less realistic, and not permit inference on the particular questions of interest for the given analysis. A more attractive alternative, may be to consider an approximation to the preferred model, so that modelling realism is maintained at the expense of some approximation error. While various posterior approximation methods are available, “likelihood-free” Bayesian methods, of which approximate Bayesian computation (ABC) is a particular case, have emerged as an effective and intuitively accessible way of performing an approximate Bayesian analysis.

In this Chapter, we aim to give an intuitive exploration of the basics of ABC methods, illustrated wherever possible by simple examples. The scope of this exploration is deliberately limited, for example, we focus only on the use of simple rejection sampling based ABC samplers, in order that this Chapter will provide an accessible introduction to a subject which is given more detailed and advanced treatments in the rest of this Handbook.







§ LIKELIHOOD-FREE INTUITION




The basic mechanism of likelihood-free methods can be fairly easily understood at an intuitive level.
For the moment, we assume that data generated under the model, y∼ p(y|θ), are discrete.
Consider the standard rejection sampling algorithm for sampling from a density f(θ):



If we specify f(θ)=π(θ|y_obs), and suppose that
the prior is used as the sampling distribution,
then the acceptance probability is proportional to the likelihood, as then f(θ)/Kg(θ)∝ p(y_obs|θ).  While direct evaluation of this acceptance probability is not available if the likelihood is computationally intractable, it is possible to stochastically determine whether or not to accept or reject a draw from the sampling density, without numerical evaluation of the acceptance probability. The following discussion assumes that the data y are discrete (this will be relaxed later).



This can be achieved by noting that the acceptance probability is proportional to the probability of generating the observed data, y_obs,  under the model p(y|θ) for a fixed parameter vector, θ. That is, suitably normalised, the likelihood function p(y|θ) can be considered as a probability mass function for the data. Put another way, for fixed θ, if we generate a dataset from the model y∼ p(y|θ), then the probability of generating our observed dataset exactly, so that y=y_obs, is precisely p(y_obs|θ). From this observation, we can use the Bernoulli event of generating y=y_obs (or not) to determine whether to accept (or reject) a draw from the sampling distribution, in lieu of directly evaluating the probability p(y_obs|θ).

This insight permits a rewriting of the simple rejection sampling algorithm, as given below.  A critical aspect of this modified algorithm is that it does not require numerical evaluation of the acceptance probability (i.e. the likelihood function). Note that if sampling is from g(θ) rather than the prior π(θ), then the acceptance probability is proportional to p(y_obs|θ)π(θ)/g(θ).  In this case, deciding whether to accept a draw from g(θ) can be split into two stages: firstly, as before, if we generate y∼ p(y|θ) such that y≠ y_obs then we reject the draw from g(θ). If however, y=y_obs, then we accept the draw from g(θ) with probability 
π(θ)/[Kg(θ)], where K≥max_θ f(θ)/g(θ). (These two steps may be interchanged so that the step with the least computational overheads is performed first.) Importance sampling versions of this and later algorithms are examined in fan+s18.








§ A PRACTICAL ILLUSTRATION: STEREOLOGICAL EXTREMES





In order to illustrate the performance of the likelihood-free rejection sampling algorithm, we perform a re-analysis of a stereological dataset with a computationally intractable model first developed by bortot+cs07.



 §.§ Background and model


Interest is in the distribution of the size of inclusions, microscopically small particles introduced during the production of steel. The steel strength is thought to be directly related to the size of the largest inclusion. 

Commonly, the sampling of inclusions involves measuring the maximum cross-sectional diameter of each observed inclusion, y_obs=(y_obs,1, …, y_obs,n)^⊤, obtained from a two-dimensional planar slice through the steel block. Each cross-sectional inclusion size is greater than some measurement threshold, y_obs,i>u.

The inferential problem is to analyse the unobserved distribution of the largest inclusion in the block, based on the information in the cross-sectional slice, y_obs. The focus on the size of the largest inclusion means that this is an extreme value variation on the standard stereological problem <cit.>.

Each observed cross-sectional inclusion diameter, y_obs,i, is associated with an unobserved true inclusion diameter V_i.
anderson+c02 proposed a mathematical model 
assuming that the inclusions were spherical with diameters V, and that their centres followed a homogeneous Poisson process with rate λ>0 in the volume of steel.
The distribution of the largest inclusion diameters, V|V>v_0
was assumed to follow a generalised Pareto distribution, with distribution function

    (V≤ v|V>v_0) = 1-[1+ξ(v-v_0)/σ]^-1/ξ_+,

for v>v_0, where [a]_+=max{0,a}, 
following standard extreme value theory arguments <cit.>. 
However, the probability of observing the cross-sectional diameter y_obs,i (where y_obs,i≤ V_i) is dependent on the value of V_i, as larger inclusion diameters give a greater chance that the inclusion will be observed in the two-dimensional planar cross-section. This means that the number of observed inclusions, n, is also a random variable.
Accordingly the parameters of the full spherical inclusion model are  θ=(λ, σ, ξ)^⊤.

anderson+c02 were able to construct a tractable likelihood function for this model by adapting the solution to Wicksell's corpuscle problem <cit.>. 
However, while their model assumptions of a Poisson process 
are not unreasonable, the assumption that the inclusions are spherical is not plausible in practice.



bortot+cs07 generalised this model to a family of ellipsoidal inclusions. 

While this model is more realistic than the spherical inclusion model, there are analytic and computational difficulties in extending likelihood-based inference to more general families of inclusion baddeley+j05,bortot+cs07. As a result ABC methods are a good candidate procedure to approximate the posterior distribution in this case.





 §.§ Analysis





For simplicity, suppose that we are interested in the spherical inclusions model, so that the true posterior distribution can be estimated directly. Suppose also that  the parameters of the generalised Pareto distribution are known to be σ=1.5 and ξ=0.1, so that interest is in the Poisson rate parameter, λ, only. In this setting, a sufficient statistic for the rate parameter is n_obs, the observed number of inclusions, so that π(θ|y_obs)=π(λ|n_obs) is the distribution of interest. 

Accordingly we can replace y_obs=n_obs in the likelihood-free rejection sampling algorithm.

For the dataset considered by bortot+cs07, n_obs=112.



Figure <ref>(a) shows scaled density estimates of π(λ|n_obs) (solid lines) obtained using the likelihood-free rejection sampling algorithm, for varying numbers of observed inclusions, n_obs=92, 102, 112, 122 and 132. As the observed number of inclusions increases, accordingly so does the location and scale of the posterior of the rate parameter. The dashed lines in Figure <ref>(a) denote the same density estimates of π(λ|n_obs), but obtained using a conditional version of the standard MCMC sampler developed by anderson+c02, which makes use of numerical evaluations of the likelihood. These estimates are known to correspond to the true posterior. The likelihood-free rejection algorithm estimates clearly coincide with the true posterior distribution.



The density estimates obtained under the likelihood-free algorithm are each based on approximately 25,000 accepted samples, obtained from 5 million draws from the U(0,100) prior. That is, the acceptance rate of the algorithm is approximately 0.5%.
This algorithm is clearly very inefficient, with the computational overheads being partially influenced by the mismatch between prior and posterior distributions, but they are primarily dominated by the probability of generating data from the model that exactly matches the observed data, n_obs. This is the price for avoiding likelihood evaluation.

On balance, the computational inefficiency is practically acceptable for this specific case. However, this raises the question of how viable this approach will be for more complex analyses, when the probability of generating data such that y=y_obs becomes even lower. Further, the acceptance probability will be exactly zero if the data generated under the model, y∼ p(y|θ), are continuous, which is likely to be the case in general.

In order to alleviate such computational overheads,
one possible variation on the likelihood-free rejection algorithm would be to adjust the potentially very low (or zero) probability requirement that y=y_obs exactly. Instead,  the acceptance criterion could require that the generated data is simply “close” to the observed data. For example, this might require that y-y_obs≤ h for some h≥ 0 and distance measure ·, such as Euclidean distance. This would also permit a relaxation of our previous assumption that data generated under the model, y∼ p(y|θ), are discrete. In this way, step <ref> of the Sampling stage of the likelihood-free rejection algorithm would become:



Of course, the output samples would no longer be draws from π(θ|y_obs) unless h=0, but will instead be draws from an approximation of π(θ|y_obs).

The logic behind this modification is that increasing h will considerably improve the acceptance rate of the algorithm. The hope is that, if h remains small, then the resulting estimate of the posterior will still be close to the true posterior. An illustration of this is shown in Figure <ref>(b), which shows density estimates obtained using the adjusted requirement that n-n_obs≤ h for h=0 (i.e. n=n_obs), 10 and 20. Computationally there is a marked improvement in algorithmic efficiency: the low 0.5% acceptance rate for h=0 increases to 10.5% and 20.5% for h=10 and 20 respectively.

However, there are now some clear deviations in the density estimate resulting from the likelihood-free algorithm, compared to the actual posterior, π(λ|n_obs) (solid lines). In fact, it is more accurate to refer to these density estimates as an approximation of the posterior. On one hand, the location and shape of the density are broadly correct, and for some applications, this level of approximation may be adequate. On the other hand, however, the scale of the approximation is clearly overestimated for larger values of h. Intuitively this makes sense: the adjusted criterion y-y_obs≤ h accepts θ∼ g(θ) draws if the generated data y is merely “close” to y_obs. As such, for many values of θ where it was previously very unlikely to generate data such that y=y_obs, it may now be possible to satisfy the more relaxed criterion. This will accordingly result in a greater range of θ values that will be accepted, and thereby increase the variability of the posterior approximation. The more relaxed the criterion (i.e. the larger the value of h), the greater the resulting variability.

It is possible to be more precise about the exact form of the posterior obtained through this adjusted procedure – this will be discussed in detail in the next Section. However, for this particular analysis, based on samples λ^(1),…,λ^(N) and datasets n^(1),…,n^(N) obtained from the likelihood-free rejection algorithm, it can be seen that as the posterior approximation is constructed from those values of λ=λ^(i) such that n^(i)-n_obs≤ h, then the posterior approximation can firstly be expressed as

    π̂(λ|n_obs)  = 1/N∑_i=1^Nδ_λ^(i)(λ)
    	   =   1/N∑_λ^(i):n^(i)-n_obs≤ hδ_λ^(i)(λ)
        =    ∑_h'=-h^*^h^*(1/N∑_λ^(i):(n^(i)-n_obs)=h'δ_λ^(i)(λ)),

where h^* is the largest integer such that h^*≤ h.
It then follows that

    lim_N→∞π̂(λ|n_obs)
    	 =
    	∑_h'=-h^*^h^*(n=n_obs+h')π(λ|n_obs+h').

That is, the “likelihood-free” approximation  of the posterior, π(θ|y_obs), is precisely an average of the individual posterior distributions π(λ|n_obs+h') for h'=-h^*,…,h^*, weighted according to (n=n_obs+h'), the probability of observing the dataset, n_obs+h', based on samples drawn from the (prior predictive) distribution p(n|λ)π(λ).
This can be loosely observed from Figure <ref>, in which the approximations for h=10  and h=20 in panel (b) respectively correspond to rough visual averages of the centre three and all five displayed posteriors in panel (a). For h=0 we obtain lim_N→∞π̂(λ|n_obs)=π(λ|n_obs) as for standard Monte Carlo algorithms.




Similar interpretations and conclusions arise when the data y are continuous, as we examine for a different model in the following Subsection. This also allows us to introduce a fundamental concept in ABC methods – the use of summary statistics.






§ A G-AND-K DISTRIBUTION ANALYSIS





The univariate g-and-k distribution is a flexible unimodal distribution that is able to describe data with significant amounts of skewness and kurtosis.
Originally developed by 
tukey77 (see also martinez+i84,hoaglin85 and
rayner+m02), the g-and-k and related distributions have been analysed in the ABC setting by peters+s06, allingham+km09, drovandi+p11 and fearnhead+p12 among others. Its density function has no closed form, but is alternatively defined through its quantile function as

    Q(q|A,B,g,k) = A + B[1+c1-exp{-gz(q)}/1+exp{-gz(q)}] (1+z(q)^2)^k z(q)

for B>0, k>-1/2, 
where z(q)=Φ^-1(q) is the q-th quantile of the standard normal distribution function. The parameter c measures overall asymmetry, and is conventionally fixed at c=0.8  (resulting in k>-1/2) <cit.>. This distribution is very flexible, with many common distributions obtained or well approximated by particular parameter settings, such as the normal distribution when g=k=0. Given θ=(A,B,g,k)^⊤, simulations z(q)∼ N(0,1) drawn from a standard normal distribution can be transformed into samples from the g-and-k distribution through equation (<ref>).
	
Figure <ref> shows a scatterplot of samples from the likelihood-free approximation of the posterior π(θ|y_obs) (grey dots), based on a simulated dataset y_obs of length n=1,000 generated from the g-and-k distribution with parameter vector θ_0=(3,1,2,0.5)^⊤. 

This analysis was based on defining 
y-y_obs = (y-y_obs)^⊤Σ̂^-1(y-y_obs)≤ h as Mahalanobis distance, with h given by the 0.005 quantile of the differences y-y_obs for i=1,…, N=100,000 Monte Carlo samples from the joint prior π(θ)=π(A)π(B)π(g)π(k)=N(1,5)× N(0.25,2)× U(0,10)× U(0,1).

The matrix Σ̂ was determined as the sample covariance matrix of y using 2,000 samples generated under the model y|θ_0 with θ=θ_0 fixed at its true value.






As is apparent from Figure <ref>, the likelihood-free approximation to π(θ|y_obs) (grey dots) is particularly poor – the true parameter vector θ_0 is not even close to  the estimated posterior samples. This outcome is a direct result of the dimension of the comparison y-y_obs. The chance of generating an n=1,000-dimensional vector y that is close to y_obs, even if θ=θ_0, is vanishingly small. 

The odds of matching y with y_obs can be increased by redefining both in terms of their order statistics, although the chances still remain extremely low (see Example 3 in Section <ref> for an illustration).

This means that h must be relatively large, which results in accepting samples θ^(i) that generate data y^(i) that are not actually close to y_obs, and thereby producing a poor approximation to π(θ|y_obs).

The obvious way to avoid this problem is to reduce the dimension of the data comparison y-y_obs. Suppose that lower dimensional statistics s=S(y) and s_obs=S(y_obs) are available, such that S(y) is sufficient for, or highly informative for θ under the model, but where (S(y))≪(y). Then the comparison y-y_obs might be replaced by s-s_obs without too much loss of information, but with the advantage that the dimension of S(y) is now much lower. That is, step 3 in the likelihood-free rejection sampling algorithm could be further replaced by:



Using this idea, drovandi+p11 suggested the statistics

    S_A   =   E_4,
    	  
    	S_B=E_6-E_2,
    	  
    	S_g=(E_6+E_2-2E_4)/S_B,
    
    	S_k   =   (E_7-E_5+E_3-E_1)/S_B

as informative for A, B, g and k respectively, so that S(y)=(S_A,S_B,S_g,S_k)^⊤, where
E_1≤ E_2≤…≤ E_8 are the octiles of y. Repeating the above g-and-k analysis but using the 4-dimensional comparison s-s_obs rather than y-y_obs (and recomputing Σ̂ and h under the same conditions), the resulting posterior samples are shown in Figure <ref> (black dots).

The difference in the quality of the approximation to π(θ|y_obs) when using S(y) rather than y, is immediately apparent. The true parameter value θ_0 is now located firmly in the centre of each pairwise posterior sample, several parameters (particularly A and g) are more precisely estimated, and evidence of dependence between parameters (as is to be expected) is now
clearly seen. 

While it is unreasonable to expect that there has been no loss of information in moving from y to S(y), clearly the overall gain in the quality of the approximation to the likelihood-free posterior
has been worth it in this case. This suggests that the use of summary statistics S(y) is a useful tool more generally in approximate Bayesian computational techniques.







§ LIKELIHOOD-FREE METHODS OR APPROXIMATE BAYESIAN COMPUTATION (ABC)?




The terms likelihood-free methods and approximate Bayesian computation are both commonly used to describe Bayesian computational methods developed for when the likelihood function is computationally intractable, or otherwise unavailable. Of course, “likelihood-free” is arguably a misnomer – in no sense is the likelihood function not involved in the analysis. It is the function used to generate the data y∼ p(y|θ), and it accordingly must exist, whether or not it can be numerically evaluated or written down. Rather, in this context, “likelihood-free” refers to any likelihood-based analysis that proceeds without direct numerical evaluation of the likelihood function. There are several techniques that could be classified according to this description.



“Approximate Bayesian computation”, commonly abbreviated to “ABC”, was first coined by  beaumont+zb02 in the context of Bayesian statistical techniques in population genetics (although see tavare17, this volume), and refers to the specific type of likelihood-free methods considered in this book. In particular, given the “approximate” in ABC, it refers to those likelihood-free methods that produce an approximation to the posterior distribution resulting from the imperfect matching of data y-y_obs or summary statistics s-s_obs.

Thus, the likelihood-free rejection algorithm described above with h=0, which only accepts samples,  θ, which have exactly reproduced the observed data y_obs, is not an ABC algorithm, as the method produces exact samples from the posterior distribution – there is no approximation. (The Monte Carlo approximation of the posterior is not considered an approximation in this sense.) It is, however, a likelihood-free method.  Whereas, the likelihood-free rejection algorithm which may accept samples if y-y_obs≤ h, for h>0, is an ABC algorithm, as the samples will be drawn from an approximation to the posterior distribution.

Similarly, when the sampler may alternatively accept samples if s-s_obs≤ h, for any h≥ 0 (including h=0), the resulting samples are also drawn from an approximate posterior distribution. As such, this is also an ABC algorithm. The only exception to this is the case where h=0 and the summary statistics are sufficient: here there is no posterior approximation – the algorithm is then likelihood-free but not an ABC method.  

With a few exceptions (such as indirect inference, see drovandi18) all of the methods considered in this book are both ABC  and (by definition) likelihood-free methods. The aim of any ABC analysis is to find a practical way of performing the Bayesian analysis, while keeping the Approximation and the Computation to a minimum.






§ THE APPROXIMATE POSTERIOR DISTRIBUTION






In contrast to the intuitive development of likelihood-free methods in the previous Sections,
we now describe the exact form of the ABC approximation to the posterior distribution that is produced from the  likelihood-free rejection algorithm.

The procedure of (i) generating θ from the sampling distribution, g(θ),  (ii) generating data, y, from the likelihood, p(y|θ), conditional on θ, and (iii) rejecting θ if y-y_obs≤ h,
is equivalent to drawing a sample (θ,y) from the joint distribution proportional to

    I(y-y_obs≤ h)p(y|θ)g(θ),

where I is the indicator function, with  I(Z)=1 if Z is true, and I(Z)=0 otherwise.
If this sample (θ,y) is then further accepted with probability proportional to π(θ)/g(θ), this implies that the likelihood-free rejection algorithm is sampling from the joint distribution proportional to

    I(y-y_obs≤ h)p(y|θ)g(θ)π(θ)/g(θ)
    
    	=
    	I(y-y_obs≤ h)p(y|θ)π(θ).

Note that if h=0, then the θ marginal of (<ref>) equals the true posterior distribution, as

    lim_h→ 0∫ I(y-y_obs≤ h)p(y|θ)π(θ) dy
    	    =    ∫δ_y_obs(y)p(y|θ)π(θ) dy
        =    
    	p(y_obs|θ)π(θ).

That is, for h=0, the likelihood-free rejection algorithm draws samples, (θ,y), for which the marginal distribution of the parameter vector is the true posterior, π(θ|y_obs). (The marginal distribution of the auxiliary dataset y is a point mass at {y=y_obs} in this case.)


It is useful in the following to generalise the above formulation slightly. In (<ref>), the indicator  term I(y-y_obs≤ h) only takes the values 0 or 1. This is useful in the sense that it allows clear “If y-y_obs≤ h then …” statements to be made in any   algorithm, which can simplify implementation. However it is intuitively wasteful of information, as it does not discriminate between those samples, θ, for which the associated dataset y exactly equals the observed dataset y_obs, and those samples, θ, for which the associated dataset is the furthest away from y_obs, i.e. y-y_obs=h. As the former case produces samples that are exact draws from the true posterior distribution, whereas the latter case does not, this produces a motivation for a more continuous scaling from 1 (when y=y_obs) to 0 (when y-y_obs is large).


This can be achieved by replacing the indicator function, I(y-y_obs≤ h), with a standard smoothing kernel function, K_h(u), with u=y-y_obs, where

    K_h(u)=1/hK(u/h).

Kernels are symmetric functions such that K(u)≥ 0 for all u, ∫ K(u)du=1, ∫ uK(u)du=0 and ∫ u^2K(u)du<∞.
Here, h>0 corresponds to the scale parameter, or “bandwidth” of the kernel function. Several common forms for kernel functions are given in Table <ref>, and these are illustrated in Figure <ref>. Following convention, we define lim_h→ 0K_h(u) as a point mass at the origin (u=0).






An alternative specification of a smoothing kernel for multivariate datasets is obtained by writing u=y-y_obs, where u=(u_1,…,u_n)^⊤, y=(y_1,…,y_n)^⊤ and y_obs=(y_obs,1,…,y_obs,n)^⊤, so that u_i=y_i-y_obs,i.  Then we can write
K_h(u)=∏_i=1^n K_h_i(u_i), where the scale parameter of each individual kernel function, K_h_i(u_i), may vary. A further, more general specification may  determine K_h(u) as a fully multivariate, smooth and symmetric function, satisfying the above moment constraints. One such example is a multivariate N(0,Σ) distribution, for some fixed covariance matrix Σ.

Substituting
the kernel function, K_h(u), into the likelihood-free rejection algorithm results in the ABC Rejection Sampling Algorithm:




In order to determine the form of the target distribution, π_ABC(θ|y_obs), of this algorithm, we can
follow the same argument as before. By (i) generating θ from the importance distribution, g(θ), (ii) generating data, y, from the likelihood, p(y|θ), conditional on θ, and then (iii) accepting the sample (θ,y) with probability proportional to K_h(y-y_obs)π(θ^(i))/g(θ^(i)), this results in samples from the joint distribution

    π_ABC(θ,y|y_obs) ∝ K_h(y-y_obs)p(y|θ)π(θ).

When K_h(u) is the uniform kernel (see Table <ref>), then (<ref>) reduces to (<ref>).

Accordingly, we define the ABC approximation to the true posterior distribution as

    π_ABC(θ|y_obs) = ∫π_ABC(θ,y|y_obs) dy,

where π_ABC(θ,y|y_obs) is given by (<ref>).

As before, as h→ 0, so that only those samples, θ, that generate data for which y=y_obs are retained,
then (<ref>) becomes

    lim_h→ 0π_ABC(θ,y|y_obs)
    	   ∝   lim_h→ 0 K_h(y-y_obs)p(y|θ)π(θ)
        =    δ_y_obs(y)p(y|θ)π(θ),

and so 
	lim_h→ 0π_ABC(θ|y_obs)∝∫δ_y_obs(y)p(y|θ)π(θ)dy
	=p(y_obs|θ)π(θ).
That is, samples from the true posterior distribution are obtained as h→ 0. However, h=0 is not a viable choice in practice, as for continuous y_obs it corresponds to an algorithm with an acceptance rate of zero. 

To see what marginal distribution the ABC rejection algorithm is sampling from for h>0 we can integrate π_ABC(θ,y|y_obs) over the  auxiliary data margin, y.





A natural question to ask is, how accurate is this approximation? Re-writing the right hand side of (<ref>) without the prior distribution, π(θ), we can similarly define the ABC approximation to the true likelihood, p(y|θ), for a fixed value of θ, as

    p_ABC(y_obs|θ) = ∫ K_h(y-y_obs)p(y|θ) dy.

In this manner, ABC can be interpreted as a regular Bayesian analysis, but with an approximated likelihood function.


 Working in the univariate case for simplicity of illustration, so that y,y_obs∈𝒴=ℝ and u=|u|, we can obtain

    p_ABC(y_obs|θ)     =    ∫ K_h(|y-y_obs|)p(y|θ)dy
        =    ∫ K(u)p(y_obs-uh|θ)du
        =    ∫ K(u)[p(y_obs|θ) - uhp'(y_obs|θ) + u^2h^2/2p”(y_obs|θ) - …] du
        =    
    	p(y_obs|θ)+ 1/2h^2p”(y_obs|θ)∫ u^2K(u)du - …

using the substitution u=(y_obs-y)/h, a Taylor expansion of p(y_obs-uh|θ) around the point y_obs, and the kernel function properties of K_h(u)=K(u/h)/h, ∫ K(u)du=1, ∫ uK(u)du=0 and K(u)=K(-u).
The above is a standard smoothing kernel density estimation expansion, and assumes that the likelihood, p(y|θ), is infinitely differentiable. As with kernel density estimation, the choice of scale parameter is more important than the choice of kernel function in terms of the quality of the approximation.


Then,  the pointwise bias in the likelihood approximation for fixed θ can be expressed as

    b_h(y|θ) := p_ABC(y|θ)-p(y|θ),

as a function of y,
which to second order can be written as

    b̂_h(y|θ) =  1/2h^2σ^2_Kp”(y|θ),

where σ^2_K=∫ u^2K(u)du is the variance of the kernel function. Accordingly,
the magnitude of the bias is reduced if h is small, corresponding to better approximations. Clearly, the second derivative of the likelihood function, p”(y|θ), is typically also unavailable if the likelihood function itself is computationally intractable.

When y,y_obs∈𝒴 is multivariate, a similar derivation to the above is available. In either case, the ABC approximation to the true posterior is defined through (<ref>).


In a similar manner, we can determine the pointwise bias in the resulting ABC posterior approximation.
From (<ref>) we can write

    b_h(y_obs|θ)π(θ)      =     p_ABC(y_obs|θ)π(θ)-p(y_obs|θ)π(θ)
        =    π_ABC(θ|y_obs)c_ABC-π(θ|y_obs)c,

where c_ABC=∫ p_ABC(y_obs|θ)π(θ)dθ>0 and c=∫ p(y_obs|θ)  π(θ)dθ>0. Rearranging (<ref>), we obtain

    a_h(θ|y_obs)     :=    π_ABC(θ|y_obs)-π(θ|y_obs)
       =   b_h(y_obs|θ)π(θ) + π(θ|y_obs)c/c_ABC-π(θ|y_obs)
        =    b_h(y_obs|θ)π(θ)/c_ABC+π(θ|y_obs)(c/c_ABC-1),

as a function of θ. As h→ 0, then b_h(y_obs|θ)→ 0 from (<ref>), and so p_ABC(y_obs|θ)→ p(y_obs|θ) pointwise, for fixed θ. Further, c/c_ABC→ 1 as h gets small, so that a_h(θ|y_obs)→ 0.







 §.§ Simple examples




In many simple cases, the ABC approximation to the posterior distribution can be derived exactly.



Example 1:

 Suppose that the observed data, y_obs, is a single draw from a univariate density function p(y|θ), and that θ is a scalar. If we consider the particular case where K_h(u) is the uniform kernel on [-h,h] (see Table <ref>), and u=|u|, then we have
 
    π_ABC(θ|y_obs)    ∝   π(θ)∫_-∞^∞ K_h(|y-y_obs|)p(y|θ)dy
        =    π(θ)/2h∫_y_obs-h^y_obs+hp(y|θ)dy
        =    π(θ)[P(y_obs+h|θ)-P(y_obs-h|θ)]/2h,

where P(y|θ)=∫_-∞^y p(z|θ)dz is the cumulative distribution function of y|θ. Noting that as
lim_h→ 0[P(y_obs+h|θ)-P(y_obs-h|θ)]/2h = p(y_obs|θ)
via l'Hopital's rule,  then π_ABC(θ|y_obs)→π(θ|y_obs) as h→ 0, as required. Also,
[P(y_obs+h|θ)-P(y_obs-h|θ)]/2h≈ 1/2h for large h, and so π_ABC(θ|y_obs)→π(θ) as h→∞.





Suppose now that p(y|θ)=θ e^-θ y, for θ, y≥ 0, is the density function of an Exp(θ) random variable, and that the prior π(θ)∝θ^α-1e^-βθ is given by a Gamma(α, β) distribution with shape and rate parameters α>0 and β>0. Then from (<ref>), and for 0<h<y_obs+β, we can directly obtain

    p_ABC(y_obs|θ)     =    1/2he^-θ y_obs(e^θ h-e^-θ h)
    b̂_h(y_obs|θ)     =    1/6h^2θ^3e^-θ y_obs
    π_ABC(θ|y_obs)    =   θ^α-1e^-θ(y_obs+β)(e^θ h-e^-θ h)
    /Γ(α)/(y_obs+β-h)^α-Γ(α)/(y_obs+β+h)^α,

where Γ(α)=∫_0^∞ z^α-1e^-zdz is the gamma function.

Figure <ref>(a) illustrates the true likelihood function, p(y|θ), (black dashed line) and  the ABC approximation to the true likelihood function, p_ABC(y|θ), (solid grey line) as a function of y for h=0.91 and θ=2. Also shown (grey dashed line), is the second order approximation to the ABC likelihood function, p(y|θ)+b̂_h(y|θ). In this case, the second order approximation provides a reasonable representation of the ABC likelihood, p_ABC(y|θ). For other choices of h and θ, the quality of this representation will vary.



The ABC approximation, π_ABC(θ|y_obs), to the true posterior π(y_obs|θ), given y_obs=2 and α=β=1.2 is shown in Figure <ref>(b) for various values of h=0.01, …, 2.7 (grey lines). The true posterior is illustrated by the black dashed line. For small h (h=0.01), π_ABC(θ|y_obs) is indistinguishable from the true posterior. As h increases, so does the scale of the approximate posterior, which begins to exhibit a large loss of precision compared to the true posterior. Both mean and mode of π_ABC(θ|y_obs) increase with h.


Finally, Figure <ref>(c) shows the resulting 
bias, a_h(θ|y_obs), in the ABC posterior approximation as a function of θ and h. Dashed and solid lines respectively show the exact bias a_h(θ|y_obs) and the second order bias â_h(θ|y_obs) (defined as a_h(θ|y_obs) in (<ref>) but with b̂_h(y|θ) substituted for b_h(y|θ)).
 Clearly, the bias in the main body of the distribution, particularly in the region around the mode, is well described by the second order approximation, â_h(θ|y_obs), whereas the bias in the distributional tails is more heavily influenced by terms of higher order than two.






























Example 2:

Suppose that the observed data, y_obs=(y_obs,1,…,y_obs,n)^⊤, are n independent draws from a univariate N(θ,σ_0^2) distribution, where the standard deviation, σ_0>0, is known.
For this model we know that p(y_obs|θ)∝ p(y̅_obs|θ), where y̅_obs=1/n∑_i y_obs,i, as the sample mean is a sufficient statistic for θ.

If we specify K_h(u) as a Gaussian N(0,h^2) kernel (see Table <ref>), then the ABC approximation to the likelihood, p(y̅_obs|θ) is given by

 
    p_ABC(y̅_obs|θ)
    	    =    ∫_-∞^∞ K_h(|y̅-y̅_obs|)p(y̅|θ)dy̅
        =    ∫_-∞^∞1/√(2π)hexp{-(y̅-y̅_obs)^2/2h^2}√(n)/√(2π)σ_0exp{-n(y̅-θ)^2/2σ_0^2} dy̅
       ∝   exp{  -(θ-y̅_obs)^2/2(σ^2_0/n+h^2)}

for h≥0.
That is, y̅_obs∼N(θ,σ_0^2/n+h^2) under the ABC approximation to the likelihood. In comparison to the true likelihood,  for which y̅_obs∼N(θ,σ_0^2/n),
the variance is inflated by h^2, the variance of the Gaussian kernel.

Accordingly, if the prior for θ is given by a N(m_0,s_0^2) distribution, where m_0 and s_0>0 are known, then

    π_ABC(θ|y_obs)=ϕ(m_0s_0^-2+y̅_obs(σ_0^2/n+h^2)^-1/s_0^-2+(σ_0^2/n+h^2)^-1,
    1/s_0^-2+(σ_0^2/n+h^2)^-1),

where ϕ(a,b^2) denotes the density of a N(a,b^2) distributed random variable.

Clearly π_ABC(θ|y_obs)→π(θ|y_obs) as h→ 0. However, the approximation will be quite reasonable if σ^2/n is the dominating component of the variance so that h is small in comparison <cit.>.

A similar result to the above is available in the case of a multivariate parameter vector, θ.





Figure <ref>(a) illustrates the resulting ABC posterior approximation π_ABC(θ|y_obs) with y̅_obs=0 when σ^2/n=1 for the improper prior given by m_0=0, s^2_0→∞, so that the true posterior distribution is N(0,1) (dashed line). The approximation is clearly quite reasonable for h=0.5 and h=0.1 as then h^2<σ^2_0/n.

Figure <ref>(b) shows the same posterior approximations  but based on a uniform kernel over [-h, h] for K_h(u), rather than the Gaussian N(0,h^2) kernel. This ABC posterior is derived from (<ref>). The resulting forms for π_ABC(θ|y_obs) are no longer within the Gaussian family for h>0, exhibit a flatter behaviour around the mean, and are more concentrated around the mean due to the compact support of the uniform kernel. The approximations with either kernel perform well for small h.







This example additionally provides some insight into the asymptotic behaviour of the ABC posterior approximation. Following standard likelihood asymptotic results, when the amount of data, n, becomes large, the true likelihood function, p(y|θ), will approximately behave as a Gaussian distribution. As most prior distributions will have little impact in this setting (they will be approximately constant over the region of high posterior density), it follows that the ABC posterior approximation, π_ABC(θ|y_obs) will follow a Gaussian distribution with a variance that is inflated by an h^2 term. Consequently, the ABC posterior approximation, π_ABC(θ|y_obs) may then in principle be improved simply by rescaling the posterior variance to remove this term <cit.>.










§ THE USE OF SUMMARY STATISTICS








 §.§ Summary statistic basics






Despite the development in the previous Section, the ABC posterior approximation π_ABC(θ|y_obs)∝∫ K_h(y-y_obs)p(y|θ)p(θ)dy is rarely used in practice. This is because, except in very specific scenarios (such as when y_obs is very low dimensional, or when the likelihood function p(y|θ) factorises into very low dimensional components), it is highly unlikely that y≈ y_obs can be generated from p(y|θ) for any choice of θ for realistic datasets. This results in the need to use a large value of the kernel scale parameter h in order to achieve viable rejection sampling algorithm acceptance rates (or a similar loss of performance in other algorithms), and in doing so produce poorer ABC posterior approximations.



In the stereological extremes analysis in Section <ref> we replaced the full dataset y_obs with a sufficient statistic n_obs for the model parameter λ when estimating π(θ|y_obs)=π(λ|n_obs). As sufficient statistics can be much lower dimensional than the full dataset, it is clear that greater approximation accuracy can be achieved for the same computational overheads when using low dimensional  statistics  (which is hinted at in the g-and-k distribution analysis in Section <ref>).

The following example, based on drovandi12, highlights the computational benefits in using lower dimensional, and less variable sufficient statistics.




Example 3:

Suppose that y=(y_1,y_2)^⊤,  where y_i∼Binomial(n,θ) with θ∼ U(0,1). Consider three possible vectors of sufficient statistics: s^1=(y_1,y_2)^⊤ is the full dataset, s^2=(y_(1),y_(2))^⊤ are the order statistics y_(1)≤ y_(2), and s^3=y_1+y_2 is the sum of the two individual values. All three vectors of statistics are sufficient for this simple model.

It is easy to compute the marginal distribution of each summary statistic p_i(s^i) = ∫_0^1p(s^i|θ)π(θ)dθ as follows:

    p_1(s^1)     =    ∫_0^1∏_i=1^2([   n; y_i ])
    				θ^y_i(1-θ)^n-y_i dθ
        =    ([   n; y_1 ])([   n; y_2 ])
    				B(y_1+y_2+1,2n-y_1-y_2+1),
    
    	p_2(s^2)     =    [2-I(y_(1)=y_(2))]∫_0^1∏_i=1^2([   n; y_i ])
    				θ^y_i(1-θ)^n-y_i dθ
        =    [2-I(y_(1)=y_(2))]
    			 ([   n; y_1 ])([   n; y_2 ])
    				B(y_1+y_2+1,2n-y_1-y_2+1),
    
    	p_3(s^3)     =    ∫_0^1([  2n; s^3 ])
    				θ^s^3(1-θ)^2n-s^3 dθ
        =    ([  2n; s^3 ])
    				B(s^3+1,2n-s^3+1)
        =     1/(2n+1),

where (a,b)=∫_0^1z^a-1(1-z)^b-1dz is the beta function. Here, p_i(s^i) is the probability of generating the vector s^i under an ABC rejection sampling algorithm with sampling distribution given by the prior, g(θ)=π(θ). That is, p_i(s_i) is the acceptance probability of the algorithm if we only accept those sufficient statistics that exactly match the observed sufficient statistics.

Suppose that we observe y_obs=(y_obs,1,y_obs,2)^⊤=(1,2)^⊤ from n=5 experiments. From the above, we have algorithm acceptance rates of:

    p_1(s^1_obs)=5/132≈ 0.038,      p_2(s^2_obs)=5/66≈ 0.076       p_3(s^3_obs)=1/11≈ 0.091,

where s^i_obs denotes the statistic s^i derived from y_obs.
The probability p_1(s^1)  is the probability of generating first y_1=1  and then y_2=2. As a result, p_1(s^1) will decrease rapidly as the length of the observed dataset y_obs increases. The probability p_2(s^2) corresponds to the probability of generating either y=(1,2)^⊤ or y=(2,1)^⊤, which are equivalent under the binomial model. Hence, s^2 has twice the probability of s^1 of occurring. Finally, the probability p_3(s^3), is the probability of generating y=(1,2)^⊤, (2,1)^⊤, (0,3)^⊤ or (3,0)^⊤. Each of these cases are indistinguishable under the assumed model, and so the event s^3 occurs with the largest probability of all.

Quite clearly, while still producing samples from the true target distribution, π(θ|y_obs), the impact on the efficiency of the sampler of the choice of sufficient statistics is considerable, even for an analysis with only two observations, y_1 and y_2.  The most efficient choice is the minimal sufficient statistic. The differences in the acceptance rates of the samplers would become even greater for larger numbers of observations, n.







While the optimally informative choice of statistic for an ABC analysis is a minimal sufficient statistic, this may still be non-viable in practice. For example, if the minimal sufficient statistic is the full dataset y_obs, sampling from π_ABC(θ|y_obs) will be highly inefficient even for moderately sized datasets. Similarly, in a scenario where the likelihood function may not be known beyond a data generation procedure, identification of any low-dimensional sufficient statistics (beyond, trivially, the full dataset y_obs) may be impossible. Further, low dimensional sufficient statistics may not even exist, depending on the model.

In general, a typical ABC analysis will involve specification of a vector of summary statistics s=S(y), where (s)≪(y). The rejection sampling algorithm with then contrast s with s_obs=S(y_obs), rather than y with y_obs. As a result, this procedure will produce samples from the distribution π_ABC(θ|s_obs) as follows:




Similar to the discussion in Section <ref>, it can be seen that the ABC posterior approximation now has the form

    π_ABC(θ|s_obs) ∝∫ K_h(s-s_obs)p(s|θ)π(θ)ds,

where p(s|θ) denotes the likelihood function of the summary statistic s=S(y) implied by p(y|θ).

(That is, p(s|θ) = ∫_𝒴δ_s(S(y))p(y|θ) dy.)

If we let h→ 0, so that only those samples, θ, that generate data for which s=s_obs are retained, then

    lim_h→ 0π_ABC(θ|s_obs)    ∝   ∫lim_h→ 0K_h(s-s_obs)p(s|θ)π(θ)ds
        =    ∫δ_s_obs(s)p(s|θ)π(θ)ds
        =     p(θ|s_obs)π(θ).

Hence, samples from the distribution π(θ|s_obs) are obtained as h→ 0. If the vector of summary statistics, s=S(y), is sufficient for the model parameters, then π(θ|s_obs)≡π(θ|y_obs), and so samples are produced from the true posterior distribution.  However, if S(y) is not sufficient – and this is typically the case in practice – then the ABC posterior approximation is given by (<ref>), where in the best scenario (i.e. as h→ 0) the approximation is given by π(θ|s_obs).




The following example illustrates the effect of using a non-sufficient summary statistic.



Example 4:

Consider again the univariate Gaussian model in Example 2. Suppose that we modify  this example <cit.>, so that the model still assumes that the observed data y_obs=(y_obs,1,…,y_obs,n)^⊤ are random draws from a univariate N(θ,σ^2_0) distribution, but where we now specify an insufficient summary statistic, s=y̅_1:n'=1/n'∑_i=1^n'y_i with n'< n.


Writing s_obs=S(y_obs), the resulting ABC approximation to the likelihood function becomes
 
    p_ABC(s_obs|θ)
    	    =    ∫ K_h(s-s_obs)p(s|θ)ds
       ∝   ∫_-∞^∞1/√(2π)hexp{-(s-s_obs)^2/2h^2}√(n')/√(2π)σ_0exp{-n'(s-θ)^2/2σ^2_0}ds
       ∝   exp{  -(θ-s_obs)^2/2(σ^2_0/ω n+h^2)},

where ω=n'/n is the proportion of the n observations used in the vector of summary statistics.

That is, s_obs∼ N(θ,σ^2/ω n + h^2).  When ω=1, then s_obs=y̅_obs is sufficient for θ and so s_obs∼ N(θ,σ^2/n + h^2) recovers the same result as Example 2.

When n'<n, so that s is no longer sufficient for θ, 
the mean of the Gaussian likelihood function is centred on the mean y̅_obs,1:n' rather than y̅_obs,1:n, but more critically
the variance of the Gaussian likelihood is σ^2/ω n + h^2. It is evident that there are now two sources of error, both of which inflate the variance of the likelihood. The first, h^2, arises through the matching of the simulated and observed data through the Gaussian kernel. The second source of error comes from the 0<ω<1 term, which can be interpreted as the degree of inefficiency of replacing y by s=S(y). That is, the use of non-sufficient statistics reduces the precision of the likelihood (and by turn, the posterior distribution) in this case.


From Example 2, it follows that when n is large and the posterior is asymptotically Gaussian, the ABC posterior approximation, π_ABC(θ|s_obs), can be improved by rescaling to remove h^2 from the posterior variance. However, correcting for the lack of sufficiency in the summary statistic, s, would require knowledge of the relative inefficiency of s over y, which may be difficult to obtain in practice.







The choice of summary statistics for an ABC analysis is a critical decision that directly affects the quality of the posterior approximation. Many approaches for determining these statistics are available, and these are reviewed in blum+nps13 and prangle17, this volume. These methods seek to trade off two aspects of the ABC posterior approximation that directly result from the choice of summary statistics. The first is that π(θ|y_obs) is approximated by π(θ|s_obs). As this represents an irrevocable potential information loss, the information content in s_obs should be high. The second aspect of the ABC posterior approximation is that the simulated and observed summary statistics are compared within a smoothing kernel K_h(s-s_obs) as part of the form of π_ABC(θ|s_obs) (<ref>). As stochastically matching s and s_obs becomes increasingly difficult as the dimension of the summary statistics increases, the dimension of s should be low.

As such, the dimension of the summary statistic should be large enough so that it contains as much information about the observed data as possible, but also low enough so that the curse-of-dimensionality of matching s and s_obs is avoided. For illustration, in Example 3, the optimum choice of summary statistic is a minimal sufficient statistic. However, for other models it may be the case that the dimension of the minimal sufficient statistic is equal to that of the original dataset. As this will cause curse-of-dimensionality problems in matching s with s_obs, it is likely that a more accurate ABC posterior approximation can be achieved by using a lower-dimensional non-sufficient statistic, rather than remaining within the class of sufficient statistics. This was indeed the case in the g-and-k distribution analysis in Section <ref>. 








 §.§ Some practical issues with summary statistics





Even with the above principles in mind, summary statistic choice remains one of the most challenging aspects of implementing ABC in practice. For instance, it is not always viable to continue to add summary statistics to s until the resulting ABC posterior approximation does not change for the worse, as is illustrated by the following example.



Example 5:

Suppose that  y=(y_1,…,y_n)^⊤ with y_i∼Poisson(λ). Combined with conjugate prior beliefs λ∼Gamma(α,β) this gives λ|y∼Gamma(α+ny̅,β+n). For this model we know that the sample mean y̅ is a sufficient statistic. However, we also know that the mean and variance of a Poisson(λ) model are both equal to λ, and so we might also expect the sample variance v^2 to also be informative for λ, although it is not sufficient. Suppose that we observe y_obs=(0,0,0,0,5)^⊤ which gives (y̅_obs,v^2_obs)=(1,5). Here, as the sample mean and variance are quite different from each other, we might expect that the Poisson model is not appropriate for these data.





Figure <ref> illustrates various ABC posterior approximations to the true target distribution (solid lines)
based on a prior with α=β=1 (dashed lines), with
K_h(u) specified as a uniform kernel over [-h,h] and u representing Euclidean distance.

The top row illustrates the resulting posterior approximations, π(λ|s_obs), when the summary statistics s are given as the sample mean y̅ (left panel), the sample standard deviation v (centre), or both (right) when the kernel scale parameter  is h=0. Using s=y̅ recovers the true posterior exactly, which is no surprise as y̅ is a sufficient statistic. Using s=v produces an informed ABC approximation, but one which is based on a variance that is consistent with a larger mean under the Poisson model. When s=(y̅,v)^⊤ then we again obtain the true posterior distribution as π(λ|y̅_obs,v_obs)≡π(λ|y̅_obs) through sufficiency, and the additional information that v brings about the sample y has no effect on the ABC estimated posterior.

The bottom row in Figure <ref> shows the same information as the top row, except that the kernel scale parameter is now non-zero (h=0.3). The posterior approximations based on s=y̅ and s=v are minor deviations away from those in the top row when h=0. This occurs as the values of λ that are able to reproduce the observed summary statistics within a non-zero tolerance h=0.3 are slightly different to those that can reproduce the summary statistics exactly. However, the third panel with s=(y̅,v)^⊤ is clearly biased to the right, with the resulting ABC posterior approximation visually appearing to be a loose average of those
distributions with s=y̅ and s=v.


This behaviour is different from when h=0. In that case, when adding more information in the vector of summary statistics in going from s=y̅ to s=(y̅,v)^⊤, the posterior approximation does not change as the summary statistic s=y̅ is sufficient and it is being matched exactly. However, when h>0, because the ABC algorithm allows a non perfect matching of the sufficient statistic y̅, it additionally allows the extra information in the sample standard deviation v to also contribute to the approximation. In this case, because the observed summary statistics y̅_obs and v_obs are inconsistent with respect to the model, this then results in a strongly biased fit when moving from s=y̅ to s=(y̅,v)^⊤. 

As such, while it may be tempting to include progressively more summary statistics into s_obs until the ABC posterior approximation does not change appreciably, the assumption that that this will provide the most accurate posterior approximation is clearly incorrect. Even if s_obs contains sufficient statistics for the model, the inclusion of further statistics can still bias the posterior approximation, particularly in the case where the observed data are inconsistent with the model.




The identification of suitable summary statistics is clearly a critical part of any analysis. Accordingly many techniques have been developed for this purpose – see e.g. blum+nps13 and prangle17 (this volume) for a detailed review and comparison of these methods. While the choice of summary statistics is itself of primary importance, it is less appreciated that the distance measure · can also have a substantial impact on ABC algorithm efficiency, and therefore the quality of the posterior approximation.


Consider the distance measure  s-s_obs=(s-s_obs)^⊤Σ^-1(s-s_obs). Here we can specify the covariance matrix Σ as the identity matrix to produce Euclidean distance, or as a diagonal matrix of non-zero weights to give weighted Euclidean distance (e.g. hamilton+crhbe05,luciani+sjft09) or as a full covariance matrix to produce Mahalanobis distance (e.g peters+fs12,erhardt+s16).

To see why standard and weighted Euclidean distance can be a poor choice, consider the setting in Figure <ref>, where candidate parameter values, θ, generating continuous bivariate statistics, s|θ, s=(s_1,s_2)^⊤, are accepted as draws from π_ABC(θ|s_obs) if s lies within a ball of radius h, centered on s_obs.   That is, K_h is the uniform kernel on [-h,h], and · denotes Euclidean distance. 



If we reasonably suppose that the elements of s may be dependent and on different scales, their true distribution under the model may be better represented by an ellipse (grey lines). As such, an efficient ABC algorithm should accept candidate draws from π_ABC(θ|s_obs) if s|θ lies within this ellipse.

Consequently, implementing a circular acceptance region (implying independence and identical scales) induces both type I (i.e. candidate samples are rejected when they should be accepted) and type II (i.e.  candidate samples are accepted when they should be rejected) errors. 


Work linking the ABC posterior with non-parametric density estimation methods (blum10; see Section <ref>) provides support for this argument.

Here, for a multivariate kernel K_H(u)=(H)^-1K(H^-1u), where K is a symmetric multivariate density function with zero mean and finite variance, a general rule of thumb is to specify the bandwidth matrix as
H∝Σ^-1/2 where Σ is the covariance matrix of the data (e.g. scott92,wand+j95). In the ABC context, this is equivalent  to defining · as Mahalanobis distance where Σ is the covariance matrix of s (or s|θ).

Note that the above argument assumes that the summaries s_1 and s_2 are both informative for the model parameter θ. For example, in the case where s_1+s_2 is uninformative, but s_1-s_2 is informative, then it is credible that the circular acceptance region could result in a more accurate ABC posterior approximation than that resulting from the elliptical region. In general, the best acceptance region is tied up with the choice of the summary statistics in a more complicated way than that presented here (see e.g. prangle17b for a discussion).



The following example illustrates the effect that different covariance matrices Σ can have on the ABC posterior approximation.




Example 6:

Suppose that the model is specified as y_1,…,y_50∼ N(θ,1), with a uniform prior θ∼ U(-5,5). Various sufficient statistics are available for this model. We consider two alternatives: s^1=(y̅_1:40,y̅_41:50)^⊤ and s^2=(y̅_1:25-y̅_26:50, y̅_26:50)^⊤ where y̅_a:b=(b-a+1)^-1∑_i=a^by_i. In each case, given the observed sufficient statistics s_obs=(0,0)^⊤, the exact posterior distribution π(θ|y_obs) is N(0,1/50) truncated to (-5,5). However, the covariance matrices of s^1 and s^2 for fixed θ are quite different
(though they do not depend on the exact value of θ),
namely

    (s^1|θ)=([ 1/40    0;    0 1/10 ]),
    	  (s^2|θ)=([ -2/25 -1/25; -1/25 -1/25 ]),

with a negative correlation between the elements of s^2 of -1/√(2)≈-0.71.

We implement ABC using the distance measure as s-s_obs=(s-s_obs)Σ^-1(s-s_obs)' and consider the impact of the choice of Σ.





We use a version of the ABC rejection sampling algorithm (see box) that maintains a sample θ^(1), …, θ^(N) of size N from the ABC posterior approximation, which progressively lowers the kernel scale parameter h until a stopping rule is satisfied. On algorithm termination, the samples are identical to those samples that would have been obtained under the standard ABC rejection sampling algorithm if it was implemented with the lowest value of h achieved under the stopping rule. This allows us to implement a rejection sampling algorithm that will terminate when a pre-specified degree of accuracy has been achieved. The (random) number of iterations obtained before algorithm termination will accordingly be an indicator of the efficiency of the model specification – in this case, the effect of different covariance matrices Σ.



Table <ref> displays the average number of data generation steps (i.e. generating y∼ p(y|θ)) in each algorithm implementation, per final accepted particle, as a function of smoothing kernel type and the form of Σ, based on 100 replicate simulations of N=500 samples. The stopping rule continued algorithm execution until an estimate of the absolute difference between empirical (F_N(θ)) and true (F(θ)) model cumulative distribution functions was below a given level. Specifically
when ∑_i=1^N|F_N(θ^(i))-F(θ^(i))|<0.01825.

In Table <ref>, the true form of Σ is given by (s^1|θ) and (s^2|θ) (<ref>), and the diagonal form refers to the matrix constructed from the diagonal elements of (s^2|θ).



The summary statistics for s=s^1 are independent, but are on different scales. Accordingly, when this difference of scale is accounted for (Σ = true), algorithm efficiency, and therefore ABC posterior approximation accuracy, is greatly improved compared to when the difference in scale is ignored (Σ = identity). The summary statistics s^2 are both negatively correlated and on different scales. As for s^1, when summary statistic scale is taken into consideration (Σ = diagonal) an improvement in algorithm efficiency and ABC posterior approximation accuracy is achieved compared to when it is ignored. However in this case, further improvements are made when the correlation between the summary statistics is also accounted for (Σ= true). These results are consistent regardless of the form of the smoothing kernel K_h. Note that the uniform kernel produces the most efficient algorithm and most accurate ABC posterior approximation, and that this steadily worsens as the form of the kernel deviates away from the uniform density, with the worst performance is obtained under the Gaussian kernel.

This approach has been implemented in practice by e.g. luciani+sjft09 and erhardt+s16, who identify some value of θ=θ^* in a high posterior density region via a pilot analysis, and then estimate (s|θ^*) based on repeated draws from p(s|θ^*).













§ AN ABC ANALYSIS IN POPULATION GENETICS




To illustrate some of the points concerning summary statistics we consider here a population genetic example, very similar to that considered in the paper by pritch99, a key paper in the development of ABC methods. In population genetics we are often confronted with sequence data (as illustrated in Table <ref>), and we wish to infer demographic parameters that may be associated with such data. The standard modelling framework that is used is Kingman's coalescent hein04, which describes the genealogical relationship of DNA sequences in a sample. The general likelihood problem that we wish to solve then can be represented as

    p(y_obs | ϕ) = ∫_H p(y_obs | H)p(H | ϕ) dH

where y_obs represents the observed set of sequences in a sample, ϕ is an unobserved vector of parameters, and H represents the unobserved genealogy history, including mutations. A common mutation model, used here, is the infinite-sites model, in which every mutation that occurs in a genealogy is unique. Typically H is high dimensional, represented as a variable-length vector of times of events in the genealogical history, and the types of events. Although the likelihood can be computed exactly for simple demographic models and small data sets hein04 it is generally more flexible to resort to Monte Carlo methods <cit.>. 


One approach is through importance sampling. Here, an instrumental distribution q_ϕ,y(H) is available that describes the distribution of all genealogical histories H that are consistent with the data y, as a function of the model parameters ϕ. The distribution q_ϕ,y(H) is easy to simulate from and has a known functional form that can be directly evaluated. It also has the property that p(y|H')=1 for H'∼ q_ϕ,y(H).

Hence, p(y_obs|ϕ) can be estimated by

    p̂(y_obs | ϕ) = 1/N∑_i=0^Np(H^(i)|ϕ)/q_ϕ,y_obs(H^(i))

where H^(i)∼ q_ϕ,y_obs(H) for i=1,…,N.




In this analysis we compare an ABC approach to the above importance sampling method that targets the true likelihood. The aim is to investigate the performance of different summary statistics on ABC inferences, using the importance sampling-based inferences as a (noisy) ground-truth. The demographic model that generates the data is one of smooth exponential expansion. In this model the current population size N_0 contracts backwards in time as N_0(t) = N_0exp(-β t) where time t is expressed in units of 2N_0 and β = 2N_0b is the growth rate in this scaled time. An additional parameter in the model is the scaled mutation rate θ_0 = 4N_0μ. 







In the ABC analysis, simulations are carried out using the ms program of  hudson02. A technical complication that needs to be accounted for when using ms is that time in this program is scaled in units of 4N_0 rather than 2N_0 that appears standardly in most treatments (e.g. hein04), and, more importantly, in the Genetree importance sampling program  <cit.> that is used for the ground-truth. The data in Table <ref> were generated using the ms command:

which simulates one instance of 20 sequences with θ=50 and α = 30, where α = β/2 (because of the different scaling of time, noted above). Assuming independent uniform priors U(0,200) for each parameter ϕ=(θ_0,α)^⊤, it is straightforward to generate particles by sampling parameter values from the prior and then compute an importance weight for each particle using an algorithm suggested by stephens00. The implementation here (described in maciuca12) is a modification of the Genetree program 
to include the Stephens and Donnelly algorithm, following deiorio04. Although the particles could be used directly for weighted density estimation, it is computationally easier to first resample them in proportion to their weights w^(i), because the distribution of weights is typically very skewed (they have high variability). For the data in Table <ref>, N=10^8 generated particles yielded an effective sample size
 (estimated by (∑_i w^(i))^2/∑_i w^(i)2) of around 300. The following analyses are based on resampling 1000 particles.

For the ABC analysis, parameter values ϕ=(θ_0,α)^⊤ are simulated from the prior, data sets are simulated using ms, and summary statistics computed. The four summary statistics examined comprise the number of segregating sites, S_0, which corresponds to the number of mutations in the genealogy under the infinite sites mutation model, the average pairwise Hamming distance between all pairs of sequences in the sample, π_0, Tajima's D, 
and Fay and Wu's H_0. 
These latter two statistics express the difference in estimates of the scaled mutation parameter θ_0, assuming a standard coalescent model (i.e. with no population growth), based on two different unbiased estimators, one of which is π_0. The average pairwise distance, π_0, is directly an estimate of θ_0 because in the standard constant size model the expected time to coalescence for a pair of sequences is 2N_0, and therefore the expected number of mutations occurring down both branches since the common ancestor is (2N_0 + 2N_0)μ. Other estimators have been developed, based on the number of segregating sites (Watterson's estimator, used in Tajima's D), or the number of segregating sites weighted by the number of times the mutant type occurs in the sample (Fu's estimator, used in Fay and Wu's H_0). Only under the standard constant size model will these estimators all have the same expectation, and therefore deviations between them can be used to identify departures from this model. Negative values of D and positive values of H_0 are expected to be found in growing populations. The output of the ms program can be piped to a program , included with ms, which computes these four summary statistics.
The observed summary statistics are:

    s_obs=(π_0,S_0,D,H_0)^⊤= (5.90, 42,-1.64,3.67)^⊤.




ABC methods were implemented by first simulating N=1,000,000 parameter values from the U(0,200) prior distributions, storing these in the file 
(in the order indicated by the key-word )  and then running the ms program with the command

The summary statistics corresponding to these simulated data were then obtained and then s-s_obs computed as Euclidean distance. The ABC posterior approximation was obtained by using a uniform kernel K_h over [-h,h] and determining the kernel scale parameter h as the value retaining the 1000 samples for which s^(i) is closest to s_obs.



The summary statistics are measured on different scales. A common practice is to centre and scale them using the standard deviation for each summary statistic sampled from the prior predictive distribution.
(However, some authors argue that the motivations for this are flawed as an arbitrary change in the prior can change the scaling of a summary statistic within the analysis. Instead, following a similar discussion to that in Example 6, the scaling should be based on (s|θ^*) for some value of θ=θ^* in the high posterior density region, rather than (s). See e.g. erhardt+s16.)
For the present analysis, the prior predictive sample standard deviations  for π_0, S_0, D and H_0 are 14.3, 69.0, 0.50 and 7.3 respectively.
In Figure <ref>  the estimated posterior distributions using both scaled and unscaled summary statistics are shown. 






Figure <ref> compares the resulting ABC posterior approximation using (a) all four summary statistics, (b) D and H_0 only, (c) π_0 and S_0 only, or (d) π_0 or S_0 alone.

The first point to note is that the data, although quite informative about θ_0 or α jointly, do not allow us to make very detailed inference about either parameter individually i.e. they are only partially identifiable in the model – at least for these data. This is the case both for the full-likelihood and ABC inferences, although the density for the full-likelihood method, as estimated by importance sampling, tends to be more localised towards the true parameter value (indicated by a +). 

When all four summary statistics are used (panel a) the 95% HPD envelope for ABC is quite similar to that for importance sampling (black line), but is shifted towards higher values of α and θ_0. Scaled or unscaled summary statistics give similar results. The ABC posterior approximation for π_0 and S_0 together (panel b) is very similar to that for the full set of summary statistics. In this case the distances for scaled and unscaled summaries are the same because S is discrete and matched exactly. This outcome perhaps indicates that one should be cautious of adding summaries such as Tajima's D because it is simply a nonlinear function of π_0 and S_0. Whereas H_0 includes additional information from the site frequency spectrum, and would be expected to be informative (positive H_0 indicates a deficit of high-frequency derived mutations compared with that expected under the standard model). Using D and H_0 together (panel c) yields a less concentrated posterior approximation. Both statistics are based on the difference of two estimators of mutation rate, and therefore it is unsurprising that θ_0 is not well localised. The posteriors based on π_0 and S_0 individually (panel d) superficially look surprisingly similar to the full-likelihood posterior. However there is much stronger support for larger values of θ_0 and α than in the importance-sampling based posterior. 





We conduct a similar analysis  with sequence data published in hammer10  from locus 9pMB8 surveyed in 11 Biaka pygmies (resulting in 22 sequences). The data are shown in Table <ref>.  Like the simulated data above, there are 42 sites that are segregating within the Biaka sample and which are compatible with the infinite sites model. The ABC simulations were performed as previously, using all four summary statistics. The observed summary statistics for these data are

    s_obs=(π_0,S_0,D,H_0)^⊤= (7.52,42,-1.35,4.0)^⊤.

The posterior computed using importance sampling was also computed as before, but required 12 × 10^8 particles to achieve a similar effective sample size to that for the previous data set. 



It is immediately apparent from Figure <ref> that the ABC posterior approximation and ground-truth posterior are very similar, unlike the previous analysis. This differing behaviour is not due to Monte Carlo error.
The result illustrates a point that outside the exponential family there is no single, low-dimensional set of summary statistics s that will be highly informative for θ, for all observed datasets. Summary statistics that work well for one dataset may perform less well on another. In the case of the two datasets considered here, it may be argued that in the latter, despite the smaller sample size, there is a stronger signal of growth in these data, which is more readily captured by the summary statistics. For the simulated data the signal is less strong, and information in other summary statistics, such as the site frequency spectrum, or higher moments of the distribution of pairwise Hamming distances, may be required for the ABC posterior to better match the true posterior.


From a computational perspective, the 10^6 ABC simulations took about 3 minutes on a desktop computer, whereas 10^8 importance sampling simulations took around 4 hours i.e. the computational effort per iteration is broadly similar for both approaches. The algorithms used in each are `similar yet different', in that they both generate genealogical trees, but in one case the tree is constrained by the data, and in the other it is independent of the data. Naively one might think that an importance sampling  algorithm should be more efficient because it always generates a tree that is compatible with the data. However, it is typically very difficult to devise an algorithm that samples trees in proportion to their conditional distribution under the model, and therefore genealogical importance sampling tends to be inefficient, as illustrated here, where 10^8 simulations only give an effective sample size of around 300. Of course, it is possible to use sequential methods, or a pseudo-marginal method to improve efficiency (andrieu+lv17,cornuet12,beaumont03), but similar approaches are available for ABC as well.





§ LEVELS OF APPROXIMATION IN ABC





The primary challenge in implementing an ABC analysis is to reduce the impact of the approximation, while restricting the required computation to acceptable levels. In effect this is the usual “more computation for more accuracy” tradeoff. It is therefore worthwhile to briefly summarise  the quality and nature of the approximations involved in any ABC analysis.
While some of these approximations are common with standard Bayesian analyses, in particular points <ref> and  <ref> below, within the ABC framework these have additional, more subtle implications.
In order, from model conception to implementation of the analysis, the ABC approximations are:


  *  All models are approximations to the real data-generation process.

While this is true for any statistical analysis, this approximation can produce an ABC-specific issue if the assumed model is not sufficiently flexible to be able to reproduce the observed summary statistics.
In this scenario
the kernel scale parameter h will necessarily be large (as all simulated data are far from the observed data), and as a consequence the quality of the ABC approximation may be low.  Further, if, for this inflexible model, the observed summary statistics contain conflicting information for a model parameter, this may cause additional bias in the posterior approximation for this parameter, as is illustrated in Example 5. In summary, this means that the more unlikely a model is to have generated the observed data, the worse the ABC approximation will be.
In general this is problematic, as it implies that routine inspection of the fitted ABC posterior may not in itself be enough to determine model adequacy, as the ABC posterior may be a poor estimate of the true posterior, and poor data generation models may appear more likely (with h>0) than they actually are (with h=0).

By extension, this also implies that posterior model probabilities of inadequate models (constructed from the normalising constant of the poorly estimated ABC posterior distribution)  may also be affected, although this has yet to be fully explored in the literature. See fearnhead18, for an exploration of related ABC asymptotics results to date, and marin+per18 for particular methods for performing ABC model choice.


  *  Use of summary statistics rather than full datasets.

The full posterior distribution π(θ|y_obs)∝ p(y_obs|θ)π(θ) is replaced by the partial posterior π(θ|s_obs)∝ p(s_obs|θ)π(θ) where s_obs=S(y_obs) is a vector of summary statistics. If S is sufficient for θ, then there is no approximation at this stage. More commonly, for non-sufficient S, there is a loss of information.


  *  Weighting of summary statistics within a region of the observed summary statistics.


The partial posterior π(θ|s_obs) is replaced by the ABC approximation to the partial posterior

    π_ABC(θ|s_obs)∝π(θ)∫ K_h(s-s_obs)p(s|θ)ds

where K_h is a standard smoothing kernel with scale parameter h≥ 0. If h=0 or in the limit as h→ 0 then there is no further approximation at this stage. In most cases however, h>0 and so ABC makes use of a kernel density estimate as an approximation to the true likelihood function. This aspect of approximation can be a particular problem in ABC when the number of model parameters θ is large, as then the vector of summary statistics, s, must be equivalently large for parameter identifiability, and hence the comparison s-s_obs will suffer from the curse of dimensionality.


  *  Approximations due to other ABC techniques.


There are a number of other ABC techniques not discussed in this Chapter that are optionally implemented in ABC analyses in order to improve some aspect of the approximations in points 1 and 2, or to achieve a greater computational performance. Many of these are discussed in later Chapters, but some common methods involve  post-processing techniques such as regression and marginal adjustments (e.g. beaumont+zb02,blum+f10,blum+nps13,blum17,nott+ofs17),
or develop alternative approximations to the intractable likelihood function, while remaining in the ABC framework, such as Expectation-Propagation ABC, synthetic likelihoods, and copula or regression-density estimation models (e.g. barthelme+c14,barthelme+cc17,wood10,price+dln17,drovandi+mr17,li+nfs15,fan+ns13,nott+ofs17).


  *  Monte Carlo error.

In common with most Bayesian analyses, performing integrations using Monte Carlo methods introduces Monte Carlo error. Typically this error may be reduced by using larger numbers of samples from the posterior, or by reducing the variability of importance weights. The same is true for an ABC analysis, although with the additional point that more posterior samples effectively allows for a lower kernel scale parameter h and consequently an improved ABC posterior approximation.

As a result,
for a fixed number of Monte Carlo samples, the choice of kernel scale parameter represents a typical bias-variance tradeoff: if h is large, more posterior draws are available, reducing variance, but at the cost of a poorer ABC approximation; if h is small, the ABC posterior approximation is improved, but Monte Carlo variance is increased.











§ INTERPRETATIONS OF ABC





There are a number of closely related ways in which ABC methods may be understood or interpreted.

The most common of these is conditional density estimation of the posterior (e.g. blum10,bonassi+yw11,nott+ofs17) in the sense usually understood in a conventional Bayesian analysis. Before observing the data, the  distribution π(θ,y)=p(y|θ)π(θ) describes prior beliefs about the model parameters and credible datasets under the model. When a dataset y_obs is observed, interest is then in the conditional distribution of θ given that y=y_obs. In the ABC setting, π(θ,y) is represented by the joint sample (θ^(i),y^(i))∼π(θ,y), i=1,…, N. Weighting the vectors θ^(i) based on the value of y^(i)-y_obs (larger weights for smaller y^(i)-y_obs), then produces an empirical  conditional density estimate of π(θ|y_obs).



Similarly, we have already discussed that the ABC approximation to the true likelihood, p_ABC(y_obs|θ), is a kernel density estimate of p(y|θ), following (<ref>) and (<ref>). This allows ABC to be considered as a regular Bayesian analysis with an approximated likelihood function.


fearnhead+p12 noted that the ABC approximation to the posterior can be considered as a continuous mixture of posterior distributions

    π_ABC(θ|y_obs)    ∝   ∫ K_h(y-y_obs)p(y|θ)π(θ) dy
        =    ∫ w(y)π(θ|y) dy

where π(θ|y)=p(y|θ)π(θ)/π(y), with weight function w(y)∝ K_h(y-y_obs)π(y). This is the continuous equivalent of equation (<ref>) obtained during the analysis of stereological extremes in Section <ref>.


While ABC is most often thought of as an approximate method, wilkinson13 pointed out that ABC methods can be considered as exact if e=y-y_obs (or e=y-y_obs) is considered as the error (either from observation error or model misspecification) obtained in fitting the model p(y|θ) to the observed data y_obs. From this perspective, the smoothing kernel K_h is simply the density function of this error, so that e∼ K_h, and h is a scale parameter to be estimated.



Finally, while ABC methods are universally used for the analysis of models with computationally intractable likelihood functions, it is often overlooked that they also provide a useful inferential mechanism for tractable models. 

As an illustration, consider a scenario where a standard Bayesian analysis is available for a complex, but incorrect model, given the observed dataset. Under this model, predictions of some particular quantity of interest, T(y), could be precise, but completely implausible due to the limitations in the model. Consider now an ABC analysis based on this model, based on matching summary statistics that include T(y). ABC methods would identify those parameter values θ that are most likely to have produced these statistics under the model. This means that predictions of T(y) under the ABC approximation now have some chance of being accurate (although they may be less precise), as the model may be able to predict the summary statistics, including T(y), even if it can't accurately predict the full dataset.

This allows ABC to be interpreted as a mechanism for fitting models based on summary statistics that may in fact be more useful than the exact inference with the full dataset. An explicit example of this in the robust model selection  context was given by li+nfs15.

Related arguments allow ABC to be thought of as a natural method to fit models when the full dataset (y_obs) is only partially observed (s_obs) and has missing data (see e.g. rodrigues+fst18). ABC methods have also been used to determine weakly informative prior distributions in a regular tractable Bayesian analysis, exploiting the mechanism of predictive data matching to identify a priori non-viable regions of the parameter space nott+dme17.








§ FURTHER READING






ABC methods have been extensively and rapidly developed since their first modern appearance in tavare+bgd97 and pritch99. Naturally a number of review articles have been written for various discipline audiences to review the techniques available at the time. While with time such reviews can rapidly become dated, they often provide useful perspectives on ABC methods as viewed at the time.  See, for example, the reviews by

beaumont10,
bertorelle+bm10,
blum+nps13,
csillery+bgf10,
sisson+f11,
marin+prr12,
turner+v12,
robert16,
erhardt+s16,
lintusaari+gdkc16
and
drovandi17.
Each of the chapters in this Handbook also makes for excellent reading and review material on focused aspects of ABC tavare17,blum17,fan+s18,prangle17,marin+per18,drovandi18,nott+ofs17,andrieu+lv17,fearnhead18,ratmann+chc18,drovandi+gkr18,wegmann18,barthelme+cc17.



Because ABC methods are now recognised as a standard Bayesian tool, their scientific reach has effectively become as extensive as standard Bayesian methods. While it is accordingly futile to exhaustively describe all areas in which ABC has applied, the below selection is provided to provide a flavour of the impact ABC methods have had. Beyond the applications in this Handbook, ABC methods have been successfully applied to applications in

α-stable models peters+sf12,
archaeology wilkinson+t09,
cell biology johnston+smbr14,vo+dpp15,vo+dps15,
coalescent models fan+k11,tavare+bgd97,
ecology <cit.>,
evolutionary history of mosquitos bennett+slmkdahihplw16,
filtering jasra+smm12,
extreme value theory erhardt+s12,erhardt+s16,
financial modelling peters+sf12,
host-parasite systems baudet+dscgms15,
HIV contact tracing blum+t10,
human evolution fagundes+rbnsbe07,
hydrological models nott+fms14,
infectious disease dynamics luciani+sjft09,aandahl+rst12,
infinite mixture models for biological signalling pathways koutroumpas+bvc16,
image analysis nott+fms14,
long range dependence in stationary processes andrade+r15,
operational risk peters+s06,
quantile distributions allingham+km09,drovandi+p11,
pathogen transmission tanaka+fls06,
phylogeography beaumont+al10,
protein networks ratmann+ahwr09,ratmann+jhsrw07,
population genetics beaumont+zb02,
psychology <cit.>,
single cell gene expression lenive+ks16,
spatial point processes <cit.>,
species migration hamilton+crhbe05,
state space models vakilzadeh+hba17,
stochastic claims reserving peters+fs12,
susceptible-infected-removed (SIR) models toni+wsis09,
trait evolution slater+hwjra12
and
wireless communications engineering peters+nsfy10.
Within this Handbook novel analyses can be found in  peters18, rodrigues+fst18, stumpf18, estoup18, holden18, wood18 and fan18.





§ CONCLUSIONS


ABC methods are based on an inherently simple mechanism – simulating data under the model of interest and comparing the output to the observed dataset. While more sophisticated ABC algorithms and techniques have subsequently been developed (and many of these are discussed in more detail in this Handbook), this core mechanic remains a constant.
It is this methodological simplicity that has made ABC methods highly accessible to researchers in across  many disciplines. We anticipate that this will continue in the future.



§ ACKNOWLEDGMENTS
 

SAS is supported by the Australian Research Council under the Discovery Project scheme (DP160102544), and the Australian Centre of Excellence in Mathematical and Statistical Frontiers (CE140100049).


chicago

empty
