{'fraction_non_alphanumeric': 0.04870434649726024, 'fraction_numerical': 0.011111838801929836, 'mean_word_length': 4.395524146054181, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multifidelity Approximate Bayesian Computation with Sequential Monte Carlo Parameter Sampling

Multifidelity approximate Bayesian computation (MF-ABC) is a likelihood-free technique for parameter inference that exploits model approximations to significantly increase the speed of ABC algorithms (Prescott and Baker, 2020). Previous work has considered MF-ABC only in the context of rejection sampling, which does not explore parameter space particularly efficiently. In this work, we integrate the multifidelity approach with the ABC sequential Monte Carlo (ABC-SMC) algorithm into a new MF-ABC-SMC algorithm. We show that the improvements generated by each of ABC-SMC and MF-ABC to the efficiency of generating Monte Carlo samples and estimates from the ABC posterior are amplified when the two techniques are used together.

Rather than calculating the likelihood of the observed data for a given parameter value, it is estimated. The model is simulated under a particular parameter value, and the simulated output is then compared with the observed data. If the simulation and observations are suitably close (according to a predetermined metric and threshold [11,14]) then, in the classical rejection sampling approach, the likelihood is approximated as 1. Otherwise, the likelihood is approximated as 0. This binary approximation is usually interpreted as an acceptance or a rejection of the parameter value input into the simulation. Using this approach, a weighted Monte Carlo sample from an approximated posterior can be built by repeatedly proposing parameters from the prior distribution, simulating the model with each parameter, and calculating the binary weight.
One widely-acknowledged weakness of the ABC approach is its heavy reliance on repeated simulation. There has been a significant amount of work dedicated to overcoming this reliance by exploiting methods for intelligent exploration of parameter space in order to reduce the number of simulations in areas of low likelihood [30]. For example, the parameters to be input into the simulation can instead be proposed from a importance distribution, rather than the prior, which we aim to construct in order to improve the algorithm's performance. One successful approach to importance sampling is known as Sequential Monte Carlo ABC (ABC-SMC), which aims to build consecutive samples using parameter proposals taken from progressively closer approximations to the posterior, parameterised by decreasing ABC thresholds [33]. Research in this area has considered how to choose the sequence of decreasing thresholds and distance metrics [8,27], and how best to evolve the parameters from one approximation to produce parameter proposals for the next [1,3,12].
Another proposed strategy for overcoming the simulation bottleneck is the multifidelity ABC (MF-ABC) approach [28]. This approach assumes that, in addition to the model under investigation (termed the high-fidelity model), there also exists a low-fidelity model, depending on the same parameters, that is significantly faster to simulate, usually at the cost of being less accurate [24]. The multifidelity approach to parameter inference uses simulations from both the low-fidelity model and high-fidelity model to approximate the likelihood. The highfidelity model is simulated as little as possible, to reduce the simulation bottleneck, but just enough to ensure that the resulting posterior estimate is suitably accurate. This technique is related to multilevel approaches to ABC [17,34,35], which use a decreasing sequence of ABC-SMC thresholds to produce coupled estimates at different levels. Other approaches that can be interpreted in the multifidelity framework include lazy ABC [26], delayed acceptance ABC [10] and early rejection [25]. In each of these, low-fidelity simulations are sometimes used to reject parameters before completing a high-fidelity simulation. Importantly, the more general MF-ABC framework in [28] allows for early acceptance as well as early rejection.
A key observation that can be made about these two techniques for improving ABC performance is that they are orthogonal, in the sense that they improve different aspects of the ABC approach. ABC-SMC considers only improving the method for proposing parameters to use in simulations and does not directly affect the binary estimate of the likelihood. In contrast, MF-ABC makes no change to the parameter proposals, but instead directly alters the method used to estimate the likelihood by using a combination of both low-fidelity and highfidelity model simulations. The complementarity of these two approaches has previously been shown in the specific context of combining delayed acceptance with SMC [10]. Thus, combining the general multifidelity framework of [28] with SMC should therefore yield significant speedup over existing methods.
1.1. Outline. In this paper we bring together these two orthogonal approaches to speeding up ABC algorithms. We will introduce a combined multifidelity sequential Monte Carlo ABC algorithm (MF-ABC-SMC). Section 2 formulates the existing ABC algorithms briefly described above, and the techniques we can use to quantify their performance. We then show how these ABC approaches can be combined in Section 3, by incorporating the multifidelity technique with the sequential importance sampling approach to ABC-SMC. In Section 4 we then fully exploit the SMC framework to optimise the multifidelity approach, producing the MF-ABC-SMC algorithm in Algorithm 4.1. This new algorithm is applied in Section 5 to a heterogeneous network of Kuramoto oscillators in a hierarchical Bayes parameter estimation task, to produce low-variance ABC posterior estimates significantly faster than the classical ABC-SMC approach. Finally, in Section 6, we discuss some important open questions for further optimising the MF-ABC-SMC algorithm.
2. Theoretical background. Assume that the model we are seeking to calibrate is a map (usually stochastic) from parameters θ, taking values in a parameter space Θ, to an output y, taking values in data space Y. We denote this map as a conditional density f (· | θ) on Y, and term the drawing of y ∼ f (· | θ) as simulating the model, with y termed a simulation.
For Bayesian inference, we furthermore assume the existence of a prior distribution π(·) on Θ, and of the observed data y obs ∈ Y that will be used to calibrate the model. The model induces the likelihood of the observed data, written L(θ) = f (y obs | θ), which is a function of θ. As described previously, the goal of Bayesian inference is to infer the posterior distribution p(θ | y obs ) on Θ, given y obs and π(·). Bayes' Theorem equates the posterior to the product of likelihood and prior, p(θ | y obs ) = 1 ζ L(θ)π(θ), where the normalisation constant, ζ, ensures that p(θ | y obs ) is a probability density that integrates to unity.

Approximate Bayesian computation. Often, the model under consideration is suf-
ficiently complicated that the likelihood cannot be calculated, necessitating a likelihood-free approach. We assume that, while the value of any f (y | θ) is not available, we are still able to simulate y ∼ f (· | θ). Let d(y, y obs ) denote a metric that quantifies how close any simulation, y, is to the observed data, y obs . For a positive threshold value > 0, we can then define a neighbourhood Ω (d, y obs ) = {y ∈ Y | d(y, y obs ) < } of model simulations that are 'close' to the data. Typically the dataset, y obs , is constant for the parameter estimation task, and the distance metric, d, is pre-determined. Hence we will often drop the (d, y obs ) dependence from our notation and simply write Ω for the -neighbourhood of y obs under the distance function d.
For a given positive distance threshold > 0, ABC replaces the exact likelihood, L(θ), with the ABC approximation to the likelihood, which is, to leading order for small , approximately proportional to L(θ). The ABC approximation to the likelihood then induces the ABC posterior, where (similarly to ζ above) the constant Z ensures that the ABC posterior is a probability distribution with unit integral.
The importance sampling ABC algorithm [23,30] presented in Algorithm 2.1 (ABC-IS) presents a simple method for drawing samples from the ABC posterior, p (θ | y obs ). In addition to the data, y obs , and the prior, π(θ), we assume an importance probability distribution q(θ) = q(θ)/Z q defined by the function q(θ), where q(θ) > 0 for all θ in the support of π.
Note that we do not assume that the normalisation constant Z q is known. For each parameter proposal θ n ∼q(·) from the importance distribution, the model y n ∼ f (· | θ n ) is simulated and a weight w n = w(θ n , y n ) is assigned, using the weighting function, to produce a weighted sample {w n , θ n }.
Note that the general stopping criterion used in Algorithm 2.1 allows for the algorithm to terminate after (for example) a fixed number of non-zero weights, a fixed number of parameter proposals, N , when a fixed budget of total computational time or memory is reached, or any other more complicated combination of such conditions. Furthermore we note that ABC-IS is easily parallelised, although care must be taken to ensure that the chosen stopping condition is correctly applied in this case [16]. If we choose the importance distributionq = π equal to the prior, then Algorithm 2.1 (ABC-IS) is known as rejection sampling, which we refer to as ABC-RS.
For any arbitrary function F (·) defined on the parameter space Θ, we can estimate the expected value of F (θ) under p (θ | y obs ), such that Although this estimate is biased (except in the ABC-RS case), it is consistent, such that the mean squared error (MSE) ofF is dominated by the variance and decays to 0 on the order 1/N .

Sequential
n }, is produced by Algorithm 2.1 using importance distributionq T and threshold T = . Note that Algorithm 2.2 requires the specification of a threshold sequence t , the stopping conditions S t , and the perturbation kernels K t (· | θ) for t = 1, . . . , T [30]. We will not re-examine these aspects of ABC-SMC in detail in this paper, and will implement ABC-SMC using established techniques to select t , S t and K t [3,8,12,33].
n } Nt n=1 using Algorithm 2.1 (ABC-IS) with importance distributionq t , neighbourhood Ω t , and stopping condition S t .

3:
Define the next importance distribution,q t+1 , proportional to using Algorithm 2.1 (ABC-IS) with importance distributionq T , neighbourhood Ω T = Ω , and stopping condition S T .
One weakness of the population Monte Carlo approach taken in Algorithm 2.2 is the O(N 2 ) cost for each run of Algorithm 2.1 (ABC-IS), where N t ∼ N is the scale of each generation's sample size. To overcome this problem, the SMC sampler was adapted in [8] to reduce this cost to O(N ). When the calculation of each w then a significant computational burden can be alleviated through using the O(N ) sampling approach. However, in many practical problems, the calculation of each weight, w n ). We will focus on the latter setting and aim to reduce the cost of ABC-SMC through reducing the cost of calculating each w (t) n using the multifidelity approach described in Subsection 2.3.

2.
3. Multifidelity ABC. Subsection 2.2 describes how the SMC strategy provides parameter proposals such that simulation time is not wasted in regions of low likelihood. An orthogonal approach to improving ABC efficiency is to avoid computationally expensive simulations, where possible, by relying on the multifidelity framework [28]. We term the model of interest, f (· | θ), which maps parameter space Θ to an output space Y, as the high-fidelity model. We now assume that, in addition, there is a low-fidelity (i.e. approximate) model, f (· | θ), of the same physical system with the same parameter space Θ. The simulations of this model are denotedỹ ∼f (· | θ), taking values in the output spaceỸ, which may differ from Y.
Importantly, we assume that the low-fidelity model is computationally cheaper, in the sense that simulationsỹ ∼f (· | θ) of the low-fidelity model incur less computational burden than simulations y ∼ f (· | θ) of the high-fidelity model.
We can also assume that the experimental observations y obs ∈ Y can be mapped to the new data space, givingỹ obs ∈Ỹ. Similarly, we define an associated region, Ω =Ω (ỹ obs ,d) = {ỹ ∈Ỹ |d(ỹ,ỹ obs ) <˜ }, that is the˜ -neighbourhood of the observed data,ỹ obs , as a subset of the output spaceỸ, defined under the distance metricd. However, in the interests of clarity, we will assume for the remainder of this article that the output spaces of each model fidelity are such thatỸ = Y.
Similarly, we assume that the observed data are such thatỹ obs = y obs , the distance metrics are such thatd = d, and the ABC thresholds are such that˜ = , so thatΩ˜ = Ω .
Then, if the low-fidelity model is a good approximation to the high-fidelity model, the outputs y andỹ will be near, in some sense, and the distances from data d(y, y obs ) and d(ỹ, y obs ) will be correlated. However, to improve this correlation we will also allow for coupling between the two models, writingf (y,ỹ | θ) as the coupled density for (y,ỹ) with marginals that coincide with the independent models f (y | θ) andf (ỹ | θ). The benefit of this approach is that, with a judicious choice of coupling, we can produce a high-fidelity simulation y ∼ f (· |ỹ, θ) conditionally on a previously-simulated low-fidelity simulation,ỹ ∼f (· | θ), where the distances are more closely correlated. Furthermore, simulations of y ∼ f (· |ỹ, θ) may be less computationally burdensome than independent simulations of y ∼ f (· | θ).
For example, suppose the high-fidelity model is a Markovian stochastic dynamical system on the time horizon t ∈ [0, T ], and that the low-fidelity model is the same system on t ∈ [0, τ ] for τ < T [26]. The low-fidelity and high-fidelity models can clearly be simulated independently on [0, τ ] and [0, T ], respectively, and the low-fidelity simulation will, on average, be less computationally expensive than the high-fidelity simulation. However, the natural coupling of the models fixes y(t) =ỹ(t) over t ∈ [0, τ ], allowing the high-fidelity model to be simulated conditional on a low-fidelity simulation. Many other possible couplings exist for different multifidelity models, often involving shared random noise processes, and methods for coupling are currently an area of active research [6,18,28].
In order to apply the multifidelity framework to ABC parameter inference, recall that each weight, w n , generated by Algorithm 2.1 requires a simulation y n ∼ f (· | θ n ) from the highfidelity model. The multifidelity approach in [28] calculates the weight w n = w(θ n ,ỹ n , u n , y n ) by replacing the weighting function in Equation (2.2) with where (ỹ, y) ∼f (·, · | θ) are coupled multifidelity simulations, where u is a unit uniform random variable, and where α(θ,ỹ) ∈ (0, 1] is a positive continuation probability. Note that, in general, we could allow the continuation probability, α(θ,ỹ, u, y), to depend on all of the stochastic variables, but we will assume α(θ,ỹ) ∈ (0, 1] to be independent of u and y. The important consequence of the multifidelity weight in Equation (2.5) is that, by the specific order in which we simulate the variables, the weight w(θ,ỹ, u, y) may be calculated without the computational cost of simulating y. Given θ, we first simulateỹ ∼f (· | θ) from the low-fidelity model. This defines the continuation probability α(θ,ỹ) ∈ (0, 1]. Second, we generate the unit uniform random variable, u. If u ≥ α(θ,ỹ) ∈ (0, 1] then we can return w(θ,ỹ, u, y) without simulating y ∼ f (· | θ,ỹ) from the coupling, thus incurring the lower computational expense of only simulating fromf . Algorithm 2.3 is an adaptation of Algorithm 2.1 that returns a weighted sample {θ n , w n } N n=1 from the ABC posterior p (y obs | θ), with weights calculated using w in Equation (2.5). As with Algorithm 2.1, in the rejection sampling case whereq(θ) = π(θ), then we refer to this algorithm as MF-ABC-RS.

9:
Update w n ← w n + [I(y n ∈ Ω ) − w n ] α(θ n ,ỹ n ). Update w n ← π(θ n ) q(θ n ) w n . 12: until S = true. Proof. We note that the density of each z = (θ,ỹ, u, y) sampled by Algorithm 2.3 is which is independent ofỹ. Therefore, for any integrable F : Θ → R, we have the identity , each sampled pair (θ n , w n ) also incurs a computational cost, which we will denote T n , producing a total computational cost of T total = n T n for the sample. We define the observed efficiency of such a sample as follows. Given a computational cost of T n for each sampled pair (w n , θ n ), the observed efficiency of a weighted Monte Carlo sample is , which is expressed in units of effective samples per time unit.
Note that the theoretical efficiency is not just a characteristic of the ABC algorithm and of the models, but also of the numerical implementation and hardware of the computers generating the simulation. For example, the theoretical efficiency of Algorithm 2.
where Z = p(y ∈ Ω ) is the normalisation constant in Equation (2.1b), and where E ν denotes expectations with respect to probability density ν. Settingq = π, the theoretical efficiency of ABC-RS is .
Hence, the theoretical efficiency of ABC-IS is improved over ABC-RS by choosing an importance distribution,q, that is more likely than π to propose θ incurring smaller simulation times and high posterior likelihoods. In the case of Algorithm 2.3 withq(θ) = π(θ) (MF-ABC-RS), this performance measure has been used to determine a good choice of continuation probability, α(θ,ỹ) [28].
By using the theoretical efficiency, ψ, as a performance metric in the remainder of this paper, we will quantify the improvement in performance over Algorithms 2.1 to 2.3 that can be achieved by combining multifidelity and SMC techniques.
3. Multifidelity ABC-SMC. There are two distinct approaches to improving the performance of ABC parameter inference specified in Section 2. ABC-SMC proposes parameters θ n from a sequence of importance distributions that progressively approximate the target ABC approximation to the posterior. In contrast, MF-ABC enables sampled parameters to be weighted without necessarily having to produce a simulation, y n , from the high-fidelity model. In this section we present the main contribution of this paper, which is to combine these orthogonal approaches into an MF-ABC-SMC algorithm.
We will replicate the procedure of extending Algorithm 2.1 (ABC-IS) into Algorithm 2.2 (ABC-SMC) in the multifidelity context, focusing first on the O(N 2 ) SMC sampler based on sequential importance sampling. Suppose that, as in Subsection 2.2, we have a decreasing sequence of thresholds, 1 > · · · > T = , inducing the neighbourhoods, Ω 1 ⊇ · · · ⊇ Ω T , and a sequence of ABC posteriors, p t (θ | y obs ). In principle, we can replace the call of instead, and the SMC algorithm would proceed in much the same way. However, the key difficulty with implementing multifidelity ABC-SMC lies in the definition of the importance distribution,q t+1 , for generation t + 1, given the Monte Carlo sample, {θ Since the weights, w n , are now calculated using the multifidelity weight in Equation there is a positive probability that there exist w n < 0. Negative values of w n are generated wheneverỹ n ∈ Ω , but where we also simulate the high-fidelity model such that y n / ∈ Ω .
This leads to two problems with the existing definition of the importance distribution in Equation (2.4). Primarily, there may exist θ in the prior support with q t+1 (θ) < 0, unless the perturbation kernels, K t , are carefully designed to avoid this. Secondly, even if the K t could be chosen to guarantee q t+1 (θ) > 0 on the prior support, it is not clear that we can easily sample from the importance distributionq t+1 ∝ q t+1 when some θ (t) n have negative weights.
n } Nt n=1 from Algorithm 2.3 (MF-ABC-IS), using the neighbourhood Ω t , continuation probability α t , importance distributionr t , and stopping condition S t .
which is defined by the weightsŵ n } can be shown to be drawn from the alternative target distribution, where the difference between the new target distribution, ρ t , and the ABC posterior, p t , is given by for t = 1, . . . , T . The function δ t implies that the new target distribution (compared to π t ) contains additional density in regions of parameter space where it is more likely thatỹ ∈ Ω t but y / ∈ Ω t ; in other words, where w (t) n < 0 is more likely.
The importance distribution in Equation (3.1) effectively makes the SMC algorithm target ρ t at each generation instead of p t . However, the weighted samples {w n from each generation's weighting function, remain from the ABC posteriors p t . Hence, at any generation (and in particular at t = T ), we can produce an estimate of a p t -integrable function F : Θ → R, such that , is a consistent Monte Carlo estimate of F under the ABC posterior. ). An alternative sampling method that is linear in N is proposed in [8], which replaces the sequential importance sampling approach described above. However, as noted in [8], both the O(N ) and O(N 2 ) sampling algorithms require O(N ) simulations of (ỹ, y) ∼f (·, · | θ) orỹ ∼f (· | θ). In this multifidelity setting, we are assuming that simulation time comprises the vast majority of the computational burden of each calculation of Equation (3.3). The benefit of the multifidelity weight is that it reduces the computational burden of generating (ỹ, y) by sometimes requiringỹ alone. In Appendix A we show that the O(N ) SMC sampling algorithm dilutes this benefit. Therefore, in this work we will focus on the sequential importance sampling SMC sampler described in Algorithm 3.1.

4.
Adaptive MF-ABC-SMC. In Algorithm 3.1 (MF-ABC-SMC-α), in addition to an assumed sequence of ABC thresholds, t , perturbation kernels, K t , and stopping conditions S t , for t = 1, . . . , T , we also assume a given sequence of continuation probabilities, α t . In this algorithm, each importance distribution,r t+1 , is determined by the output at generation t. Methodologies for adaptively choosing the perturbation kernels, K t+1 , and the ABC thresholds, t+1 , based on the preceding generations' samples, have been explored in previous work [8,12]. In this section, we will consider the adaptive approach to choosing each generation's continuation probability, α t+1 , based on the simulation output of generation t. For simplicity, we will constrain the search for optimal α(θ,ỹ) to the piecewise constant for the constants η 1 , η 2 ∈ (0, 1]. Here, η 1 is the probability of generating y after a 'positive' lowfidelity simulation (whereỹ ∈ Ω ) and η 2 is the probability of generating y after a 'negative' low-fidelity simulation (whereỹ / ∈ Ω ). The goal of this section is to specify the values of the two parameters, η 1 and η 2 , that will give the largest theoretical efficiency, ψ. In previous work, we have derived the optimal values of η 1 and η 2 to use in the special case of Algorithm 2.3 (MF-ABC-RS) corresponding to rejection sampling, where the importance distributionq = π equal to the prior distribution [28]. We can now extend this analysis by finding optimal values of η 1 and η 2 to use in the more general case of Algorithm 2.3 (MF-ABC-IS). The key to this optimisation is the following lemma, which describes how the efficiency of Algorithm 2.3 varies with the continuation probabilities used. The lemma assumes a given importance distribution, q(θ), defined as the normalisation of the known non-negative function q(θ).
where the denominator is expressed as a function of (η 1 , η 2 ) such that The coefficients in ψ are given by the integrals then the minimum value of φ(η 1 , η 2 ) in Equation (4.2), and the optimal value of (η 1 , η 2 ) in the entire positive quadrant, are given bȳ Define the two functions, for x > 0, Then the minimum value of φ on the boundary, ∂H, of H is attained at the minimum of   1). The coefficient W can be written W = ZE π (p /q), as a scaling of the prior expectation of the ratio p /q. Larger values of W require q to be more concentrated (relative to π) in regions of high posterior density, which is a well-known characteristic of importance sampling [23]. Thus, neither Z nor W relate specifically to the multifidelity approach and are properties of ABC importance sampling.
The coefficientT lo represents the average time taken to simulateỹ from the low-fidelity model. The other two time-based coefficients,T hi,p andT hi,n , represent the average time taken to complete the simulation of (ỹ, y) ∼f from the coupling, conditional on whether or notỹ ∈ Ω . These coefficients therefore determine how much time can be saved by avoiding expensive simulations and stopping after generatingỹ ∼f (· | θ).
The key determinants of the success of the multifidelity technique are W fp and W fn and their tradeoff between the high-fidelity simulation costs,T hi,p andT hi,n . Equation (4.2) implies that the marginal cost to the efficiency of decreasing η 1 and η 2 is smaller for smaller values of W fp and W fn . These coefficients can be written as the two prior expectations, Thus, they are scalings of the probabilities of a false positive (whereỹ ∈ Ω but y / ∈ Ω ) and a false negative (whereỹ / ∈ Ω but y ∈ Ω ), respectively. Hence, small values of W fp and W fn correspond to at least one of the following cases. First, if the low-fidelity and high-fidelity  and t n = T (y n ); the distances from data,d n = d(ỹ n , y obs ) and d n = d(y n , y obs ); the importance densities q n = q(θ n ); and the continuation probabilities α n = α(θ n ,ỹ n ).
We now consider finding the optimal continuation probability, α , to be used in a new run of Algorithm 2.3, with the new importance distribution, q (θ), and the new ABC threshold, .

5:
Update the estimates in Equation (4.6) with the values stored at step 2, using importance distribution r t+1 (θ) and ABC threshold t+1 .

7:
Set α t+1 (ỹ, θ) = η 1 I(ỹ ∈ Ω ) + η 2 I(ỹ / ∈ Ω ). Thus, there is extra cost at each generation on the order O(N 2 ). However, in common with Note 3.1, we will assume that this cost is justified by our aim to reduce the large simulation burden that dominates the algorithm's run time.

Example: Kuramoto Oscillator Network.
To demonstrate the multifidelity and SMC approaches to parameter inference, we will infer the parameters of a Kuramoto oscillator model on a complete network, with stochastic heterogeneity in each node's intrinsic frequency. In Subsection 5.1 we consider the performance of the previously developed ABC algorithms introduced in Section 2 (ABC-RS, ABC-SMC and MF-ABC-RS) and demonstrate the orthogonal ways in which the SMC and multifidelity techniques improve performance. In Subsection 5.2 we apply the adaptive algorithm, Algorithm 4.1 (MF-ABC-SMC) to demonstrate that the efficiency of parameter estimation is significantly improved by combining the multifidelity and SMC approaches. The algorithms have been implemented in Julia [4] and the source code can be found at github.com/tpprescott/mf-abc-smc. The Kuramoto oscillator model is defined on a complete network of M nodes, where each node, i, has a dynamically evolving phase value, φ i , determined by the ordinary differential equation for i = 1, . . . , M . Each constant ω i , the intrinsic angular velocity, is an independent draw from a Cauchy distribution with median ω 0 and dispersion parameter γ. In addition to these two parameters, we have an interconnection strength K. Simulations of the ODE system are run over a fixed time interval t ∈ [0, T ], and we will assume fixed initial conditions φ i (0) = 0 for all i.
The multifidelity approach makes use of a low-dimensional approximation of the coupled oscillator dynamics, as described in [13,21,22]. The approximation is based on tracking the Daido order parameters, which are a set of complex-valued representations of the high- for positive integers n and the imaginary unit i. A system of coupled ODEs can be generated for the set of Z n . Under the assumption that Z n (t) = Z 1 (t) n , known as the Ott-Antonsen ansatz [13,21], the system can be reduced to a single ODE for Z 1 , which is known as the Kuramoto parameter. This complex-valued trajectory is usually represented by two real trajectories, corresponding to its magnitude R(t) = Z 1 (t)) and phase Φ(t) = arg(Z 1 (t)).
The approximation of the M -dimensional ODE system in Equation In colour are five trajectories for R(t) and Φ(t) for the high-fidelity model in Equation (5.1). Stochasticity arises from sampling ωi ∼ Cauchy(ω0, γ) for i = 1, . . . , 256. In black is the deterministic trajectory from the low-fidelity model in Equation (5.2). All simulations were completed using parameters K = 2, ω0 = π/3 and γ = 0.1. Note that the phase plot "unwraps" the trajectories of Φ(t) = arg(Z1(t)) to avoid 2πdiscontinuities.
is thus given by the two-dimensional ODE systeṁ with initial conditions (R(0),Φ(0)) = (1, 0), which directly simulates the low-dimensional representation of the M -dimensional state vector. The goal of this example is to infer the parameters (K, ω 0 , γ) based on synthetic data, generated by simulating a system of M = 256 oscillators with random angular velocities ω i over t ∈ (0, 30]. We record the trajectories R obs (t) and Φ obs (t) of the magnitude and phase of the Kuramoto parameter. The parameter values used to generate these data are (K = 2, ω 0 = π/3, γ = 0.1). The likelihood of the observed data under the model in Equation (5.1) with stochastic parameters ω i ∼ Cauchy(ω 0 , γ) is unavailable, and we must therefore resort to ABC inference, requiring repeated simulation.
Example trajectories of the high-fidelity and low-fidelity models in Equations (5.1) and (5.2) are given in Figure 5.1. The trajectories R obs (t), Φ obs (t), R(t), Φ(t),R(t) andΦ(t) on t ∈ [0, 30] are infinite dimensional. In order to easily compare trajectories, we will select a finite number of informative summary statistics from the trajectories, guided by the approximated system in Equation (5.2). We take Note that the low-fidelity model is deterministic. Therefore there is no meaningful definition of a coupling between the models at each fidelity: any simulation y ∼ f (· |ỹ, (K, ω 0 , γ)) = f (· | (K, ω 0 , γ)) from the high-fidelity model will be independent ofỹ.

Existing ABC algorithms.
We set independent uniform priors on [1,3] Table 5.1 Comparing existing ABC algorithms for sampling p0.5((K, ω0, γ) | y obs ) based on a uniform prior, using = 0.5. The stopping condition for Algorithms 2.1 and 2.3 is N = 6000. The threshold schedule for Algorithm 2. for the performance of each of these algorithms, quantified in terms of ESS, total simulation time, and observed efficiency (i.e. the ratio of the first two). Note that the ESS of the sample from Algorithm 2.2 (ABC-SMC) depends on only N 4 = 1500 weights, w (4) n , corresponding to the final generation. However, we will measure the observed efficiency by using the total time to simulate, which includes the total simulation time of the preceding generations. . This condition reflects a specification that we need each generation's sample to be, in some sense, 'good enough' to produce a reliable importance distribution that can be used in the next generation. Finally, we specified the perturbation kernels K t (· | (K, ω 0 , γ) (t) n ) at each generation to be Gaussians centred on the parameter value (K, ω 0 , γ) (t) n . The covariance matrices are diagonal matrices diag(σ and similarly for σ  γ . These perturbation kernels implement a typical choice for the covariance of using twice the empirical variance of the observed parameter values [3,12]. Note that we use this definition for the multifidelity case also, where weights w (t) n may be negative, since we are using the ρ t (as defined in Equation (3.2)) and not the p t as the target distributions. Further to these common inputs, the parameters ρ 1 and ρ 2 are the only additional algorithm parameters we need to specify to implement Algorithm 4.1 (MF-ABC-SMC). We set lower bounds of ρ 1 = ρ 2 = 0.01 on the allowed continuation probabilities, with the aim to limit the variability of w     that MF-ABC-SMC is reliably significantly less computationally expensive than ABC-SMC to produce comparable samples, with the average simulation times reflecting the 60% saving identified in Table 5.2.
Recall our initial observation that the average simulation time of the low-fidelity model, 10 µs, is approximately 2% of that of the high-fidelity model, 520 µs, with expectations taken over the prior. Given this initial difference in each model's average simulation time, the observed 148% increase in efficiency from Algorithm 2.2 to Algorithm 4.1 is determined by a number of other factors specific to the characteristics of SMC sampling, which we now explore.

5.2.2.
MF-ABC-SMC is more effective in early generations. In producing the samples summarised in Table 5  Finally, note that the sequence of t was chosen to be equal across ABC-SMC and MF-ABC-SMC to allow for a direct comparison between the efficiencies at each generation. This constraint produced an overall 60% saving in simulation time, and a 148% increase in efficiency. However, in practice, we should also aim to adaptively choose the sequence t (for either algorithm) to optimise performance. In the following subsection, we describe how incorporating the adaptive selection of t based on the preceding generations' output allows for a better comparison between the two algorithms. (MF-ABC-SMC) requires a decreasing sequence of t values to be pre-specified. Often, appro-priate values of t cannot be determined before any simulations have been generated. If the sequence of t decreases too slowly then the algorithm will take a long time to reach the target posterior; too quickly, and acceptance rates will be too low. As a result, rather than specifying a sequence of thresholds a priori, previous work in the SMC context [8] has explored choosing t+1 by predicting its effect on the ESS of that generation. In the spirit of this approach, we adaptively choose the sequence t by predicting its effect on the efficiency of that generation, as defined in Lemma 4.1.  1; ,q), also varies with . We assume that at each generation we have a target efficiency, ψ t , that enables a given ESS to be generated with a known computational budget.
This procedure produces a sequence of ABC thresholds designed such that each generation's efficiency is maintained at a target level. It alternates between finding continuation probabilities that maximise the efficiency at the preceding threshold, t , and then finding a value of t+1 < t such that the predicted efficiency matches the target. If t+1 = t then the target efficiency, ψ t+1 , needs to be reviewed as it is not achievable. Note that for the case of ABC-SMC, which is equivalent to MF-ABC-SMC with fixed continuation probabilities η 1 = η 2 = 1, we skip the optimisation in step a and use η 1 = η 2 = 1 in step b.
The strength of the multifidelity approach in this context is that the additional degrees of freedom afforded by the continuation probabilities allows for the ABC thresholds to decrease more quickly, while maintaining a target efficiency at each generation. This benefit is depicted in Figure 5.5. Each curve in the left-hand plot is the predicted efficiency ψ(1, 1; ,q t+1 ) as a function of . Each curve in the right-hand plot is the predicted efficiency ψ(η 1 , η 2 ; ,r t+1 ), again as a function of , where we have found optimal continuation probabilities, (η 1 , η 2 ).
For each algorithm, we choose the target efficiency at each generation t > 1 to equal the efficiency observed in generation 1. In our example, each generation has a stopping condition of ESS ≥ 400. As a result, by choosing to specify a constant target efficiency equal to the observed efficiency of generation t = 1, we effectively impose a constant target simulation budget for each generation. Since (for the example in Figure 5.5) we have specified four generations for each run of the adaptive versions of Algorithms 2.2 and 4.1, we have thus specified a fixed, equal total simulation time for each algorithm. In this setting, the results in Figure 5.5 show that the bias incurred by using an ABC approximation to the posterior with threshold > 0 is vastly reduced by implementing the multifidelity approach to SMC. Using MF-ABC-SMC, the sample produced in generation 4 is from p 0.11 (θ | y obs ), while ABC-SMC can only produce a sample from p 0.51 (θ | y obs ) for a similar computational cost. Thus, by incorporating the multifidelity approach into a method for adaptively selecting ABC thresholds, we can afford to allow the sequence t of MF-ABC-SMC thresholds to decrease much more rapidly, at no cost to the efficiency of the algorithm.
6. Discussion and conclusions. In this work we have examined how to integrate two approaches to overcoming the computational bottleneck of repeated simulation within ABC: the SMC technique for producing parameter proposals, and the multifidelity technique for using low-fidelity models to reduce overall simulation time. By combining these approaches, we have produced the MF-ABC-SMC algorithm in Algorithm 4.1.
The results in Section 5 demonstrate that the efficiency of sampling from the ABC posterior (measured as the ratio of the ESS to simulation time) can be significantly improved by using Some of these questions are specific to the implementation of multifidelity approaches. However, others arise from a re-evaluation of SMC implementation strategies in this new context. Below, we consider these two classes of question in turn.

Multifidelity implementation.
The key to maximising the benefit of MF-ABC-SMC is the ability to set a continuation probability based on the simulations generated during preceding generations. The estimates in Equation (4.6) of the quantities in Equation (4.3) are natural Monte Carlo approximations to the required integrals. In the SMC context, when generating the continuation probability for generation t + 1, each of these estimates could actually be constructed using the parameter proposals, importance weights, simulations, distances, and continuation probabilities of any (or all) generations 1 ≤ s ≤ t, not just generation t. Future work should clarify how best to combine many generations' samples into estimates of Equation (4.3), and the potential for improvement that might arise from this.
Another question arises when we break the assumption made at the start of Subsection 2.3, and no longer assume that the output spaces of each model fidelity are such thatỸ = Y. In general, the observed data,ỹ obs = y obs , the distance metrics,d(ỹ,ỹ obs ) = d(y, y obs ), and the thresholds,˜ = , may all be distinct. In this case, any estimate,w, of I(y ∈ Ω ) can be used in place of I(ỹ ∈ Ω ) to give a multifidelity acceptance weight of the form w(θ,ỹ, u, y) =w(θ,ỹ) + I(u < α(θ,ỹ)) α(θ,ỹ) (I(y ∈ Ω ) −w(θ,ỹ)) .
Note that, in the case of equal output spaces, we might consider usingw(y n ) = I(ỹ ∈ Ω˜ ) for distinct thresholds,˜ = . However, in general,w(θ,ỹ) may also encompass completely different output spaces based on distinct modelling frameworks for the same system (albeit with the same parameter space). In a similar way to evolving acceptance probabilities across generations, we could also evolve the estimatew across generations, using the information gathered from repeated simulation of both high-fidelity and low-fidelity models to better approximate I(y ∈ Ω ) and thus reduce our reliance on the high-fidelity model.
The form of continuation probability, α(θ,ỹ), defined in Equation (4.1) implements a multifidelity algorithm that provides a single continuation probability forỹ ∈ Ω and another forỹ / ∈ Ω , independently of the parameter value. There may be significant improvements to the multifidelity approach available through making α(θ,ỹ) depend more generally on θ andỹ.
For example, it is probable that there should be less need to simulate y after generatingỹ such that d(ỹ, y obs ) , than if d(ỹ, y obs ) = 1.01 . However, with α as defined in Equation (4.1), these two cases are treated equally. There is likely to be significant potential for improved performance from exploring less simplistic forms for the continuation probability.
Finally, as has been noted in previous work on multifidelity ABC [28], there is much potential in being able to use multiple low-fidelity models, beyond a single low-fidelity approximation. If there exist multiple low-fidelity models, different generations of MF-ABC-SMC may allow us to progressively focus on using the most efficient model, and identify the specific regions of parameter space for which one model or another may bring most benefit for parameter estimation.

SMC implementation. Previous research into the implementation of ABC-SMC has
ensured that the importance distribution formed from the preceding generation, as given in Equation (2.4), is optimal, by choice of the perturbation kernels, K t [12]. This has typically been treated as a requirement to trade-off a wide exploration of parameter space against a high acceptance rate. We have replaced the acceptance rate by the theoretical efficiency as the quantification of an ABC algorithm's performance. Therefore, since we now explicitly include the simulation time in the definition of the algorithm's performance, the perturbation kernels that optimise the tradeoff between efficiency and exploration may be reformulated and hence different optima used.
In Subsection 5.2 we applied a widely-used strategy for determining the perturbation kernels [3]. This strategy has been justified only in the context of positive weights and an importance distribution that approximates the ABC posterior. However, the importance distribution used in MF-ABC-SMC potentially makes this choice of perturbation kernels suboptimal. It remains to extend existing results on the optimality of perturbation kernels, such as those in [12], to apply to importance distributions of the form in Equation (3.1) that approx-imate the alternative target distribution in Equation (3.2). There is therefore justification for reopening the question of specifying optimal perturbation kernels for SMC, in the context of both including simulation time in the performance tradeoff and for dealing with multifidelity samples with negative weights.
We have restricted the choice of SMC sampler to the O(N 2 ) PMC sampler, due to the effect of the O(N ) sampler in diluting the benefits of the multifidelity approach, as discussed in Section 3 and Appendix A. Although this choice can be justified when simulation times dominate the algorithm, future work in this area should seek to optimise a multifidelity approach in the context of the more efficient sampling technique. Given that the observed effect of multifidelity ABC is to significantly reduce the simulation time per parameter proposal, and thereby make the simulation cost at each iteration much less dominant, this extension will be necessary to ensure optimal performance of MF-ABC-SMC in larger sample sizes.
In Subsection 5.2.3 we described how to adapt Algorithms 2.2 and 4.1 to implement an adaptive sequence of t as an approach to minimising bias for a fixed computational budget.
The strategy we used was to choose t to maintain an efficiency as close as possible to a target, set to equal the observed efficiency of the first generation. Further work in this area should investigate the use of more sophisticated strategies for choosing each t . This question relates closely to the sequence of stopping criteria. In Subsection 5.2, we constrained the effective sample size at each generation to be at least 400 to ensure a relatively low variance in each generation's sample, but this choice was made arbitrarily. Future work should therefore consider how to best achieve the ultimate goal of the SMC algorithm: a sample from a final generation with minimal bias relative to the true posterior, with a small variance, and constructed quickly. This goal should be achieved by optimising the complex, interdependent choices of stopping criteria, ABC thresholds, continuation probabilities and perturbation kernels. n , a perturbation to z is proposed with density The proposal is accepted with a Metropolis-Hastings acceptance probability based on the ratioŵ , for the non-negative multifidelity weight w t (z) = I(ỹ ∈ Ω t ) + I(u < α t ) α t (I(y ∈ Ω t ) − I(ỹ ∈ Ω t )) .
The MCMC step thus produces a sample point for the new generation, z The benefit of the multifidelity weight is that, for u < α t , we only need to generateỹ rather than (ỹ, y), avoiding a significant computational cost. However, the weight update step given in Equation (A.1) requires two weight calculations dependent on z (t+1) n , each of which may require the expensive computation of (ỹ, y) ∼f (·, · | θ). Indeed, the cost saving of needing to generateỹ ∼f (· | θ) alone for the weight update step exists only if for S 1 (i.e. equal steady states, and thus equal values for γ/K) can be distinguished by the speed at which they reach their steady state, which we will infer through S 3 . Note that the sampling point, T 1/2 , used in S 3 is chosen to be relevant to the observed data specifically, and aims to distinguish any simulated trajectories from R obs and Ω obs in particular. Thus we select S 3 to identify the scale of γ and K. Hence, we assume that these three summary statistics will be sufficient to identify the parameters.   ∈ Ω t ) used to produce the multifidelity sample for each of the 50 replicates of generations 2, 4, 6 and 8.