ABC Samplers
Y. Fan S. A. Sisson


Byzantine-Resilient Distributed Observers for LTI Systems footnoteinfo
    [
    December 30, 2023
======================================================================







§ INTRODUCTION


Approximate Bayesian computation (ABC) is a phrase that describes a collection of methods and algorithms designed to perform a Bayesian analysis using an approximation to the true posterior distribution, when the likelihood function implied by the data generating process is computationally intractable.

For observed data y_obs∈𝒴, the likelihood function p(y|θ) depends on a vector of model parameters θ∈Θ, from which prior beliefs π(θ) may be updated into posterior beliefs π(θ|y_obs)∝ p(y_obs|θ)π(θ) via Bayes' theorem.

In the standard ABC framework (see e.g. sisson+fb17, this volume) the ABC approximation to π(θ|y_obs) is given by

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

where K_h(u)=K(u/h)/h is a standard kernel density function with scale parameter h>0, · is an appropriate distance metric (e.g. Euclidean or Mahalanobis distance), p(s|θ)
is the (intractable) likelihood function of the low-dimensional vector of summary statistics s=S(y) implied by p(y|θ), and s_obs=S(y_obs). 
Defining K_h(s-s_obs)→δ_s_obs(s) as h→ 0, 
where δ_Z(z) denotes the Dirac measure, defined as δ_Z(z)=1 if z ∈ Z 
and δ_Z(z)=0 otherwise,
then as a result 

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

Accordingly, π_ABC(θ|s_obs) provides an approximation to the partial posterior π(θ|s_obs), which becomes more accurate as h gets small.
If the summary statistics s are sufficient for θ, then π(θ|s_obs) will equal π(θ|y_obs), and so, for small h,
the ABC posterior approximation π_ABC(θ|s_obs) will be a good approximation of the true posterior. If either s is not sufficient, or h is not small, then the ABC posterior approximation π_ABC(θ|s_obs) will be of the form (<ref>). 

In terms of drawing samples from the approximate posterior π_ABC(θ|s_obs), the choice of summary statistics s=S(y) is typically  considered known, and interest is then in sampling from π_ABC(θ|s_obs), for a specific and low value of h, as efficiently as possible. The more efficient the simulation procedure, the further h can be lowered within the sampling framework, resulting in samples from a more accurate approximation of π(θ|s_obs).

In this chapter we survey the various forms of ABC algorithms that have been developed to sample from π_ABC(θ|s_obs). These have broadly followed the familiar Monte Carlo classes of algorithms, including rejection and importance sampling, Markov chain Monte Carlo (MCMC), and sequential Monte Carlo (SMC) based algorithms. While each of these classes have  their ABC-specific implementations and  characteristics, in general they target the joint distribution of parameter vector θ and summary statistic s given by

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

By noting that (<ref>) is obtained from (<ref>) by integrating over s (i.e. π_ABC(θ|s_obs)=∫π_ABC(θ,s|s_obs)ds), samples from π_ABC(θ|s_obs) can  be obtained by first drawing samples from (<ref>) and then discarding the marginal s values.

An alternative, but  related Monte Carlo approach is based on sampling from (<ref>) directly, by 
obtaining an unbiased and non-negative estimate of the ABC posterior distribution function, such as

    π̂_ABC(θ|s_obs) ∝π(θ)/T∑_t=1^T K_h(s(t)-s_obs),

where s(1),…,s(T)∼ p(s|θ) are samples from the intractable model given θ, and then using this estimate in place of π_ABC(θ|s_obs) within a standard Monte Carlo algorithm (e.g. delmoral+dj12). This approach falls within the family of pseudo-marginal Monte Carlo methods <cit.>, a more general class of likelihood-free samplers that has gained popularity outside of the ABC setting. For a detailed and in-depth discussion of the connections between ABC-MCMC algorithms and pseudo-marginal MCMC methods see andrieu+sv17 (this volume).






§ REJECTION AND IMPORTANCE SAMPLING










 §.§ Rejection sampling




The earliest ABC samplers (e.g. tavare+bgd97,pritchard+spf99) were basic rejection sampling algorithms.

Under the standard rejection sampling framework (e.g.  rip87,liu01) interest is in 
 
obtaining samples from some target distribution f(θ)=Z^-1f̃(θ), which is known up to a normalising constant Z=∫_-∞^∞f̃(θ)dθ. The standard rejection sampling algorithm obtains draws θ'∼ g(θ) from a sampling density g(θ) from which it is trivial to sample,   such that 

    f̃(θ)≤  M g(θ),

for all θ and some positive constant M>0.

The draws θ' are then accepted as  independent samples from the target density f(θ) with probability f̃(θ')/Mg(θ'). 

To see that the above procedure is correct, for simplicity consider the case in which θ is univariate, with the extension to multivariate θ being straightforward.  Define A as the event that a sample θ' from g(θ) is accepted. Then, the overall acceptance rate of the algorithm is

    (A) = ∫f̃(θ)/Mg(θ)g(θ)dθ= 1/M∫f̃(θ) dθ = Z/M,

and hence the distribution of accepted draws is

    (θ'≤θ | A) = (θ'≤θ, A)/(A)=∫_-∞^θf̃(θ')/Mg(θ')g(θ')dθ'/P(A)= F(θ),

as required, where F(θ)=∫_-∞^θ f(z)dz is the distribution function associated with f(θ). 

The efficiency of the algorithm is associated with the value of M, with smaller  values of M (subject to f̃(θ)≤  M g(θ), ∀θ) corresponding to more efficient samplers. That is, for fixed g(θ) the optimum choice is M=max_θf̃(θ)/g(θ). Good choice of the sampling distribution g(θ), e.g. to approximate f(θ), can result in smaller values of M.


The ABC version of the rejection sampler was discussed in sisson+fb17 (this volume), which we reproduce here as Algorithm 1.





Originally developed by pritchard+spf99 following earlier ideas by tavare+bgd97, the ABC rejection sampling algorithm is typically described heuristically as follows: for the candidate parameter vector θ'∼ g(θ), a dataset y' is generated from the (intractable) generative model and summary statistics s'=S(y') computed.
If the simulated and observed datasets are similar (in some manner), so that s'≈ s_obs, then θ' could credibly have generated the observed data under the given model, and so θ' is retained and forms part of the sample from the ABC posterior distribution π(θ |s_obs). 
Conversely, if s' and s_obs are dissimilar, then θ' is unlikely to have generated the observed data for this model, and so θ' is discarded. 

The parameter vectors accepted under this approach offer support for s_obs under the model, and so may be considered to be drawn approximately from the posterior distribution π(θ|s_obs). In this manner, the evaluation of the likelihood p(y_obs | θ'), essential to most Bayesian posterior simulation methods,  is replaced by an evaluation of the proximity of summaries of a simulated dataset s' to the observed summaries s_obs.



More precisely, this algorithm targets π_ABC(θ,s|s_obs) given by (<ref>), the joint distribution of parameter vector and summary statistic given s_obs. Accordingly the sampling distribution is also defined on this space as g(θ,s) = p(s|θ)g(θ), and the acceptance probability of the vector (θ,s) is then given by

    π_ABC(θ,s|s_obs)/M g(θ,s)∝K_h(s-s_obs)p(s|θ)π(θ)/Mp(s|θ)g(θ)
    	=
    	K_h(s-s_obs)π(θ)/Mg(θ).

The normalising constant M is similarly given by

    M≥max_s,  θK_h( s- s_obs)p(s|θ)π(θ)/p(s|θ)g(θ)=K_h(0)max_θπ(θ)/g(θ)

with max_sK_h(s-s_obs)=K_h(0) resulting from  the zero-mean, symmetry and (typically) unimodal characteristics of standard kernel density functions. Accordingly, the construction of the target and sampling distributions on the joint space (θ,s) results in the form of the acceptance probability and normalisation constant M being free of intractable likelihood terms.

An example implementation of this algorithm is given in sisson+fb17 (this volume).







 §.§ Importance sampling




One down side of rejection sampling is the need to determine a near optimal value for the normalising constant M in order to produce an efficient algorithm. Importance sampling is a procedure that, rather than calculating acceptance probabilities, avoids this by alternatively assigning the draw θ'∼ g(θ) an (importance) weight w(θ')= f(θ')/g(θ'). The weighted vector θ' is then a draw from f(θ),  and desired expectations under the target distribution f are computed as weighted expectations under the importance sampling density g.

To see this, suppose that we are interested in estimating the expectation

    E_f[h(θ)] = ∫ h(θ)f(θ)dθ.
 
By defining w(θ)=f(θ)/g(θ) we have

    E_g[w(θ)h(θ)] = ∫ w(θ)h(θ)g(θ)dθ 
    	= ∫ h(θ)f(θ)dθ = E_f[h(θ)].

In this manner we can then estimate the expectation as

    E_f[h(θ)] ≈1/N∑_i=1^Nw^(i)h(θ^(i)),

where w^(i)=w(θ^(i)), and where θ^(i)∼ g(θ) are draws from g. In the more typical case where the target distribution is unnormalised, so that f(θ)=Z^-1f̃(θ),  we can work with f̃(θ) only by defining w̃(θ)=f̃(θ)/g(θ) and then noting that

    E_g[w̃(θ)] = ∫w̃(θ)g(θ)dθ = ∫f̃(θ)dθ = Z ≈1/N∑_i=1^Nw̃^(i),

for θ^(i)∼ g(θ), where w̃^(i)=w̃(θ^(i)). As a result, the expectation E_f[h(θ)] may be approximated as

    E_f[h(θ)]     =    ∫ h(θ)f(θ) = 1/Z∫ h(θ)f̃(θ)
        =    1/Z∫w̃(θ)h(θ)g(θ)dθ = 
    	E_g[w̃(θ)h(θ)]/E_g[w̃(θ)]
       ≈   1/N∑_i=1^Nw̃^(i)h(θ^(i))/1/N∑_i=1^Nw̃^(i)
    	= ∑_i=1^N W^(i)h(θ^(i)),

for θ^(i)∼ g(θ), where W^(i) = w̃^(i)/∑_j=1^N w̃^(j) denotes normalised weights. This approximation is not unbiased due to the biased estimator of 1/Z, although the bias becomes small as N becomes large.





From an ABC perspective, importance sampling works much the same as rejection sampling. The
target distribution is π_ABC(θ,s|s_obs), and the importance distribution on joint parameter value and summary statistics space is g(θ,s)=p(s|θ)g(θ). As a result, the (unnormalised) importance weights are computed as

    π_ABC(θ,s|s_obs)/g(θ,s)∝K_h(s-s_obs)p(s|θ)π(θ)/p(s|θ)g(θ) =  K_h(s-s_obs)π(θ)/g(θ) := w̃(θ),

which is again free of intractable likelihood terms. The full ABC importance sampling algorithm is given in Algorithm 2.




As with rejection sampling, the 
choice of the (marginal) importance distribution g(θ) is crucial to the efficiency of the algorithm. 

In standard importance sampling, if g(θ)∝f̃(θ) then w̃(θ)∝ 1. In this case, there is no variation in the importance weights, and each sample θ^(i) contributes equally when computing posterior expectations via (<ref>). However, if g(θ) is different to f̃(θ) then the variability in w̃(θ) from θ means that some samples θ^(i) will contribute more than others in this computation. In extreme cases, Monte Carlo estimates of expectations can be highly variable when they are dominated by a small number of θ^(i) with relatively large weights. This is known as sample degeneracy. Accordingly, for importance sampling algorithms, the focus is on reducing the variability of w̃(θ) over θ.


A common measure of the degree of
sample degeneracy is the effective sample size (ESS) liucw98,liu01, 
estimated as

    ESS = (∑_i=1^N[W^(i)]^2)^-1, 
    	   1≤ ESS≤ N,

which is computed using the normalised weights  W^(i) = w̃^(i)/∑_j=1^N w̃^(j).

The ESS is an estimate of the effective number of equally weighted θ^(i) in a given weighted sample, which can be loosely interpreted as the information content.
When g(θ)∝f̃(θ) so that we have samples directly from f(θ),  then W^(i)=1/N and ESS=N.  However, when there is severe particle degeneracy in the extreme case where  W^(1)=1 and W^(i)=0 for i=2,…,N, then ESS=1.



Specifically in the ABC framework where w̃(θ^(i))=K_h(s^(i)-s_obs)π(θ^(i))/g(θ^(i)), if g is diffuse compared to the (marginal) target distribution π_ABC(θ|s_obs), samples θ^(i)∼ g(θ) from regions of  low posterior density will tend to generate summary statistics s^(i) that are very far from the observed statistics s_obs, and this will produce low weights w^(i) compared to samples θ^(i) in regions of high posterior density. This is the same as for the standard importance sampling case. However, in ABC importance sampling, an additional factor is that the importance weight w̃(θ) is a function of the kernel function K_h(s-s_obs), which contains the stochastic term s. This has some implications, which are also relevant for sequential Monte Carlo-based ABC samplers, discussed in Section <ref>.


When K_h has non-compact support, such as when K_h(u)=ϕ(u; 0, h^2), where ϕ(x; μ,σ^2) denotes the Gaussian density function with mean μ and variance σ^2, the importance weight w̃^(i) is guaranteed to be non-zero for each i. However the resulting importance weight can be highly variable, depending on whether s^(i) is close to or far from s_obs. This typically produces samples (θ^(i), w̃^(i)) with low effective sample sizes.

If K_h has a compact support (and this is typical in most ABC implementations), then w̃^(i)=0 is 
likely for small h, even when θ^(i) is in a high posterior density region. This means that Algorithm 2 will return many (θ^(i),w̃^(i)) for which the weight is exactly zero, resulting in low effective sample sizes, and maybe even complete algorithm failure if w̃^(i)=0 for all i=1,…,N. As a result, a common variation of Algorithm 2 is to repeat steps 1–4 for each i, until a non-zero weight has been generated. 

This effectively introduces a rejection sampling step within the importance sampling algorithm. This idea (c.f. fearnhead+p12) can be used to improve the ESS for ABC importance sampling algorithms, regardless of the choice of the kernel, by modifying step 4 in Algorithm 2 to be



When K_h has compact support, this ensures that steps 1–4 of Algorithm 2 are repeated until K_h(s^(i)-s_obs) is non-zero. When K_h has non-compact support, this offers some control over the variability of the weights, as only samples for which s^(i) is reasonably close to s_obs are likely to be accepted.







Under Algorithm 3, in the particular case of when K_h is the uniform kernel on [-h,h], if in addition g(θ)∝π(θ) so that the importance distribution is proportional to the prior, then w̃^(i)∝ 1 for any i. This results in W^(i)=1/N and ESS=N and the ABC importance sampling algorithm
effectively reduces to the ABC rejection sampling algorithm, but without the need to compute the normalising constant M. This setup is very common in practice as it removes the need to compute importance weights, and to worry about algorithm performance with respect to effective sample size, which is always maximised. However, in this case, algorithm performance is dominated by the number of times steps 1–4 are repeated before a sample θ^(i) is accepted. In general the efficiency of Algorithm 3 is a combination of the resulting effective sample size and the number of repetitions of the sampling steps 1–4.









 §.§ Importance/rejection sampler variants


There are many variants on ABC importance and rejection samplers. A few of these are detailed below, chosen either because of their popularity, or because of their links with particular ABC samplers discussed in later Sections.





  §.§.§ Rejection control importance sampling






 liucw98 developed a general importance-rejection algorithm technique known as rejection control, with the aim of reducing the number of θ^(i) samples that are produced with very small weights in an importance sampler. This method was exploited within an ABC sequential Monte Carlo framework by sisson+ft07 and peters+fs12 (see Section <ref>), however it may also be implemented directly within an ABC importance sampler as outlined below.
  
Suppose that a weighted sample (θ^(i),w̃^(i)) is drawn from f(θ) using an importance sampling algorithm. In order to control the size of the importance weight, w̃^(i) is compared to some pre-specified threshold value c>0. If w̃^(i)>c then the weight is considered sufficiently large, and the sample θ^(i) is accepted. However if w̃^(i)<c then θ^(i) is probabilistically rejected, with a higher rejection rate for lower w̃^(i). In this manner, the variability of the accepted importance weights can be reduced.

In particular, each sample θ^(i) 
is accepted with probability

    r^(i)=min{1, w̃^(i)/c},

which results in the automatic acceptance of samples for which w̃^(i)>c and an acceptance probability of w̃^(i)/c otherwise. This means that larger c results in less variable weights, although at the price of more rejections.

The accepted samples are then draws from the modified importance sampling distribution

    g^*(θ) = M^-1min{1, w̃(θ)/c} g(θ)

where  w̃(θ)=f̃(θ)/g(θ), and with normalising constant M=∫min{1, w̃(θ)/c} g(θ) dθ. As a result, setting

    w̃^*(θ)=f̃(θ)/g^*(θ)    w̃^*(i)=f̃(θ^(i))/g^*(θ^(i)) = M w̃^(i)/r^(i)

means that the samples (θ^(i),w̃^*(i)) will be weighted samples from f(θ) but with the property that 

    Var_g^*[f̃(θ)/g^*(θ)] ≤ Var_g[f̃(θ)/g(θ)].

That is, the rejection control algorithm can reduce the variance of the importance weights <cit.>.

While it may be difficult to evaluate M = E_g[min{1, w̃(θ)/c}] analytically, it may be estimated from the samples (θ^(i),w̃^(i)) via

    M̂≈1/N∑_i=1^N min{1, w̃^(i)/c}.

If an estimate of M is not required, its computation can be avoided for importance sampling purposes by calculating the normalised weights W^*(i) = w̃^*(i)/∑_j=1^Nw̃^*(j), as the M term then cancels in numerator and denominator.



As with ABC rejection sampling (Algorithm 2), the ABC implementation of rejection importance control  targets π_ABC(θ,s|s_obs) resulting in a weight calculation of w̃^(i)=K_h(s^(i)-s_obs)π(θ^(i))/g(θ^(i)). The full algorithm is given in Algorithm 4.


Note that while Algorithm 4 requires pre-specification of the rejection threshold c, a suitable value may be practically difficult to determine in advance. As such, Algorithm 4 may be alternatively executed by first implementing steps 1–3 only for i=1,…,N, and then specifying c as some quantile of the resulting empirical distribution of w̃^(1),…,w̃^(N). Following this, Algorithm 4 may then continue implementation from step 4 onwards for each i=1,…,N (e.g. peters+fs12).







As with ABC importance/rejection sampling (Algorithm 3), when K_h has compact support, rejection control will replace those samples for which the simulated and observed summary statistics are too far apart, resulting in w̃^(i)=0. More generally, however, rejection control provides much greater control over the variability of the weights regardless of K_h, producing more uniform weights for larger c. The price for this control is the  greater number of rejections induced as c increases peters+fs12.






  §.§.§ k-nearest neighbour ABC importance sampling





While most published descriptions of importance and rejection sampling ABC algorithms follow the format given in Algorithms 1–4, in practice it is not uncommon to deviate from these and implement a slight variation. The reason for this is that Algorithms 1–4 require pre-specification of the kernel scale parameter h>0, without which importance weights cannot be calculated and accept/reject decisions cannot be made. In reality, as the scale of the distances s^(i)-s_obs is unlikely to be known in advance, it is difficult to pre-determine a suitable value for h.

Algorithm 5 presents a variation on the ABC importance sampler of Algorithm 3 that avoids pre-specification of h. Here a large number N' of (θ^(i),s^(i)) pairs are generated from the importance sampling distribution p(s|θ)g(θ). These are the only samples that will be used in the algorithm, so the computational overheads are fixed at N' draws from the model, unlike Algorithm 3 in which the number of draws is random and unknown in advance. The N samples for which s^(i) is closest to s_obs (as measured by ·) are then identified, and h determined to be the smallest possible value so that only these N samples have non-zero weights (assuming a kernel K_h with compact support). Once h is fixed, the importance weights can be calculated as before, and the N samples (θ^(i),w̃^(i)) with non-zero w̃^(i) are returned as weighted samples from π_ABC(θ|s_obs). 



This approach is explicitly used in e.g. beaumont+zb02 and blum+nps13, and implicitly in many other ABC implementations. The differences between Algorithms 3 and 5 may seem small – if the value of h
determined in Algorithm 5 was used in Algorithm 3, then (assuming the same pseudo-random numbers used in the appropriate places) the resulting draws from π_ABC(θ|s_obs) would be identical. However, Algorithm 5 is based on a k-nearest neighbour algorithm for density estimation of the ABC likelihood function, and so possesses very different theoretical properties compared to Algorithm 3. This k-nearest neighbour approach is discussed and analysed in detail in the rejection sampling context by biau+cg15.





  §.§.§ ABC rejection sampling with stopping rule


A version of ABC rejection sampling (Algorithm 1) which similarly does not require pre-specification of the kernel scale parameter h>0 is presented in Example 6 (Section 7.2) of sisson+fb17 (this volume). We do not reproduce this algorithm here for brevity.

The algorithm identifies the smallest value of h needed to accept exactly N samples before some stopping rule is achieved. This stopping rule could be based on an overall computational budget (such as using exactly N' total draws from p(s|θ)), or on some perceived level of accuracy of the resulting ABC posterior approximation.

If the stopping rule is based on an overall computational budget of exactly N' draws from p(s|θ) (and again the same pseudo-random numbers), this algorithm will produce exactly the same final samples from π_ABC(θ|s_obs) as Algorithm 1, were the ABC rejection sampler to adopt the identified choice of h. Of course, the advantage here is that the value of h is automatically determined.





  §.§.§ Rejection-based ABC algorithms for expensive simulators



It is not uncommon for the data generation step y^(i)∼ p(y|θ^(i)) in ABC algorithms to be expensive, and thereby dominate the computational overheads of the algorithms. While there are a few principled ways to mitigate this (see discussion of prangle+ek17 and everitt+r17 in Section <ref>), within rejection-based ABC algorithms it is sometimes possible to reject a proposed sampler θ^(i) before generating the data y^(i)∼ p(y|θ^(i)). 

To see this, note that e.g. step 4 of Algorithm 3


  4. With probability K_h(s^(i)-s_obs)/K_h(0) set w̃^(i)=π(θ^(i))/g(θ^(i)), else go to 1.

can be alternatively implemented as 


  4. With probability π(θ^(i))/g(θ^(i)) set w̃^(i)=K_h(s^(i)-s_obs)/K_h(0), else go to 1.

This means that steps 2 and 3 of Algorithm 3 (generate y^(i)∼ p(y|θ^(i)) and compute s^(i)=S(y^(i))) need not be performed until the event in step 4 with probability π(θ^(i))/g(θ^(i)) has occurred.
This allows for a possible early rejection of θ^(i) before any data generation needs to take place. (Note that if g(θ)=π(θ) there is no benefit to be gained.) This modification trades some computational savings for weights w̃^(i) constructed from different terms, and thereby  having different variance properties.

This idea, which is a standard technique in standard sequential Monte Carlo samplers (e.g. delmoral06), can be implemented in any rejection-based ABC algorithm, including ABC-MCMC and ABC-SMC samplers (Sections <ref> and <ref>).






  §.§.§ Marginal ABC samplers




Until now we have presented ABC algorithms as producing samples (θ^(i),s^(i)) exactly from the joint distribution π_ABC(θ,s|s_obs)∝ K_h(s-s_obs)p(s|θ)π(θ).
As a result, samples θ^(i) from the ABC approximation to the posterior π(θ|y_obs) given by

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

may be obtained by marginalising over the realised s^(i).

An alternative approach to construct an ABC algorithm could be to directly target the (marginal) posterior π_ABC(θ|s_obs) rather than the joint posterior π_ABC(θ,s|s_obs). This approach becomes apparent when noting that π_ABC(θ|s_obs) can be estimated pointwise (up to proportionality), for fixed θ, as

    ∫ K_h(s-s_obs)p(s|θ)π(θ) ds≈π(θ)/T∑_t=1^TK_h(s(t)-s_obs) := π̂_ABC(θ|s_obs),

where s(1),…,s(T)∼ p(s|θ) are T independent draws of summary statistics from the partial likelihood p(s|θ) for a given θ. This Monte Carlo estimate of π_ABC(θ|s_obs) is unbiased up to proportionality (in that E_s|θ[π̂_ABC(θ|s_obs)]∝π_ABC(θ|s_obs)), and so π̂_ABC(θ|s_obs) may be used in place of π_ABC(θ|s_obs) in a standard rejection or importance sampler which targets π_ABC(θ|s_obs). Using this substitution will produce a random, estimated acceptance probability or importance weight. However, because it is also unbiased (up to proportionality), the resulting target distribution will remain the same as if the exact weight had been used i.e. π_ABC(θ|s_obs), although the sampler weights/acceptances will become more variable. This is the so-called marginal ABC sampler (e.g. marjoram+mpt03, 	reeves+p05, sisson+ft07, ratmann+ahwr09, toni+wsis09, peters+fs12, among others).

As with standard Monte Carlo estimates, the number of Monte Carlo draws T affects the variability of the ABC posterior estimator. bornn+psw17 explore the question of how many draws, T, produces the most efficient overall sampler in the context of ABC rejection and Markov chain Monte Carlo algorithms. If T is large, the estimate of π_ABC(θ|s_obs) is accurate and so the acceptance probability is accurate but at the cost of many Monte Carlo draws, however if T is small, the acceptance probability is highly variable but is much cheaper to evaluate. When using a uniform kernel K_h, bornn+psw17  conclude that in fact, T=1 is the most efficient, as (loosely) the combination of T draws used to accept one θ^(i) could be better used to accept up to T different θ^(i)'s, each using one Monte Carlo draw per ABC posterior estimate.

The idea of the marginal ABC sampler is closely related to the construction of the more recently developed pseudo-marginal sampler 
<cit.>, a more general class of likelihood-free sampler that has gained popularity outside of the ABC setting. 

Here, rather than treating π̂_ABC(θ|s_obs) as an unbiased estimate of π_ABC(θ|s_obs) in an algorithm that targets π_ABC(θ|s_obs), an alternative joint posterior distribution can be constructed

    π_ABC(θ,s(1),…,s(T)|s_obs)∝[1/T∑_t=1^TK_h(s(t)-s_obs)]
    	[∏_t=1^Tp(s(t)|θ)]
    	π(θ),

which is defined over the joint posterior of θ and all T summary statistic replicates (e.g. delmoral+dj12,sisson+f11), where T=1 gives the usual ABC joint posterior π_ABC(θ,s|s_obs). A useful property of this form of joint posterior is that the θ-marginal distribution is the same for any value of T, and in particular

    ∫…∫π_ABC(θ,s(1),…,s(T)|s_obs) ds(1)… ds(T)
    	=
    	π_ABC(θ|s_obs).

This means that any sampler targeting π_ABC(θ,s(1),…,s(T)|s_obs) can produce samples from π_ABC(θ|s_obs). 

Consider now an importance sampler targeting π_ABC(θ,s(1),…,s(T)|s_obs) with the importance sampling density

    g(θ,s(1),…,s(T)) = g(θ)∏_t=1^Tp(s(t)|θ).

The resulting importance weight is

    [                                 π_ABC(θ,s(1),…,s(T)|s_obs)/g(θ,s(1),…,s(T)) ∝[1/T∑_t=1^TK_h(s(t)-s_obs)]
    	[∏_t=1^Tp(s(t)|θ)]
    	π(θ)/g(θ)∏_t=1^Tp(s(t)|θ);                                                                                                                                =
    	π̂_ABC(θ|s_obs)/g(θ). ]

This means that any marginal ABC sampler targeting π_ABC(θ|s_obs) through the unbiased estimate of the ABC posterior given by π̂_ABC(θ|s_obs) is directly equivalent to an exact algorithm targeting π_ABC(θ,s(1),…,s(T)|s_obs). That is, all marginal ABC samplers are justified by their equivalent joint space ABC algorithm.

This idea also extends to using unbiased approximations of posterior distributions within MCMC samplers (see next Section), where the technique has expanded beyond ABC algorithms to more general target distributions. Here it is more generally known as pseudo-marginal Monte Carlo methods.

See andrieu+sv17 (this volume)
for a more detailed discussion of the connections between ABC marginal samplers and pseudo-marginal MCMC methods.









§ MARKOV CHAIN MONTE CARLO METHODS








Markov chain Monte Carlo (MCMC) methods are a highly accessible class of algorithms for obtaining samples from complex distributions (e.g. brooksgjm11). By constructing a Markov chain with the target distribution of interest as its limiting distribution, following chain convergence,
a realised random sample path from this chain will behave like a (serially correlated) sample from the target distribution,. Their strong performance and simplicity of implementation has made MCMC algorithms the dominant Monte Carlo method for the past two decades brooksgjm11.

As such, it is only natural that MCMC-based ABC algorithms have been developed.







 §.§ ABC MCMC samplers


The Metropolis-Hastings algorithm is the most popular class of MCMC algorithm. Given the current chain state θ^(i), the next value in the sequence is obtain by sampling a candidate value θ' from a proposal distribution θ' ∼ g(θ^(i), θ)=g(θ|θ^(i)), which is then accepted
with probability a(θ,θ')=min{1, f(θ')g(θ', θ^(i)) /f(θ)g(θ^(i), θ')} so that θ^(i+1)=θ', or otherwise rejected so that θ^(i+1)=θ^(i). Under this mechanism the target distribution is f(θ), and there is great flexibility in the choice of the proposal distribution g. 

An implementation of this sampler in the ABC setting is given in Algorithm 6.

ABC MCMC algorithms were originally developed by  marjoram+mpt03. See e.g.   bortot+cs07, wegmann+le09, ratmann+ahwr09, sisson+f11 and andrieu+sv17 (this volume)
for more discussion on ABC MCMC samplers.




As with ABC importance and rejection samplers, the target distribution of ABC MCMC algorithms is the joint ABC posterior π_ABC(θ,s|s_obs). On this space the proposal distribution becomes

    g[(θ,s),(θ',s')] =g(θ,θ') p(s'|θ'),

and as a result the acceptance probability of the proposed move from (θ^(i), s^(i)) to (θ', s')∼ g[(θ^(i),s^(i)),(θ',s')] becomes
a[(θ^(i),s^(i)),(θ',s')] = min{1,α[(θ^(i),s^(i)),(θ',s')]}, where

    α[(θ^(i),s^(i)),(θ',s')]
    	   =   π_ABC(θ',s'|s_obs)g[(θ',s'),(θ^(i),s^(i))]/π_ABC(θ^(i),s^(i)|s_obs)g[(θ^(i),s^(i)),(θ',s')]
       =   K_h(s'-s_obs)p(s'|θ')π(θ')/K_h(s^(i)-s_obs)p(s^(i)|θ^(i))π(θ^(i))g(θ',θ^(i)) p(s^(i)|θ^(i))/g(θ^(i),θ') p(s'|θ')
       =   K_h(s'-s_obs)π(θ')/K_h(s^(i)-s_obs)π(θ^(i))g(θ',θ^(i))/g(θ^(i),θ'),


which is free of intractable likelihood terms, p(s|θ), and so may be directly evaluated.







Algorithm 6 satisfies the detailed balance (time reversibility) condition with respect to π_ABC(θ,s|s_obs), which ensures that π_ABC(θ,s|s_obs) is the stationary distribution of the Markov chain. Detailed balance states that

    π_ABC(θ, s|s_obs)P[(θ, s),(θ', s')]  = π_ABC(θ',s'|s_obs)P[(θ', s'),(θ, s)],

where
the Metropolis-Hastings transition kernel P is given by

    P[(θ, s),(θ', s')]=g[(θ, s),(θ', s')]a[(θ, s),(θ', s')].

Assuming that (without loss of generality) a[(θ', s'),(θ, s)]=min{1,α[(θ',s'),(θ,s)]}=1 (and so a[(θ, s),(θ', s')]=α[(θ,s),(θ',s')]), the
detailed balance condition  is satisfied since

    π_ABC(θ, s|s_obs)P[(θ, s),(θ', s')]
        =   π_ABC(θ, s|s_obs)g[(θ, s),(θ', s')]α[(θ, s),(θ', s')]
       =   K_h(s-s_obs||)p(s|θ)π(θ)/zg(θ,θ')p(s'|θ')K_h(s'-s_obs)π(θ')g(θ',θ)/K_h(s-s_obs)π(θ)g(θ,θ')
       =   K_h(s'-s_obs)p(s'|θ')π(θ')/zg(θ',θ)p(s|θ)
       =   π_ABC(θ',s'|s_obs)P[(θ', s'),(θ, s)],

where z=∫∫ K_h(s-s_obs)p(s|θ)π(θ)dsdθ is the normalisation constant of π_ABC(θ,s|s_obs) (e.g. sisson+f11).


marjoram+mpt03 found that the ABC MCMC algorithm offered an improved acceptance rate over rejection sampling-based ABC algorithms with the same scale parameter h, although at the price of serial correlation in the Markov chain sample path θ^(1),…,θ^(N).

Thus, for kernels K_h with compact support, the same mechanism that causes many rejections or zero weights in ABC rejection and importance samplers, now results in many rejected proposals in the ABC MCMC algorithm. The difference here is that the chain simply remains at the current state θ^(i) for long periods of time, giving additional posterior weight to θ^(i).

Techniques for improving the performance of standard MCMC algorithms may also be applied to ABC MCMC samplers. However there is one feature of ABC MCMC that is different to that of the standard algorithm, that is particularly acute when using kernel functions K_h with compact support.


Consider a proposed move from θ^(i) to θ'. In standard MCMC, the acceptance probability is based on the relative density of the posterior evaluated at θ' compared to that evaluated at θ^(i). In ABC MCMC the density of the posterior at θ' is determined through the ability of the model to generate a summary statistic s'∼ p(s|θ') that is close to s_obs as measured through K_h(s'-s_obs). That is, to move to θ', a summary statistic s' must be generated that is close enough to s_obs. This is the standard ABC mechanism. However, the result of this for the ABC MCMC algorithm is that it means that the acceptance rate of the sampler is directly related to the value of the (intractable) likelihood function evaluated at θ'. As a result, the sampler may mix rapidly in regions of high posterior density, but will have much worse mixing in regions of relatively low posterior density sisson+ft07. For this reason, ABC-MCMC samplers can often get stuck in regions of low posterior density for long periods time, effectively producing convergence issues for the algorithm. 

This effect is more pronounced when the kernel K_h has compact support, such as the uniform kernel on [-h,h] which is endemic in ABC implementations, although it is still present for kernels defined on the real line, such as the Gaussian density kernel. 

In a study of  `sojourn time' within ABC MCMC samplers (that is, the number of consecutive iterations in the sampler in which a univariate parameter θ remained above some high threshold), sisson+f11 found empirically that samplers with uniform kernels had a substantially higher expected sojourn time than samplers with Gaussian kernels, indicating that the latter had superior chain mixing in distributional tails. Despite this, ABC MCMC samplers are routinely implemented with uniform kernels K_h. 



Chain mixing can be improved by alternatively targeting the joint posterior distribution π_ABC(θ,s(1),…,s(T)|s_obs) given by (<ref>). Under this (pseudo) marginal sampler framework (Section <ref>) as T→∞, the mixing properties of the ABC MCMC  approach that of the equivalent standard MCMC sampler directly targeting π_ABC(θ|s_obs) (if it would be possible to numerically evaluate the density function). sisson+f11 empirically demonstrated this improvement, as measured in sojourn times, as T increases. Of course, this improvement of chain mixing comes at the price of overall sampler performance as the computational overheads of generating s(1),…,s(T) for large T would be extremely high. The results of bornn+psw17, that T=1 is the optimum efficiency choice for uniform kernels K_h, also hold for ABC MCMC samplers.







 §.§ Augmented space ABC-MCMC samplers





The standard ABC MCMC sampler (Algorithm 6) requires pre-specification of the kernel scale parameter h. As with ABC rejection and importance samplers, there are a number of ways in which lack of knowledge of a suitable value for the kernel scale parameter can be incorporated into the basic algorithm. Most of these methods also attempt to improve chain mixing over the standard algorithm which uses a fixed, low value of h.

At the very simplest level,  this could involve adaptively adjusting h as a function of s-s_obs at the current and proposed  states of the chain, and either allow h to slowly reduce to some target value to improve convergence at the start of the sampler (e.g. ratmann+jhsrw07, sisson+f11, p.325), or adaptively choose h to achieve some pre-determined overall sampler acceptance probability.


Augmenting the dimension of the target distribution is a common strategy to improve the performance of Monte Carlo algorithms.

In order to help the ABC MCMC sampler escape from regions of low posterior density,  bortot+cs07 proposed augmenting the joint ABC posterior π_ABC(θ,s|s_obs) to additionally include the kernel bandwidth h,  treating this as an unknown additional parameter. The resulting joint posterior distribution is given by

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

and the resulting ABC approximation to the partial posterior π(θ|s_obs) is then given by

    π̇_ABC(θ|s_obs)=∫∫π_ABC(θ,s,h|s_obs)dsdh

where h>0. 

Here h is treated as a tempering parameter in the manner of simulated tempering <cit.>, with larger and smaller values respectively corresponding to “hot” and “cold” tempered posterior distributions. 

Larger values of h increase the scale of the kernel density function K_h, under which the sampler is more likely to accept proposed moves and thereby alleviating the sampler's mixing problems, although at the price of a less accurate posterior approximation. Lower values of h produce a more accurate posterior approximation, but will induce slower chain mixing.

The density π(h) is a pseudo-prior, which serves to influence the mixing of the sampler through the tempered distributions.



 



Note that the augmented space ABC posterior approximation π̇_ABC(θ|s_obs) given by (<ref>) will in general be different to that of π_ABC(θ|s_obs) as the latter contains a fixed value of h, whereas the former integrates over the uncertainty inherent in this parameter. 

Rather than use (<ref>) as the final ABC approximation to π(θ|s_obs), bortot+cs07 chose to remove those samples (θ^(i),s^(i),h^(i)) for which h^(i) was considered too large to come from a good approximation to π(θ|s_obs). In particular, they examined the distribution of θ^(i)|h^(i)≤ h^*, aiming to choose the largest value of h^* such that the distribution of θ^(i)|h^(i)≤ h^* did not change if h^* was reduced further. The resulting ABC posterior approximation is therefore given by

    π̈_ABC(θ|s_obs)=∫_0^h^*∫π_ABC(θ,s,h|s_obs)dsdh.


This approach effectively permits an a posteriori evaluation of an appropriate value h^* such that  the approximation π̈_ABC(θ|s_obs) is as close as possible (subject to Monte Carlo variability) to the true posterior π(θ|s_obs).

A similar idea was explored by baragatti+gp13 in an ABC version of the parallel tempering algorithm of geyer+t95. Here M>1 parallel ABC MCMC chains are implemented with different kernel density scale parameters h_M<h_M-1<…<h_1, with state transitions allowed between chains so that the states of the more rapidy mixing chains (with higher h values) can propagate down to the more slowly mixing chains (with lower h). The final ABC posterior approximation is the output from the chain with h=h_M.

A related augmented space ABC sampler based on the equi-energy MCMC sampler of kou+zw06 could similarly be implemented.




ratmann+ahwr09 take the auxiliary space ABC sampler of bortot+cs07 beyond the solely mechanical question of improving Markov chain mixing, and towards estimation of the distribution of s-s_obs under the model. This is more in line with the ABC approximation π̇_ABC(θ|s_obs) given by (<ref>), and the interpretation of the ABC approximation to π(θ|s_obs) as an exact model in the presence of model error due to wilkinson13. It additionally allows an assessment of model adequacy. Instead of comparing s to s_obs through K_h(s-s_obs) with a single h, ratmann+ahwr09 alternatively make the comparison independently and univariately for each of the q summary statistics in s=(s_1,…,s_q)^⊤  via K_h_r(τ_r-|s_r-s_obs,r|) for r=1,…,q. Here, τ_r is the parameter denoting the true but unknown discrepancy between the r-th summary statistics of s and s_obs, i.e. |s_r-s_obs,r|, and so if τ_r=0 then the model can adequately explain the observed data as described through the r-th summary statistic. The full model has a joint target distribution of

    π_ABC(θ,s(1),…,s(T),τ|s_obs) 
       ∝   min_r[1/Th_r∑_t=1^T K_h_r(τ_r-|s_r(t)-s_obs,r|)]
    	[∏_t=1^T p(s(t)|θ)]
    	π(θ)π(τ),

based on T samples s(1),…,s(T)∼ p(s|θ), where s_r(t) is the r-th element of s(t), and π(τ)=∏_r=1^qπ(τ_r). 

The minimum over the univariate density estimates aims to focus the model on the most conservative estimate of model adequacy, while also reducing computation over τ to its univariate margins.

Here interest is in the posterior distribution of τ in order to determine model adequacy (i.e. if the posterior marginal distribution of π_ABC(τ_r|s_obs) is centered on 0), whereas the margin specific kernel scale parameters h_r are determined via standard kernel density estimation arguments over the observed sample |s_r(t)-s_obs,r| for t=1,…,T.









 §.§ Other ABC MCMC samplers




The field of MCMC research with tractable target distributions is fairly mature, and it is not difficult to imagine that many known techniques can be directly applied to ABC MCMC algorithms to improve their performance. Different forms of algorithms include 

Hamiltonian Monte Carlo ABC samplers meeds+lw15 which use a moderate number of simulations under the intractable model to produce an ABC estimate of the otherwise intractable gradient of the potential energy function,

multi-try Metropolis ABC (aandahlThesis, kobayashi+k15) which uses multiple proposals to choose from at each stage of the sampler to ensure improved mixing and acceptance rates,

in addition to the various augmented space samplers discussed in the previous Section bortot+cs07,ratmann+ahwr09,baragatti+gp13. 

Of course, transdimensional ABC MCMC samplers can also be implemented for multi-model posterior inference. 

General improvements in efficiency can be obtained by using quasi Monte Carlo ABC methods to form efficient proposal distributions cabras+nr15.

In a similar manner, neal12 developed a coupled ABC MCMC sampler which uses the same random numbers to generate the summary statistics for different parameter values, and showed this algorithm to be more efficient than the standard ABC MCMC sampler.




Within the standard ABC MCMC sampler, kousathanas+lhqfw16 proposed using a subset s̃⊆ s of the vector of summary statistics within the acceptance probability when updating a subset of the model parameters conditional on the rest. Here the idea was to reduce the dimension of the comparison s̃-s̃_obs within the kernel K_h to increase the efficiency and mixing of the algorithm. rodrigues17 developed a related algorithm based on the Gibbs sampler.



lee+l14 present an analysis of the variance bounding and geometric ergodicity properties of three reversible kernels used for ABC MCMC, previously suggested by lee+ad12, which are based on the uniform kernel K_h.
Given that current state of the chain is θ^(i) and a proposed new state is drawn from θ'∼ g(θ^(i),θ), the following algorithms were examined (where I(·) denotes the indicator function):


  * Method 1:
Draw s'(1),…,s'(T)∼ p(s|θ').

Accept the move θ^(i+1)=θ' (and s(t)=s'(t) ∀ t) with probability

    min{1,[∑_t=1^T I(s'(t)-s_obs≤ h)]π(θ')g(θ',θ^(i))/[∑_t=1^TI(s(t)-s_obs≤ h)]π(θ^(i))g(θ^(i),θ')}

else reject and set θ^(i+1)=θ^(i).

  * Method 2:
Draw s(1),…,s(T-1)∼ p(s|θ^(i)) and s'(1),…,s'(T)∼ p(s|θ').

Accept the move θ^(i+1)=θ' with probability

    min{1,[∑_t=1^T I(s'(t)-s_obs≤ h)]π(θ')g(θ',θ^(i))/[1+∑_t=1^T-1I(s(t)-s_obs≤ h)]π(θ^(i))g(θ^(i),θ')}

else reject and set θ^(i+1)=θ^(i).

  * Method 3:
Reject the move and set θ^(i+1)=θ^(i) with probability

    1-min{1,π(θ')g(θ',θ^(i))/π(θ^(i))g(θ^(i),θ)}.

For T=1,2,… 
draw s(T)∼ p(s|θ^(i)) and s'(T)∼ p(s|θ')
until ∑_t=1^TI(s(t)-s_obs≤ h)+I(s'(t)-s_obs≤ h)≥ 1.


If I(s'(T)-s_obs≤ h)=1 then set θ^(i+1)=θ' else set θ^(i+1)=θ^(i).


Method 1 is the acceptance probability constructed from the standard Monte Carlo estimate of the ABC posterior π̂_ABC(θ|s_obs) using a fixed number, T, of summary statistic draws, as described in Section  <ref>. Method 2 is the same as Method 1, except that T-1 of the summary statistics of the current chain state s|θ^(i) are regenerated anew in the denominator of the acceptance probability. The idea here is to help the Markov chain escape regions of low posterior probability more easily than under Method 1, at the cost of higher computation. Method 3 produces a random number of summary statistic generations, with computation increasing until either s|θ^(i) or s'|θ' is sufficiently close to s_obs.

Under some technical conditions, lee+l14 conclude that Methods 1 and 2 cannot be variance bounding, and that Method 3 (as with the standard Metropolis-Hastings algorithm if it were analytically tractable) can be both variance bounding and geometrically ergodic. Overall these results, 
in addition to other methods for constructing estimates of intractable likelihoods (e.g.  buchholz+c17),
are very interesting from the perspective of future simulation-based algorithm design.









§ SEQUENTIAL MONTE CARLO SAMPLING






It can be difficult to design an importance sampling density g(θ) that is able to efficiently place a large number of samples in regions of high posterior density. Sequential Monte Carlo (SMC) and sequential importance sampling (SIS) algorithms are designed to overcome this difficulty by constructing a sequence of slowly changing intermediary distributions f_m(θ), m=0,…,M, where f_0(θ)=g(θ) is the initial importance sampling distribution, and f_M(θ)=f(θ) is the target distribution of interest. A population of particles (i.e samples θ^(i), i=1,…,N) is then propagated between these distributions, in sequence, so that f_1(θ),…,f_M-1(θ) act as an efficient importance sampling bridge between g(θ) and f(θ). There are a number of techniques available for specification of the intermediary distributions (e.g. geyer+t95, delmoral06). 

There is a rich literature on the construction of efficient SMC and SIS algorithms. See e.g. liucw98, gilksb01, neal01, doucetfg01 and chopin02 among others. These algorithms invariably involve some combination of three main ideas.

Given a weighted sample (θ^(1)_m-1,w^(1)_m-1),…,(θ^(N)_m-1,w^(N)_m-1) from intermediary distribution f_m-1(θ), the reweighting step propagates the particles to the next intermediary distribution f_m(θ). This could involve a simple importance reweighting, or something more involved if hybrid importance/rejection schemes are employed (e.g. liucw98).

Depending on the efficiency of the transitions between f_m-1(θ) and f_m(θ), the variability of the importance weights w_m^(i) could be very high, with some particles having very small weights, and others having very large weights – commonly known as particle degeneracy. This can be measured through the effective sample size (<ref>) (liucw98, liu01). The resampling step is designed to replenish the particle population by resampling the particles from their empirical distribution (θ^(1)_m,w^(1)_m),…,(θ^(N)_m,w^(N)_m). In this manner, particles with low weights in regions of low density will likely be discarded in favour of particles with higher weights in regions of higher density. Following resampling, the effective sample size will be reset to N as each weight will then be set to w^(i)_m=1/N. Resampling should not occur too frequently. A common criterion is to resample when the effective sample size falls below a pre-specified threshold, typically E=N/2. See e.g. douccm05 for a review and comparison of various resampling methods. 

Finally, the move step aims to both move the particles to regions of high probability, and increase the particle diversity in the population. The latter is important since, particularly after resampling, particles with high weights can be replicated in the sample. Any transition kernel F_m(θ,θ') can be used for the move step, although an MCMC kernel is a common choice (e.g. gilksb01) as it results in the importance weight being unchanged, although there is also the chance that the proposed move is rejected. Other kernels, such as F_m(θ,θ')=ϕ(θ'; θ,σ_m^2) to add a random normal scatter to the particles, will require the importance weights to be modified. See e.g. delmoral06 for discussion on different forms of the move kernel.




 §.§ Sequential importance sampling



In the ABC framework a natural choice for the sequence of intermediary distributions is

    f_m(θ) =π_ABC,h_m(θ,s|s_obs)∝ K_h_m(s-s_obs)p(s|θ)π(θ),

for m=0,…,M, indexed by the kernel scale parameter, where the sequence h_0≥ h_1≥…≥ h_M is a monotonic decreasing sequence. Accordingly, each successive distribution with decreasing h_m, will be less diffuse and a closer approximation to π(θ|s_obs) sisson+ft07.

A sequential importance sampling version of the ABC rejection control importance sampler (Algorithm 4) is given in Algorithm 7.



This algorithm is a particular version of the sampler proposed by peters+fs12 (see also sisson+ft07) who incorporated the partial rejection control mechanism of liucw98 and liu01 into the SMC sampler framework. When applied in the ABC setting, rejection control provides one means of controlling the otherwise highly variable particle weights. 

As with the ABC rejection control importance sampler (Algorithm 4), samples from an importance sampling distribution g_m(θ), constructed from the samples from the previous population targeting f_m-1(θ), are combined with the rejection control mechanism in order to target f_m(θ).




The initial sampling distribution g(θ) can be any importance sampling density, as with standard importance sampling algorithms.  There are a number of adaptive ways to construct the subsequent importance distributions g_m(θ) for m=1,…,M, based on the population of samples from the previous intermediary distribution (θ^(1)_m-1,w^*(1)_m-1),…,(θ^(N)_m-1,w^*(N)_m-1).

The simplest of these is to  specify g_m(θ) as some standard parametric family, such as the multivariate Normal distribution, with parameters estimated from the previous particle population (e.g. chopin02).

Another option is to construct a kernel density estimate of the distribution of the previous particle population
g_m(θ) = ∑_i=1^N W^*(i)_m-1F_m(θ^(i)_m-1,θ) where W_m-1^*(i)=w_m-1^*(i)/∑_j=1^Nw_m-1^*(j), and F_m(θ,θ') is some forward mutation kernel describing the probability of moving from θ to θ', such as F_m(θ,θ') = ϕ(θ'; θ,Σ_m), the multivariate normal density function centred at θ and with covariance matrix Σ_m delmoral06,beaumont+cmr09,peters+fs12. 

These and other possibilities may also be constructed by first reweighting the draws from the previous population f_m-1(θ) so that they target f_m(θ).




If the kernel K_h has compact support then step <ref> of Algorithm 7 will automatically reject any θ^(i)_m for which K_h_m(s^(i)_m-s_obs)=0. (This also happens for Algorithm 4.) This practical outcome occurs for most ABC SIS and SMC algorithms used in practice, as use of the uniform kernel is predominant (e.g. sisson+ft07, toni+wsis09, beaumont+cmr09, delmoral+dj12), although the rejection of θ^(i) is sometimes hard coded as in Algorithm 3, rather than being part of a more sophisticated importance weight variance control mechanism, such as rejection control. 



In the limit as rejection thresholds c_m→ 0 for m=1,…,M (and defining 0/0:=1), the rejection control mechanism will allow all particles to proceed to the next stage of the algorithm. Therefore c_m→0 represents a standard sequential importance sampler that will likely result in the collapse of the particle population (i.e. all weights w^*(i)_m=0) in the ABC setting, for low h_m.

However, non-zero rejection control thresholds c_m permit a finer scale control over the importance weights w^(i)_m beyond distinguishing between zero and non-zero weights, with larger c_m resulting in more similar weights with less variability, though at the price of higher computation through more rejections. In this manner, rejection control provides one way in which ABC SMC algorithms may be implemented with kernels K_h that are non-uniform, or have non-compact support, without which the effective sample size of the sampler would deteriorate almost immediately for low h_m peters+fs12.













As with the ABC rejection control importance sampler (Algorithm 4), suitable rejection thresholds may be dynamically determined during algorithm run-time by, for each m, first implementing steps <ref> and <ref> for i=1,…,N, specifying c_m as some function (such as a quantile) of the empirical distribution of the realised w^(1)_m,…,w^(N)_m, and then continuing Algorithm 7 from step <ref> onwards for each i=1,…,N peters+fs12.







The sequence of scale parameters h_0≥ h_1≥…≥ h_M in Algorithm 7 has been presented as requiring  pre-specification in order to implement the sampler. However, as with any annealing-type algorithm, identifying an efficient sequence is a challenging problem. Fortunately, as with the automatic determination of the rejection control thresholds c_m, choice of the scale parameters can also be automated, and one such method to achieve this is discussed in the next Section. To initialise the algorithm efficiently, setting h_0=∞ would result in all particles θ_0^(1),…,θ_0^(N) having relatively similar weights w_0^(i)=π(θ^(i)_0)/g(θ^(i)_0), as a function of the prior and initial sampling distributions.




 §.§ Sequential Monte Carlo samplers




An alternative representation of population based algorithms is the sequential Monte Carlo sampler delmoral06. Here, the particles are defined on the space of the path that each particle will take through the sequence of distributions f_0(θ),…,f_M(θ).  Hence, if θ^(i)_m∈Θ, then the path of particle i through the first m distributions is given by θ_1:m^(i) = (θ_1^(i),…,θ_m^(i))∈Θ^m for m=1,…,M. SMC samplers explicitly implement each of the reweighting, resampling and move steps, and at their most general level have sophisticated implementations (e.g. delmoral06). A number of SMC samplers have been developed in the ABC framework (see sisson+ft07, toni+wsis09, beaumont+cmr09, drovandi2011estimation, delmoral+dj12). 

Algorithm 8 presents a generalisation (to general kernels K_h) of the adaptive ABC SMC sampler of delmoral+dj12. 


This algorithm provides an alternative method to rejection control  to avoid the collapse of the particle population, for an arbitrary choice of kernel K_h, by making particular sampler design choices.

Firstly, the probability of generating particles θ_m^(i) with identically zero weights w^(i)_m=0 is reduced by

increasing the number of summary statistics drawn to T, thereby targeting the joint distribution π_ABC(θ,s(1),…,s(T)|s_obs) as described in Section <ref>, although at the price of greater computation. Within the scope of an ABC SMC sampler that makes use of MCMC kernels within the move step (as with Algorithm 8), the alternative algorithms analysed by lee+l14 (see Section <ref>) could also be implemented (e.g. bernron+jgr17).


In combination with the increased number of summary statistic replicates, Algorithm 8 directly controls the degree of particle degeneracy in moving from distribution f_m-1(θ) to f_m(θ). In particular, the next kernel scale parameter h_m<h_m-1 is chosen as the value which results in the effective sample size following the reweighting step, being reduced by a user specified proportion, α. In this manner, the sample degeneracy will reduce in a controlled manner at each iteration, and the sequence of h_m will adaptively reduce at exactly the rate needed to achieve this. When the effective sample size is reduced below some value E, resampling occurs and resets the effective sample size back to N, and the process repeats. As a result, resampling repeatedly occurs automatically after a fixed number of reweighting steps, as determined by α.

This algorithm requires a stopping rule to terminate. If left to continue, h_m would eventually  reduce very slowly, which is an indication that the sampler can no longer efficiently move the particles around the parameter space. delmoral+dj12 argue that this identifies natural values of h_M that should then be adopted. In particular, they terminate their algorithm when the MCMC move rate drops below 1.5%, which then determines the final value of h_M. Alternative strategies to adaptively choose the kernel scale parameter sequence have been proposed by drovandi2011estimation,  drovandi2011likelihood, silk+fs13 and daly+ghc17.


SMC algorithms provide many easy opportunities for sampler adaptation, unlike MCMC samplers which are constrained by the need to maintain the target distribution of the chain. For example, within ABC SMC algorithms, 

prangle17 adaptively learns the relative weightings of the summary statistics within the distance function s-s_obs to improve efficiency, 

bonassi+w15 construct adaptive move proposal kernels g_m(θ_m^(i),θ)=∑_i=1^Nν^(i)_m-1F_m(θ_m-1^(i),θ) based on weighting components of g_m, via ν^(i)_m-1, based on the proximity of s^(i)_m-1 to s_obs,

and filippi+bcs13 develop a different method of adaptively constructing the sequence of intermediary distributions, f_m(θ),  based on Kullback-Leibler divergences between successive distributions.


Other ideas can be incorporated within ABC SMC algorithms in particular settings, or can use the ideas from ABC SMC algorithms to tackle problems related to posterior simulation.
For example, prangle+ek17 use ideas from rare event modelling to improve sampler efficiency within ABC SMC algorithms.

When simulation from the model p(s|θ) is expensive, everitt+r17 first use a cheap approximate simulator within an ABC SMC algorithm to rule out unlikely areas of the parameter space, so that expensive computation with the full simulator is avoided until absolutely necessary.

jasra+smm12 implement an ABC approximation within an SMC algorithm to perform filtering for a hidden Markov model.

dean+sjp14 and yildirim+sdj15 use ABC SMC methods for optimisation purposes (with a different sequence of intermediary distributions), so as to derive maximum (intractable) likelihood estimators for hidden Markov models.









§ DISCUSSION







ABC samplers have proved to be highly accessible and simple to implement, and it is this that has driven the popularity and spread of ABC methods more generally.

Multi-model versions of each of these algorithms are available (e.g. toni+wsis09, chkerbtii+ccb15) or can be easily constructed, with ABC posterior model probabilities and Bayes factors being determined by the relative values of 1/N∑_i=1^NK_h(s^(i)-s_obs) under each model.

Although here the user needs to clearly understand the ideas of summary statistic informativeness for model choice (e.g. marin+prr14) and the problems involved in computing Bayes factors as h→ 0 martin+frr17.



Improvements to general ABC samplers include increasing algorithmic efficiency by using quasi Monte Carlo methods <cit.>, and the use of multi-level rejection sampling warne+bs17 (see jasra+jnst17 for the SMC version) for variance reduction.

The lazy ABC method of prangle16 states that it may be possible to terminate expensive simulations s∼ p(s|θ) early, if it is also possible to calculate the probability that the full simulation when run to completion would have been rejected.

Diagnostics to determine whether the kernel scale parameter h_m is sufficiently low that π_ABC(θ|s_obs) is indistinguishable from π(θ|s_obs) were developed by prangle+bps14.

There are many related results on the rate of convergence of ABC algorithms as measured through the mean squared error of point estimates blum10,fearnhead+p12,calvet+c15,biau+cg15,barber+vm15.

ABC samplers have also allowed previously unclear links to other algorithms to become better understood – for example, nott+mt12 have reinterpreted the Kalman filter as an ABC algorithm, and drovandi18 (this volume) has comprehensively described the links between ABC and indirect inference.



A number of algorithms related to ABC methods have emerged, including 
Bayesian empirical likelihoods mengersen+pr13 and bootstrap likelihoods zhu+ml16, the synthetic likelihood (wood10, drovandi+gkr18, this volume), the expectation-propagation ABC algorithm (barthelmec14,barthelme+cc17, this volume), 
albert+ks15's particle-based simulated annealing algorithm, 

and forneron+n16's
 optimisation-based likelihood free importance sampling algorithm.

Perhaps the biggest offshoot of ABC samplers is the more general pseudo-marginal Monte Carlo method <cit.>, which implements exact Monte Carlo simulation with an unbiased estimate of the target distribution, of which ABC is a particular case. See andrieu+sv17 (this volume) for an ABC-centred exploration of these methods.







§ ACKNOWLEDGEMENTS


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

