New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Sequential Monte Carlo without likelihoods

Edited by Michael S. Waterman, University of Southern California, Los Angeles, CA, and approved December 4, 2006 (received for review August 19, 2006)
This article has a Correction. Please see:
 Correction for Sisson et al., Sequential Monte Carlo without likelihoods  September 02, 2009
Abstract
Recent new methods in Bayesian simulation have provided ways of evaluating posterior distributions in the presence of analytically or computationally intractable likelihood functions. Despite representing a substantial methodological advance, existing methods based on rejection sampling or Markov chain Monte Carlo can be highly inefficient and accordingly require far more iterations than may be practical to implement. Here we propose a sequential Monte Carlo sampler that convincingly overcomes these inefficiencies. We demonstrate its implementation through an epidemiological study of the transmission rate of tuberculosis.
 approximate Bayesian computation
 Bayesian inference
 importance sampling
 intractable likelihoods
 tuberculosis
Termed approximate Bayesian computation (ABC), recent new methods in Bayesian inference have provided ways of evaluating posterior distributions when the likelihood function is analytically or computationally intractable (1–4). ABC algorithms represent a substantial methodological advance because they now admit realistic inference on problems that were considered intractable only a few years ago. The rapidly increasing use of these methods has found application in a diverse range of fields, including molecular genetics (5), ecology (6), epidemiology (7), evolutionary biology (8, 9) and extreme value theory (1).
Given a likelihood function, f(x_{0}θ), and a prior distribution π(θ) on the parameter space, Θ, interest is in the posterior distribution f(θx_{0}) ∝ f(x_{0}θ)π(θ), the probability distribution of the parameters having observed the data, x_{0} (10, 11).
To avoid directly evaluating the likelihood, all ABC algorithms incorporate the following procedure to obtain a random sample from the posterior distribution. For a candidate parameter vector θ* drawn from some density, a simulated data set x* is generated from the likelihood function f(xθ*) conditioned on θ*. This vector is then accepted if simulated and observed data are sufficiently “close.” Here, closeness is achieved if a vector of summary statistics S(·) calculated for the simulated and observed data are within a fixed tolerance (ϵ) of each other according to a distance function ρ (e.g., Euclidean distance). In this manner, ABC methods sample from the joint distribution f(θ, xρ(S(x), S(x_{0})) ≤ ϵ), where interest is usually in the marginal f(θρ(S(x), S(x_{0})) ≤ ϵ). The algorithms work by accepting a value θ with an average probability of Pr(ρ(S(x), S(x_{0})) ≤ ϵθ). If the summary statistics S(·) are nearsufficient and ϵ is small then f(θρ(S(x), S(x_{0})) ≤ ϵ) should be a reasonable approximation to f(θx_{0}).
Existing ABC methods for obtaining samples from the posterior distribution either involve rejection sampling (3, 4, 12) or Markov chain Monte Carlo (MCMC) simulation (1, 2). Both of these classes of methods can be inefficient. The ABC rejection sampler proceeds as follows
ABCREJ Algorithm
REJ1. Generate a candidate value θ* ∼ π(θ) from the prior.

REJ2. Generate a data set x* ∼ f(xθ*).

REJ3. Accept θ* if ρ(S(x*), S(x_{0})) ≤ ϵ.

REJ4. If rejected, go to REJ1.
Each accepted vector represents an independent draw from f(θρ(S(x), S(x_{0})) ≤ ϵ). Acceptance rates for algorithm ABCREJ can be very low as candidate parameter vectors are generated from the prior π(θ), which may be diffuse with respect to the posterior. Accordingly, Marjoram et al. (2) proposed to embed the likelihoodfree simulation method within the well known MCMC framework. This algorithm proceeds as follows
ABCMCMC Algorithm
MC1. Initialize θ_{1}, i = 1.

MC2. Generate a candidate value θ* ∼ q(θθ_{i}), where q is some proposal density.

MC3. Generate a data set x* ∼ f(xθ*).

MC4. Set θ_{i}_{+1} = θ* with probability otherwise set θ_{i+1} = θ_{i}.

MC5. If i < N, increment i = i + 1 and go to MC2.
Here 1(A) = 1 if A is true, and 0 otherwise. The candidate vector is generated from an arbitrary proposal density q(θ·) and accepted with the usual Metropolis–Hastings acceptance probability. The (intractable) likelihood ratio is now coarsely approximated by 1 if simulated and observed data are sufficiently “close,” and 0 otherwise. Algorithm ABCMCMC generates a sequence of serially and highly correlated samples from f(θρ(S(x), S(x_{0})) ≤ ϵ). Determination of the chain length, N, is therefore obtained through a careful assessment of convergence (13) and consideration of the chain's ability to explore parameter space (i.e., chain mixing).
When the prior and posterior are dissimilar, algorithm ABCMCMC delivers substantial increases in acceptance rates over algorithm ABCREJ [Marjoram et al. (2) report 0.2% acceptance rates over 0.0008% in a simple coalescent tree analysis], although at the price of generating dependent samples. However, because acceptance rates for ABC samplers are directly proportional to the likelihood; if the ABCMCMC sampler enters an area of relatively low probability with a poor proposal mechanism, the efficiency of the algorithm is strongly reduced because it then becomes difficult to move anywhere with a reasonable chance of acceptance, and so the sampler “sticks” in that part of the state space for long periods of time. This is illustrated in the following toy example.
Toy Example
As an illustration, suppose that the posterior of interest is given by the mixture of two normal distributions where φ(μ, σ^{2}) is the density function of a N(μ, σ^{2}) random variable. Here, the second component implies large regions of relatively low probability with respect to the lower variance first component (Figure 1Lower). In the ABC setting, this posterior may be realized by drawing x = (x_{1},…, x_{100}), x_{i} ∼ N(θ, 1) and by specifying where x̄ = 1/100 Σ x_{i} denotes the sample mean.
With a prior of π(θ) ∼ U(−10, 10) and ϵ = 0.025, the ABCREJ algorithm required a mean of 400.806 datageneration steps (simulations from the likelihood) for each accepted realization based on 1,000 accepted realizations. With an acceptance rate of <0.25%, this is highly inefficient.
In contrast, Fig. 1 shows the result of implementing the ABCMCMC algorithm initialized at θ_{0} = 0, with N = 10,000 iterations and with proposals generated via the random walk q(θθ_{t}) ∼ N(θ_{t}, 0.15^{2}). When the sampler is within the highdensity region, transitions between different parameter values are frequent (acceptance rate ≈5%). However, when the sampler moves outside of this region, the frequency of transitions drops markedly, especially so for the extended period at ≈5,000–9,000 iterations. Of course, the samples should visit the distributional tails and will do so for a number of iterations proportional to the posterior probability. However, as is evident from the histogram of the sampler output, the number of further iterations required so that the realizations in the tail are in proportion to the target distribution will be far in excess of the initial 10,000.
In the search for improved ABC methods, poor Markov chain performance may be improved by embedding the target distribution as a marginal distribution within a larger family of related distributions among which it is far easier for the sampler to move around. This was essentially the approach adapted for ABC by Bortot et al. (1), although such algorithms are wasteful by construction in that samples from the auxiliary distributions are not used in the final analysis.
As an alternative, we propose to improve upon simple rejection sampling by adopting a sequential Monte Carlo (SMC)based simulation approach. Here, a full population of parameters θ^{(1)},…, θ^{(N)} (termed particles) is propagated from an initial, userspecified distribution, through a sequence of intermediary distributions, until it ultimately represents a sample from the target distribution. SMC methods can be considered an extension of importance sampling. We will demonstrate that the SMC approach can substantially outperform both MCMC and rejection sampling in the ABC framework.
Advantages of the SMC approach are

Like rejection sampling, the sampler will never become “stuck” in regions of low probability.

Unlike rejection sampling, severe inefficiencies generated by mismatch of (initial) sampling and target distributions are avoided.

Particles that represent the target posterior poorly are eliminated in favor of those particles that represent the posterior well.

The populationbased nature of the sampler means that complex (e.g., multimodal) posteriors may be explored more efficiently.

Samples are obtained from a number of distributions with differing tolerances (ϵ). This permits an a posteriori, or dynamic, examination of the robustness of the posterior to this choice.

Unlike MCMC, SMC particles are uncorrelated and do not require the determination of a burnin period or assessment of convergence.
Disadvantages of the SMC approach are considered in Discussion.
Here, we propose an ABC sampler based on SMC simulation. We initially outline a generic SMC algorithm before exploiting ideas based on PRC to derive a more efficient algorithm for the ABC setting. Finally, we demonstrate its utility with regard to the toy example considered above and in a reexamination of a previously implemented analysis of the transmission rate of tuberculosis.
Methods
SMC Without Likelihoods.
We wish to sample N particles θ^{(1)},…, θ^{(N)} from the posterior f(θρ(S(x), S(x_{0})) ≤ ϵ), for observed data x_{0}, and for unknown parameter vector θ ∈ Θ. We assume that the summary statistics S(·), the tolerance ϵ, and the distance function ρ are known.
Let θ _{1}^{(1)},…, θ_{1}^{(N)} ∼ μ_{1}(θ) be an initial population from a known density, μ_{1}, from which direct sampling is possible, and by f_{T}(θ) = f(θρ(S(x), S(x_{0})) ≤ ϵ) the target distribution (this notation will become clear shortly). Standard importance sampling would then indicate how well each particle θ_{1}^{(i)} adheres to f_{T}(θ) by specifying the importance weight, W_{T}^{(i)} = f_{T}(θ _{1}^{(i)})/μ_{1}(θ_{1}^{(i)}), it should receive in the full population of N particles. The effectiveness of importance sampling is sensitive to the choice of sampling distribution, μ_{1}. The prior π(θ) is often used for this purpose. Importance sampling can be highly inefficient if μ_{1} is diffuse relative to f_{T} and can fail completely in the case of sampling and target distribution mismatch.
The idea behind sequential sampling methods is to avoid the potential disparity between μ_{1} and f_{T} by specifying a sequence of intermediary distributions f_{1},…, f_{T−1}, such that they evolve gradually from the initial distribution towards the target distribution. For example, one can choose a geometric path specification where f_{t}(θ) = f_{T}(θ)^{φt} μ_{1}(θ)^{1−φt} with 0 < φ_{1} < … < φ_{T} = 1 (14, 15). Hence, it is possible to move smoothly and effectively in sequence from μ_{1} to f_{T} using repeated importance sampling, generating a series of particle populations {θ_{t}^{(i)}} = {θ_{t}^{(i)} : i = 1,…, N}, for t = 1, … T. That is, sequential methods proceed by moving and reweighting the particles according to how well they adhere to each successive distribution, f_{t}.
In the ABC setting, we may naturally define the sequence of distributions f_{1}…, f_{T} as Here, _{x(1)},…, x_{(Bt)} are B_{t} data sets generated under a fixed parameter vector, x_{(b)} ∼ f(xθ), and {ϵ_{t}} is a strictly decreasing sequence of tolerances. The nested family of distributions generated by varying ϵ (continuously) was considered by Bortot et al. (1) in their augmented state space ABC algorithm. By specifying ϵ_{t} < ϵ_{t−1}, we ensure that the likely range of parameter values in each progressive distribution is a subset of the one it precedes. This is a desirable property for our sampling distributions. By specifying ϵ_{T} = ϵ, we realize the final particle population {θ_{T}^{(i)}} as a weighted sample from the target distribution. Setting B_{t} = B and ϵ_{t} = ϵ for all t reduces to the likelihood specified by Marjoram et al. (2), and further B = 1 to the likelihood adopted in algorithm ABCMCMC. The target distribution, f_{T}, is specified by ϵ_{T} = ϵ.
PRC.
In ABC samplers, B_{t} = 1 is the most computationally efficient specification, in that some action occurs (e.g., a realization or move proposal is accepted) each time a nonzero likelihood is generated. However, because the particle weight under SMC methods is directly proportional to the likelihood, there is a large probability that the likelihood, and therefore the particle weight, will be zero, thereby rendering the particle useless. Fortunately, the idea of PRC (see chapters 2 and 3 of ref. 16) permits the repeated resampling (and moving) of particles from the previous population to replace those particles with zero weight. PRC continues until N particles with nonzero weight are obtained. See Appendix for further details.
At each step, SMC methods move each particle according to a Markov kernel K_{t} to improve particle dispersion. This induces a particle approximation to the importance sampling distribution μ_{t}(θ_{t}) = ∫_{Θ} μ_{t−1}(θ_{t−1})K_{t}(θ_{t}θ_{t−1})dθ_{t−1}, for populations t = 2,…, T. Choices of K_{t} include a standard smoothing kernel (e.g., Gaussian) or a Metropolis–Hastings accept/reject step. The kernel accordingly enters the particle weight calculation.
A recent innovation in SMC methods has been the introduction of a backward Markov kernel L_{t−1} with density L_{t−1}(θ_{t−1}θ_{t}) into the weight calculation (17). The backward kernel relates to a timereversed SMC sampler with the same target marginal distribution as the (forwardtime) SMC sampler induced by K_{t}. Because only specification of K_{t} is required in order to implement an SMC sampler, the backward kernel is essentially arbitrary. The kernel L_{t−1} may therefore be optimized to minimize the variance of the weights induced by the importance sampling distribution μ_{t} (through K_{t}). This is difficult in general, however, so simplified forms are often adopted. See ref. 17 for further discussion.
SMC algorithms measure the degree of sample degeneracy within each population through the effective sample size (ESS). ESS calculates the equivalent number of random samples required to obtain an estimate, such that its Monte Carlo variation is equal to that of the N weighted particles. This may be estimated as 1 ≤ [Σ_{i=1}^{N} (W_{t}^{(i)})^{2}]^{−1} ≤ N for each t (16, 18). Sample degeneracy can occur through sampling and target distribution mismatch when a small number of particles have very large weights. Through a resampling step, particles with larger weights become better represented in the resampled population than those with smaller weights. Those particles with sufficiently small weights, which poorly approximate f_{t}, may be eliminated. The resampling threshold, E, is commonly taken to be N/2.
Combining PRC with SMC, we obtain the following (ABCPRC) algorithm
ABCPRC Algorithm
PRC1. Initialize ϵ_{1},…, ϵ_{T}, and specify initial sampling distribution μ_{1}.
Set population indicator t = 1.

PRC2. Set particle indicator i = 1.

PRC2.1. If t = 1, sample θ** ∼ μ_{1}(θ) independently from μ_{1}.
If t > 1, sample θ* from the previous population {θ_{t−1}^{(i)}} with weights {W_{t−1}^{(i)}}, and perturb the particle to θ** ∼

K_{t}(θθ*) according to a Markov transition kernel K_{t}. Generate a data set x** ∼ f(xθ**).
If ρ(S(x**), S(x_{0})) ≥ ϵ_{t}, then go to PRC2.1.

PRC2.2. Set where π(θ) denotes the prior distribution for θ, and L_{t−1} is a backward transition kernel.

If i < N, increment i = i + 1 and go to PRC2.1.

PRC3. Normalize the weights so that Σ_{i=1}^{N} W_{t}^{(i)} = 1.
If ESS = [Σ_{i=1}^{N}(W_{t}^{(i)})^{2}]^{−1} < E then resample with replacement, the particles {θ_{t}^{(i)}} with weights {W_{t}^{(i)}} to obtain a new population {θ_{t}^{(i)}}, and set weights {W_{t}^{(i)} = 1/N}.

PRC4. If t < T, increment t = t + 1 and go to PRC2.
Samples {θ_{T}^{(i)}} are weighted samples from the posterior distribution f(θρ(S(x), S(x_{0})) ≤ ϵ). The validity of this algorithm is derived by construction from the validity of the combination of general SMC methods and the PRC process (see Appendix). Algorithm ABCPRC corresponds to algorithm ABCREJ for the special case when T = 1 and μ_{1}(θ) = π(θ).
For the remainder of this article, we consider K_{t}(θ_{t}θ_{t−1}) = L_{t−1}(θ_{t−1}θ_{t}) as a Gaussian kernel with common variance (following ref. 19), which we have found to perform adequately. For discussions on closer to optimal choices of L_{t−1}, see ref. 17, and for applications of SMC and more technical proofs of the SMC algorithm's validity, see refs. 16, 17, and 20–24. Finally, we note that if K_{t}(θ_{t}θ_{t−1}) = L_{t−1}(θ_{t−1}θ_{t}), μ_{1}(θ) = π(θ) and the prior π(θ) ∝ 1 over Θ, then all weights are equal throughout the sampling process and may safely be ignored [in addition to ignoring all population (PRC3) resampling steps].
Results
Toy Example (Revisited).
We now implement algorithm ABCPRC in the mixture of normals posterior considered earlier. We generate a sample of N = 1,000 particles by considering a sequence of three distributions f_{1}, f_{2}, and f_{3} defined by Eq. 1 with ϵ_{1} = 2, ϵ_{2} = 0.5 and ϵ_{3} = 0.025, and with prior distribution π(θ) ∼ U(−10, 10). We specify μ_{1}(θ) = π(θ) and K_{t}(θ_{t}θ_{t−1}) = L_{t−1}(θ_{t−1}θ_{t}) as a Gaussian random walk so that all weights are equal.
The initial (μ_{1}) population and histograms of f_{1} to f_{3} are illustrated in Figure 2. The movement in distribution towards the target distribution is a clear progression, with the final sample adhering remarkably well to the target distribution, especially in the tails where the ABCMCMC algorithm performed particularly poorly (Figure 1).
ABC algorithms may be intuitively compared through the number of likelihood “evaluations,” that is, the number of datageneration steps. Table 1 enumerates the mean number of datageneration steps required to move a particle between two successive populations. As the tolerance reduces with each successive distribution, f_{t}, the number of datageneration steps we expect increases. This effect is partially offset by the degree of similarity between population distribution f_{t} and its induced sampling distribution μ_{t}. The total number of datasimulation steps in the ABCPRC algorithm was 75,895. This is more than the illustrated 10,000 for the ABCMCMC algorithm (Figure 1), but this latter simulation requires a substantially longer run before we can be satisfied that a representative sample has been drawn. Accordingly, the ABCPRC algorithm is far more efficient for this case.
In contrast, using the ABCREJ algorithm results in utilizing >400 datasimulation steps per final particle (Table 1). Here, there is a clear advantage in adopting a series of intermediary distributions between μ_{1} and the target distribution. Finally, as an indication of the maximum possible efficiency of ABC samplers for this example, performing rejection sampling with sampling distribution equal to the posterior distribution requires >21 datageneration steps per final particle. Note that each particle must still satisfy steps REJ2 and REJ3, so we do not automatically accept every particle proposed. This gives algorithm ABCPRC 28% of maximum possible efficiency in this case, compared with only 5% for rejection sampling.
Analysis of Tuberculosis Transmission Rates.
We now reimplement an analysis of tuberculosis transmission rates originally investigated using algorithm ABCMCMC (7). The aim of this study was to estimate three compound parameters of biological interest, namely, the reproductive value, the doubling time, and the net transmission rate. The data used come from a study of tuberculosis isolates collected in San Francisco during the early 1990s by Small et al. (25). These consist of 473 isolates genotypes using the IS6110 marker; the resulting DNA fingerprints can be grouped into 326 distinct genotypes as follows where n^{k} indicates there were k clusters of size n. The ABCMCMC method was used in conjunction with a stochastic model of transmission, which is an extension of a birth and death process to include mutation of the marker. Simulating samples from this model allows comparisons with the actual data through two summary statistics: g, the number of distinct genotypes in the sample, and H, the gene diversity. An informative prior was used for the mutation rate taken from published estimates of the rate for IS6110. Further details can be found in ref. 7.
Tanaka et al. (7) previously implemented the ABCMCMC algorithm with tolerance ϵ = 0.0025. Three Markov chains with an average acceptance rate of ≈0.3% were thinned and combined to form the final sample, utilizing >2.5 million datageneration steps. (Actually, more were used, as one chain became “stuck” in a distributional tail for most of its length, as illustrated in Fig. 1, and had to be rerun.)
We illustrate algorithm ABCPRC with a sequence of 10 distributions, defined by ϵ_{1} = 1 and for t = 2,…, 9, ϵ_{t} = ½(3ϵ_{t−1} − ϵ_{10}) is taken to be halfway between the previous tolerance and the target of ϵ_{10} = 0.0025. Ten distributions were selected so that successive distributions were reasonably similar. We adopt K_{t}(θ_{t}θ_{t−1}) = L_{t−1}(θ_{t−1}θ_{t}) as the ABCMCMC Gaussian random walk proposal of Tanaka et al. (7) with a slightly larger step for the mutation parameter.
Based on a population size of N = 1,000, Fig. 3 illustrates smoothed posterior distributions of the quantities of interest: (Figure 3 Left) the net transmission rate (α − δ) is the rate of increase of the number of infectious cases in the population; the doubling time [log(2)/(α − δ)] is the required duration for the number of infectious cases in the population to double; (Figure 3 Right) the reproductive value (α/δ) is the expected number of new infectious cases produced by a single infectious case while the primary case is still infectious. As is evident, the results of the ABCMCMC and ABCPRC algorithms are essentially indistinguishable.
Relative algorithm efficiencies can again be measured by the mean number of datageneration steps per final particle. Table 2 lists the number of datageneration instances in ABCPRC and ABCREJ algorithms. For algorithm ABCREJ, this amounts to a mean of 7,206.3 datageneration steps per particle. In contrast, algorithm ABCPRC yields a mean of 1,421.3 datageneration steps per particle, >5 times more efficient.
Comparisons to the original ABCMCMC analysis of Tanaka et al. (7) can also be made in terms of the number of data generation steps required to generate one uncorrelated particle. Here, the Markov nature of the sampler and the very low acceptance rates induce a strongly dependent sequence. Thinning this sequence so that there were no significant (marginal) autocorrelations above lag 10 resulted in using 8,834 data generations steps per realization. Repeating this so that there were no significant autocorrelations at any lag yielded 67 uncorrelated particles, corresponding to ≈27,313 datageneration steps per final realization. By this measure, algorithm ABCPRC is ≈20 times more efficient than the MCMC implementation.
Discussion
Likelihoodfree samplers for Bayesian computation are growing in importance, particularly in population genetics and epidemiology, so it is crucial that efficient and accessible algorithms are made available to the practitioner. Existing MCMC algorithms exhibit an inherent inefficiency in their construction, whereas rejection methods are wasteful when sampling and target distributions are mismatched. Through an SMC approach, we may circumvent these problems and generate improved sampling, particularly in distributional tails, while achieving large computational savings.
Evaluations of certain userspecified aspects of algorithm ABCPRC have not been presented, although these have been studied for SMC algorithms in general, and the necessity of their specification could be considered a disadvantage of the SMC approach. The incorporation of measures other than effective sample size to determine the optimal resampling time is given by Chen et al. (26) and forward and backward kernel choice by Del Moral et al. (17). Jasra et al. (27) give a study of various tolerance schedules and the number of distributions, f_{1},…, f_{T}. It seems credible that the tolerance schedule and distribution number could be determined dynamically based on onestepahead estimates of distributional change, f_{t−1} → f_{t}, and the required computation (number of datageneration steps). This could be considered one method of selecting the final tolerance ϵ_{T}.
Acknowledgments
We thank two anonymous referees whose suggestions have strongly improved the final form of this article. This work was supported by the Australian Research Council through the Discovery Grant scheme (DP0664970 and DP0556732) and by the Faculty of Science, University of New South Wales.
Appendix
We briefly justify the use of partial rejection control in deriving algorithm ABCPRC. Following ref. 17, a generic sequential Monte Carlo algorithm is implemented as follows
SMC Algorithm
SMC1. Identify the sequence of distributions f_{1},…, f_{T}, where f_{T} corresponds to the target distribution, and initial sampling distribution μ_{1}.
Set population indicator t = 1.

SMC2. Set particle indicator i = 1.

SMC2.1. If t = 1, sample θ_{t}^{(i)} ∼ μ_{1}(θ) independently from μ_{1}.
If t > 1, perturb each particle to θ_{t}^{(i)} ∼ K_{t}(θθ_{t−1}^{(i)}) according to a Markov transition kernel K_{t}.

SMC2.2. Evaluate weights W_{t}^{(i)} for each particle according to where f_{t} denotes the intermediate distribution at step t and L_{t−1} is a backward transition kernel.
If i < N, increment i = i + 1 and go to SMC2.1.

SMC3. Normalize the weights so that Σ_{i=1}^{N} W_{t}^{(i)} = 1.
If ESS = Σ_{i=1}^{N} (W_{t}^{(i)})^{2}]^{−1} < E, then resample with replacement, the particles {θ_{t}^{(i)}} with weights {W_{t}^{(i)}} to obtain a new population {θ_{t}^{(i)}}, and set weights {W_{t}^{(i)} = 1/N}.

SMC4. If t < T, increment t = t + 1 and go to SMC2.
Algorithm SMC can be justified intuitively by considering the final weight of particle θ_{T}^{(i)}, assuming no weight normalization for clarity The ratio f_{T}(θ_{T}^{(i)})/μ_{1}(θ_{1}^{(i)} is immediately identifiable as the standard importance sampling weight with μ_{1} as the sampling distribution. The product of kernel ratios term evaluates the ratio of probabilities of moving from θ_{T}^{(i)} → θ_{1}^{(i)} (numerator) and from θ_{1}^{(i)} → θ_{T}^{(i)} (denominator).
Recognizing that many particles of small weight will have minimal impact on the final population, f_{T}, (e.g., they may be lost in the resampling step SMC3), partial rejection control aims to remove them at an earlier stage according to the following scheme (see chapters 2 and 3 of ref. 16).
Given a particle population {θ_{t}^{(i)}} with weights {W_{t}^{(i)}}, a small threshold, c, is specified such that all particles with weights greater than c remain unchanged. For those particles with weights smaller than c, there is a chance [with probability min (1, W_{t}^{(i)}/c)] that these particles also remain unchanged, otherwise they are replaced by a particle from the previous population {θ_{t−1}^{(i)}} chosen according to weights {W_{t−1}^{(i)}}. This particle is then propagated from distribution f_{t−1} to f_{t} (via K_{t}) as before, where its weight is then compared to the threshold, c, once more. This process is repeated until all particles have passed the threshold. PRC is performed within SMC algorithms before the population resampling step (SMC3). See ref. 16 for a justification of this approach.
For a particle population {θ_{t}^{(i)}} with weights {W_{t}^{(i)}} and t > 1, this process is given in algorithmic form by
PRC Algorithm
A1. Set threshold value c > 0 and particle indicator i = 1.

A2. With probability min{1, W_{t}^{(i)}/c}, set weight W_{t}^{(i)} = max{W_{t}^{(i)}, c} and go to A4.
Otherwise, go to A3.

A3. Sample a new particle, θ*, from {θ_{t−1}^{(i)}} with probability proportional to {W_{t−1}^{(i)}}.
Perturb the particle to θ_{t}^{(i)} ∼ K_{t}(θθ*) according to a Markov transition kernel K_{t}.
Set where W̄_{t−1} = 1/N Σ_{j=1}^{N} W_{t−1}^{(j)} and L_{t−1} is a backward transition kernel.
Go to A2.

A4. If i < N, increment i = i + 1 and go to A2.

A5. Normalize weights according to Σ_{i–1}^{N} W_{t}^{(i)} = 1.
PRC may benefit the SMC algorithm in the ABC setting as follows: The minimum computational specification for the sequence {B_{t}} in the posterior (Eq. 1) is B_{t} = 1 for all t. In this setting, large numbers of particles will have identically zero weight, as ρ(S(x), S(x_{0})) > ϵ_{t} occurs with high probability. Suppose we then implement the PRC algorithm for some c > 0 such that only identically zero weights are smaller than c. This will remove those particles for which ρ(S(x), S(x_{0})) > ϵ_{t} and replace them with particles for which ρ(S(x), S(x_{0})) ≤ ϵ_{t}, which then belong to f_{t}.
This process is equivalent to deriving the entire population {θ_{t}^{(i)}} one particle at a time, by taking random draws from the previous population, perturbing the particle, and accepting the particle if ρ(S(x), S(x_{0})) ≤ ϵ_{t}. That is, incorporating the datageneration process, we can replace step SMC2 in algorithm SMC above with step PRC2 in algorithm ABCPRC. Accordingly, we are able to maximize algorithm efficiency in that we obtain a new particle on each occasion for which ρ(S(x), S(x_{0})) ≤ ϵ_{t}, rather than wasting extra datageneration steps in overevaluating likelihood values.
Footnotes
 ↵^{‡}To whom correspondence should be addressed. Email: scott.sisson{at}unsw.edu.au

Author contributions: S.A.S. designed research; S.A.S. and Y.F. performed research; S.A.S. and M.M.T. contributed new reagents/analytic tools; Y.F. and M.M.T. analyzed data; and S.A.S., Y.F., and M.M.T. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS direct submission.
Abbreviations
 ABC,
 approximate Bayesian computation;
 MCMC,
 Markov chain Monte Carlo;
 PRC,
 partial rejection control;
 SMC,
 sequential Monte Carlo;
 ESS,
 effective sample size.
 Received August 19, 2006.
 © 2007 by The National Academy of Sciences of the USA
References
 ↵
 Bortot P,
 Coles SG,
 Sisson SA
 ↵
 Marjoram P,
 Molitor J,
 Plagnol V,
 Tavaré S
 ↵
 Beaumont MA,
 Zhang W,
 Balding DJ
 ↵
 ↵
 ↵
 Butler A,
 Glasbey C,
 Allcroft A,
 Wanless S
 ↵
 Tanaka MM,
 Francis AR,
 Luciani F,
 Sisson SA
 ↵
 Leman SC,
 Chen Y,
 Stajich JE,
 Noor MAF,
 Uyenoyama MK
 ↵
 Thornton K,
 Andolfatto P
 ↵
 Robert C,
 Casella G
 ↵
 Gilks WR,
 Richardson S,
 Spiegelhalter DJ
 ↵
 Tavaré S,
 Balding DJ,
 Griffiths RC,
 Donnelly P
 ↵
 ↵
 ↵
 ↵
 Liu JS
 ↵
 ↵
 ↵
 ↵
 Doucet A,
 de Freitas N,
 Gordon N
 ↵
 Del Moral P
 ↵
 ↵
 ↵
 Peters GW
 ↵
 ↵
 ↵
 Jasra A,
 Stephens DA,
 Holmes CC