{'fraction_non_alphanumeric': 0.04761829000540936, 'fraction_numerical': 0.012696089349922042, 'mean_word_length': 4.466602887458689, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient Multifidelity Likelihood-Free Bayesian Inference with Adaptive Computational Resource Allocation

Likelihood-free Bayesian inference algorithms are popular methods for calibrating the parameters of complex, stochastic models, required when the likelihood of the observed data is intractable. These algorithms characteristically rely heavily on repeated model simulations. However, whenever the computational cost of simulation is even moderately expensive, the significant burden incurred by likelihood-free algorithms leaves them unviable in many practical applications. The multifidelity approach has been introduced (originally in the context of approximate Bayesian computation) to reduce the simulation burden of likelihood-free inference without loss of accuracy, by using the information provided by simulating computationally cheap, approximate models in place of the model of interest. The first contribution of this work is to demonstrate that multifidelity techniques can be applied in the general likelihood-free Bayesian inference setting. Analytical results on the optimal allocation of computational resources to simulations at different levels of fidelity are derived, and subsequently implemented practically. We provide an adaptive multifidelity likelihood-free inference algorithm that learns the relationships between models at different fidelities and adapts resource allocation accordingly, and demonstrate that this algorithm produces posterior estimates with near-optimal efficiency.

understand the relationship between model parameters and behaviour. In time-sensitive and costaware applications, the typical computational burden of such simulation-based methods makes them impractical. Multifidelity methods, reviewed by Peherstorfer et al. [2016Peherstorfer et al. [ , 2018, are a family of approaches that exploit information gathered from simulations, not only of a single model of interest, but also of additional approximate or surrogate models. In this article, the term model refers to the underlying mathematical abstraction of a system in combination with the computer code used to implement simulations. Thus, 'model approximation' may refer to mathematical simplifications and/or approximations in numerical methods. The fundamental challenge when implementing multifidelity techniques is the allocation of computational resources between different models, for the purposes of balancing a characteristic trade-off between maintaining accuracy and saving computational burden.
In this work, we consider a specific outer-loop application that arises in Bayesian statistics, the goal of which is to calibrate a parametrised model against observed data. Bayesian inference uses the likelihood of the observed data to update a prior distribution on the model parameters into a posterior distribution, according to Bayes's rule. In the situation where the likelihood of the data cannot be calculated, we rely on so-called likelihood-free methods that provide estimates of the likelihood by comparing model simulations to data. For example, approximate Bayesian computation (ABC) is a widely-known likelihood-free inference technique [Sisson et al., 2020, Sunnåker et al., 2013, where the likelihood is typically estimated as a binary value, recording whether or not the distance between a simulation and the observed data falls within a given threshold. Other likelihood-free methods are also available, such as pseudo-marginal methods and Bayesian synthetic likelihoods (BSL). In this work, we develop a generalised likelihood-free framework for which ABC, pseudo-marginal and BSL can be expressed as specific cases, as described in Section 2.1.
The significant cost of likelihood-free inference has motivated several successful proposals for improving the efficiency of likelihood-free samplers, such as (in the specific context of ABC) ABC-MCMC [Marjoram et al., 2003] and ABC-SMC [Sisson et al., 2007, Toni et al., 2009, Del Moral et al., 2011. These approaches aim to efficiently explore parameter space by avoiding the proposal of low-likelihood parameters, reducing the required number of expensive simulations required and reducing the ABC rejection rate. However, an 'orthogonal' technique for improving the efficiency of likelihood-free inference is to instead ensure that each simulation-based likelihood estimate is, on average, less computationally expensive to generate.
In previous work, Baker [2020, 2021] investigated multifidelity approaches to likelihoodfree Bayesian inference [Cranmer et al., 2020], with a specific focus on ABC [Sisson et al., 2020, Sunnåker et al., 2013. Suppose that there exists a low-fidelity approximation to the parametrised model of interest, and that the approximation is relatively cheap to simulate. Monte Carlo estimates of the posterior distribution, with respect to the likelihood of the original high-fidelity model, can be constructed using the simulation outputs of the low-fidelity approximation. Prescott and Baker [2020] showed that using the low-fidelity approximation introduces no further bias, so long as, for any parameter proposal, there is a positive probability of simulating the high-fidelity model to check and potentially correct a low-fidelity likelihood estimate. The key to the success of the multifidelity ABC (MF-ABC) approach is to choose this positive probability to be suitably small, thereby simulating the original model as little as possible, while ensuring it is large enough that the variance of the resulting Monte Carlo estimate is suitably small. The result of the multifidelity approach is to reduce the expected cost of estimating the likelihood for each parameter proposal in any Monte Carlo sampling algorithm. In subsequent work, Prescott and Baker [2021] showed that this approach integrates with sequential Monte Carlo (SMC) sampling for efficient parameter space exploration [Toni et al., 2009, Del Moral et al., 2011, Drovandi and Pettitt, 2011. Thus, the synergistic effect of combining multifidelity and with SMC to improve the efficiency of ABC has been demonstrated.
Multifidelity ABC can be compared with previous techniques for exploiting model approximation in ABC, such as Preconditioning ABC [Warne et al., 2021a], Lazy ABC (LZ-ABC) [Prangle, 2016], and Delayed Acceptance ABC (DA-ABC) [Christen andFox, 2005, Everitt andRowińska, 2021]. Similarly to sampling techniques such as SMC, the preconditioning approach seeks to explore parameter space more efficiently, by proposing parameters for high-fidelity simulation with greater low-fidelity posterior mass. In contrast, each of MF-ABC, LZ-ABC, and DA-ABC seeks to make each parameter proposal quicker to evaluate, on average, by using the output of the low-fidelity simulation to directly decide whether to simulate the high-fidelity model. In both LZ-ABC and DA-ABC, a parameter proposal is either (a) rejected early, based on the simulated output of the low-fidelity model, or (b) sent to a high-fidelity simulation, to make a final decision on ABC acceptance or rejection. The distinctive aspect of MF-ABC is that step (a) is different; it is not necessary to reject early to avoid high-fidelity simulation. Instead the low-fidelity simulation can be used to make the accept/reject decision directly. In both DA-ABC and MF-ABC, the decision between (a) or (b) is based solely on whether the low-fidelity simulation would be accepted or rejected. In contrast, LZ-ABC allows for a much more generic decision of whether to simulate the high-fidelity model, requiring an extensive exploration of practical tuning methods.
In this paper, we will show that the multifidelity approach can be applied to any simulationbased likelihood-free inference methodology, including but not limited to ABC. We achieve this by developing a generalised framework for likelihood-free inference, and deriving a multifidelity method to operate in this framework. A successful multifidelity likelihood-free inference algorithm requires us to determine how many simulations of the high-fidelity model to perform, based on the parameter value and the simulated output of the low-fidelity model. We provide theoretical results and practical, automated tuning methods to allocate computational resources between two models, designed to optimise the performance of multifidelity likelihood-free importance sampling.

Outline
In Section 2 we review likelihood-free Bayesian inference through constructing a generalised likelihoodfree framework. Section 3 shows how the theory underpinning MF-ABC, as introduced by Prescott and Baker [2020], can be applied in this general likelihood-free Bayesian context. We analyse the performance of the resulting multifidelity likelihood-free importance sampling algorithm. The main results of this paper are set out in Section 3.3, in which we determine the optimal allocation of computational resources between the two models to achieve the best possible performance of multifidelity inference. Section 4 explores how to practically allocate computation between model fidelities, by adaptively evolving the allocation in response to learned relationships between simulations at each fidelity across parameter space. We illustrate adaptive multifidelity inference by applying the algorithm to a simple biochemical network motif in Section 5. We show that, using a low-fidelity Michaelis-Menten approximation together with the exact model (both simulated using the exact algorithm of Gillespie [1977]) our adaptive implementation of multifidelity likelihoodfree inference can achieve a quantifiable speed-up in constructing posterior estimates to a specified variance and with no additional bias. Code for this example, developed in Julia 1.6.2 [Bezanson et al., 2017], is available at github.com/tpprescott/mf-lf. Finally, in Section 6 we discuss how greater improvements may be achieved for more challenging inference tasks.

Likelihood-free inference
We consider a stochastic model of the data generating process, defined by a distribution with parametrised probability density function, f (· | θ), where the parameter vector θ takes values in a parameter space Θ. For any θ, the model induces a probability density, denoted f (y | θ), on observable outputs, with y taking values in an output space Y. We note that the model is usually implemented in computer code to allow simulation, through which outputs y ∈ Y can be generated. We write y ∼ f (· | θ) to denote simulation of the model f given parameter values θ. Taking the experimentally observed data y 0 ∈ Y, we define the likelihood function to be a function of θ using the density, L(θ) = f (y 0 | θ), of the observed data under this model. Bayesian inference updates prior knowledge of the parameter values, θ ∈ Θ, which we encode in a prior distribution with density π(θ). The information provided by the experimental data, encoded in the likelihood function, L(θ), is combined with the prior using Bayes' rule to form a posterior distribution, with density where Z = L(θ)π(θ) dθ normalises π(· | y 0 ) to be a probability distribution on Θ. For a given, arbitrary, integrable function G : Θ → R, we take the goal of the inference task as the production of a Monte Carlo estimate of the posterior expectation, conditioned on the observed data.

Approximating the likelihood with simulation
In most practical settings, models tend to be sufficiently complicated that calculating L(θ) = f (y 0 | θ) for θ ∈ Θ is intractable. In this case, we exploit the ability to produce independent simulations from the model, In the following, we will slightly abuse notation by using the shorthand y ∼ f (· | θ) to represent K independent, identically distributed draws from the parametrised distribution f (· | θ).
Given the observed data, y 0 , we can define a real-valued function, referred to as a likelihood-free weighting, which varies over the joint space of parameter values and simulation outputs. Here, ω is a function of a parameter value, θ, and a vector, y, of stochastic simulations. For a fixed θ, we can take conditional expectations of ω(θ, y) over the probability density of simulations, f (y | θ), to define an approximate likelihood function, where L ω (θ) is assumed to be approximately equal to the modelled likelihood function, L(θ), up to a constant of proportionality. Note that, given θ, the random value of the likelihood-free weighting, ω(θ, y), determined by K stochastic simulations, y ∼ f (· | θ), is a Monte Carlo estimate of the approximate likelihood function, L ω (θ).
The approximate likelihood function is used to define the likelihood-free approximation to the posterior, where the normalisation constant Z ω = L ω (θ)π(θ) dθ ensures that π ω is a probability distribution. The likelihood-free approximation to the posterior, π ω (θ | y 0 ), subsequently induces a biased estimate ofḠ, given byḠ In this situation, the success of likelihood-free inference depends on ensuring that the likelihoodfree weighting, ω(θ, y) is chosen such that the squared difference (Ḡ −Ḡ ω ) 2 between the posterior expectation,Ḡ, and its likelihood-free approximation,Ḡ ω , is as small as possible.

Example: ABC
Approximate Bayesian computation (ABC) is a widely-used example of likelihood-free inference, where is the random fraction of K simulations, y k , which are within a distance of of the observed data, as measured by the metric d. Sisson et al. [2007Sisson et al. [ , 2020 note that it is standard practice to choose K = 1. Taking the conditional expectation of ω ABC , given θ, this likelihood-free weighting induces the ABC approximation to the likelihood, for any K. Under appropriate choices of d and , the approximate likelihood function, L ABC (θ), may be considered approximately proportional to the likelihood, L(θ).

Likelihood-free importance sampling
A simple approach to estimating the likelihood-free approximate posterior mean,Ḡ ω , is to use importance sampling. We assume that parameter proposals θ i ∼ q(·), for i = 1, . . . , N , can be sampled from a given importance distribution, the support of which must include the prior support.
In practice, we need only know importance density values, q(θ), up to a multiplicative constant, but for simplicity we assume that q(θ) is known. We also assume that we have access to the prior probability density, π(θ).
The likelihood-free importance sampling algorithm is described in Algorithm 1. This algorithm requires the specification of an importance distribution, q, and a likelihood-free weighting, ω(θ, y), with conditional expectation, L ω (θ) = E(ω | θ). The output of Algorithm 1,Ĝ, is an estimate of the likelihood-free approximate posterior expectation,Ḡ ω = E πω (G | y 0 ). In Theorem 1, we prove the standard result thatĜ is a consistent estimate ofḠ ω , and quantify the dominant behaviour of the mean squared error in the limit of large sample sizes, N → ∞. For notational simplicity, we define the function ∆(θ) = G(θ) −Ḡ ω to recentre G at the approximate posterior mean, and denote the Monte Carlo error betweenĜ andḠ ω as the estimated mean value of ∆, denoted bŷ ∆ =Ĝ −Ḡ ω .
Theorem 1. For the weighted sample values (θ i , w i ) produced in each iteration of Algorithm 1, let w denote the random value of the weight w i , and let ∆ denote the random value of ∆(θ i ). The mean squared error (MSE) of the output,Ĝ, of Algorithm 1, as an estimator for the approximate posterior expectation,Ḡ ω , is given to leading order by Thus,Ĝ is a consistent estimator ofḠ ω .
We write the function Φ(r, s) = (r/s −Ḡ ω ) 2 , and note that the MSE is the expected value of the function∆ 2 = Φ(R, S). Using the delta method, we take expectations of the second-order Taylor expansion of Φ(R, S) about (µ R , µ S ) = (E(R), E(S)), to give Taking expectations with respect to N independent draws of (θ, y) with density f (y | θ)q(θ), it is straightforward to write where we recall that E(w) = Z ω = L ω (θ)π(θ) dθ is the normalising constant in Equation (4b).
We substitute these expectations into the Taylor expansion of E(∆ 2 ), noting that the leading-order term, Φ(µ R , µ S ), is zero. Thus, we can write the dominant behaviour of the MSE as as N → ∞. Substituting into this expression the definitions of R and S as summations of N independent identically distributed random variables, we have and Equation (6) follows, on noting that E(w∆) = 0 and that E(w) = Z ω .
Theorem 1 determines the leading-order behaviour of the MSE of the output of Algorithm 1 in terms of sample size. We can also quantify the performance of this algorithm in terms of how the MSE decreases with increasing the overall computational budget.
Corollary 2. Let the computational cost of each iteration of Algorithm 1 be denoted by the random variable C. The leading order behaviour of the MSE ofĜ as an estimate ofḠ ω is as the total simulation budget C tot → ∞.
Proof. As the given computational budget increases, C tot → ∞, the Monte Carlo sample size that can be produced in that budget increases on the order of N ∼ C tot /E(C). On substituting this expression into Equation (6), the result follows.
We can use the leading-order coefficient of 1/C tot in Equation (7) to quantify the performance of likelihood-free importance sampling. Importantly, this expression explicitly depends on the expected computational cost, C, of each iteration of Algorithm 1. In the importance sampling context, the optimal importance distribution q should seek to minimise this coefficient, by trading off a preference for parameter values with lower computational burden against ensuring small variability in the weighted errors, w∆. However, for simplicity, we will assume in this paper that q is fixed. Instead, we seek ways to use model approximations to directly reduce the leading-order coefficient in Equation (7), based on the identified trade-off between decreasing computational burden, C, and controlling the variance of the weighted error, w∆.

Multifidelity inference
In Corollary 2, the performance of Algorithm 1 is quantified explicitly in terms of how the Monte Carlo error between the estimate,Ĝ, and the approximated posterior mean,Ḡ ω , decays with increasing computational budget, C tot . It initially appears reasonable to conclude that the linear dependence of the performance on the expected iteration time, E(C), implies that if we can speed up the simulation step of Algorithm 1, then we can significantly reduce the MSE for a given computational budget.
Suppose that there exists an alternative model that we can use in Algorithm 1 in place of the original model, f (· | θ), such that the expected computation time for each iteration, E(C), is significantly reduced. There are two important issues that prevent this being a viable option for improving the efficiency of likelihood-free inference. The first problem is that we need to be able to quantify the effect of the alternative model on the ratio E(w 2 ∆ 2 )/E(w) 2 to ensure that the overall performance of the algorithm is improved. It is not sufficient to show that the computational burden of each iteration is reduced, if too many more iterations are subsequently required to achieve a specified MSE.
The second problem arises from the observation that the limiting value ofĜ, as output from Algorithm 1, isḠ ω , with residual bias, recalling thatḠ ω is the approximate posterior expectation induced by L ω (θ) = E(ω | θ), and where the approximand,Ḡ, is the posterior expectation induced by the likelihood, L(θ) = f (y 0 | θ). We will identify this limiting residual squared bias, (Ḡ ω −Ḡ) 2 , as the fidelity of the model/likelihoodfree weighting pair. We emphasise here that the fidelity depends both on the model and the likelihood-free weighting used in Algorithm 1, and is contextual to the target function, G. For a given posterior mean,Ḡ, a model and likelihood-free weighting pair for which the value of (Ḡ ω −Ḡ) 2 is small is termed high-fidelity, while larger values of (Ḡ ω −Ḡ) 2 are termed low-fidelity. Thus, if we use an alternative model in place of f in Algorithm 1, the model (and likelihood-free weighting) may be too low-fidelity, in the sense of having too large a residual squared bias versus the posterior expectation of interest,Ḡ.
The multifidelity framework overcomes both these problems, by removing the need for a binary choice between the expensive model of interest and its cheaper alternative. Instead, we carry out likelihood-free inference using information from both models. We introduce the multifidelity likelihood-free importance sampling algorithm in Section 3.1. In Section 3.2, we show that multifidelity likelihood-free importance sampling loses no fidelity versus high-fidelity importance sampling. Section 3.3 contains the main analytical results of this paper, in which we explore the conditions under which multifidelity inference can improve the performance of likelihood-free importance sampling with finite computational budgets, as quantified by the leading-order characterisation of the MSE given in Equation (7).

Multifidelity likelihood-free importance sampling
We denote the high-fidelity model and likelihood-free weighting as f hi and ω hi , respectively. The likelihood under the high-fidelity model is denoted L hi (θ) = f hi (y 0 | θ), and is assumed to be intractable. Following the notation introduced in Equation (4), the high-fidelity pair f hi and ω hi induce the approximate likelihood, L hi (θ) = E(ω hi | θ) and the corresponding likelihood-free approximation to the posterior expectation,Ḡ hi . We further assume that simulating each y hi ∼ f hi (· | θ) is computationally expensive. This computational expense motivates the use of an approximate, low-fidelity model and likelihood-free weighting, denoted f lo and ω lo , respectively, inducing the approximate likelihood, L lo (θ) = E(ω lo | θ), and corresponding likelihood-free approximation to the posterior expectation,Ḡ lo . We note that the low-fidelity model, f lo , induces its own likelihood, L lo = f lo (y 0 | θ), but assume that this remains intractable, requiring instead the simulation-based Bayesian approach. However, we assume that simulations of the low-fidelity model, y lo ∼ f lo (· | θ), are significantly cheaper to produce compared to simulations of the high-fidelity model.
Given the models f lo and f hi , we will term the joint distribution f mf (y lo , y hi | θ) a multifidelity model when f mf has marginals equal to the low-and high-fidelity densities, f lo (y lo | θ) and f hi (y hi | θ). The models may be conditionally independent, such that f mf (y lo , y hi | θ) = f lo (y lo | θ)f hi (y hi | θ), in which case simulations at each model fidelity can be carried out independently given θ. Furthermore, if the simulations are conditionally independent, this means that the resulting likelihood-free weights, ω lo (θ, y lo ) and ω hi (θ, y hi ), are also conditionally independent.
However, in the more general definition of the multifidelity model as a joint distribution, we allow for coupling between the two fidelities. Conditioned on the low-fidelity simulations, y lo , and on parameter values, θ, we can produce a coupled simulation, y hi , from the density Coupling imposes correlations between the resulting likelihood-free weights, ω lo (θ, y lo ) and ω hi (θ, y hi ), which thus allows evaluated values of ω lo to provide more information about unknown values of ω hi , thereby acting as a variance reduction technique [Owen, 2013]. Given the (coupled) multifidelity model, we can calculate a multifidelity likelihood-free weighting as follows.
Definition 3. Let M be any non-negative integer-valued random variable, with conditional probability mass function p(· | θ, y lo ), and with a positive conditional mean, Given a parameter value, θ, we define z = (y lo , y hi,1 , y hi,2 , . . . , y hi,m ), (8a) noting that each y hi,i may be coupled to the low-fidelity simulation y lo ∼ f lo (· | θ). We combine Equations (8a) to (8d) to write the density of z as φ(z | θ). We further define the multifidelity likelihood-free weighting function, as the low-fidelity likelihood-free weighting, corrected by a randomly drawn number, M = m, of conditionally independent high-fidelity likelihood-free weightings. Taking expectations over z, we write as the multifidelity approximation to the likelihood.
Given M = m, only m replicates of y hi,i ∼ f hi (· | θ, y lo ) need to be simulated for ω mf (θ, z) to be evaluated. Thus, whenever m = 0, this means that no high-fidelity simulations need to be completed for ω mf (θ, z) to be calculated, removing the high-fidelity simulation cost from that iteration. Algorithm 2 presents the adaptation of the basic importance sampling method of Algorithm 1 to incorporate the multifidelity weighting function. The simulation step, y ∼ f (· | θ), in Algorithm 1 is replaced by the MF-Simulate function in Algorithm 2.

Accuracy of multifidelity inference
We observe that using f hi and ω hi in Algorithm 1 produces an estimate of the high-fidelity approximate posterior expectation,Ḡ hi . In Proposition 4, we show that the multifidelity approximate likelihood, L mf (θ) = E(ω mf (θ, z) | θ), is equal to the high-fidelity approximate likelihood, L hi (θ) = E(ω hi (θ, y hi ) | θ). As a result, Algorithm 2 also produces a consistent estimate of the high-fidelity approximate posterior expectation,Ḡ hi .
Proposition 4. The multifidelity approximation to the likelihood, L mf (θ) = E(ω mf | θ), is equal to the high-fidelity approximation to the likelihood, L hi (θ) = E(ω hi | θ). Therefore, the estimatê G mf produced by Algorithm 2 is a consistent estimate of the high-fidelity approximate posterior Algorithm 2 Multifidelity likelihood-free importance sampling.
For ω mf in Equation (8e), calculate the weight Proof. We take the expectation of ω mf conditional on (θ, y lo , m), to find Further taking expectations over the random integer M , which has conditional expected value µ(θ, y lo ), gives Further taking expectations with respect to y lo , it follows that L ω mf (θ) = L ω hi (θ). Therefore, the likelihood-free approximate posteriors, π ω mf = π ω hi , are equal and thusĜ mf is a consistent estimate ofḠ ω mf =Ḡ ω hi , as required.
In the limit of infinite computational budgets, the estimate produced by multifidelity importance sampling, in Algorithm 2, is as accurate as the estimate produced by high-fidelity importance sampling, in Algorithm 1 using f hi and ω hi . However, we still need to show that the performance of Algorithm 2 exceeds that of Algorithm 1 in the practical context of limited computational budgets. In Section 3.3, we introduce a method to quantify the performance of Algorithms 1 and 2 and show that the performance of multifidelity inference is strongly determined by the distribution of M , the random number of high-fidelity simulations required at each iteration.

Comparing performance
Corollary 2 determines the leading-order behaviour of the MSE of the output of Algorithm 1 as the computational budget increases. A similar result applies to the output of Algorithm 2. We compare two settings: first, using Algorithm 1 with the high-fidelity model, f hi , and likelihood-free weighting, ω hi . Each iteration has computational cost denoted C hi , and produces a weighted Monte Carlo sample with weights w i as independent draws of the random variable w hi . The output of Algorithm 1 is denotedĜ hi , with Monte Carlo error∆ hi =Ĝ hi −Ḡ hi . The MSE for Algorithm 1 has leading-order behaviour as the total simulation budget C tot → ∞, where ∆(θ) = G(θ) −Ḡ hi .
Second, we use Algorithm 2 with the multifidelity model f mf and likelihood-free weighting, ω mf . Each iteration has computational cost denoted C mf , and produces a weighted Monte Carlo sample with weights w i as independent draws of the random variable w mf . The output of Algorithm 2 isĜ mf , with Monte Carlo error∆ mf =Ĝ mf −Ḡ hi . The MSE for Algorithm 2 has leading-order behaviour as the total simulation budget C tot → ∞, where again ∆(θ) = G(θ) −Ḡ hi .
The main result of the paper is given in Theorem 5. For concreteness, we assume that the random variable, M , determining the required number of high-fidelity simulations in each iteration of Algorithm 2, is Poisson distributed, conditional on the parameter value and low-fidelity simulation output. We show that, for a given multifidelity model and likelihood-free weightings, the mean function, µ, for M determines the performance of Algorithm 2 relative to Algorithm 1.
Theorem 5. Assume that the random number of high-fidelity simulations, M , required in each iteration of Algorithm 2 is Poisson distributed with conditional mean µ(θ, y lo ). Let c hi (θ) [respectively, c lo (θ) and c hi (θ, y lo )] be the expected time taken to simulate y hi ∼ f hi (· | θ) [respectively, to simulate y lo ∼ f lo (· | θ) and to produce the coupled high-fidelity simulation y hi ∼ f hi (· | θ, y lo )].
Further, assume that the computational cost of each iteration of Algorithm 1 and Algorithm 2 can be approximated by the dominant cost of simulation alone, neglecting the costs of the other calculations.
The performance of Algorithm 2, quantified in Equation (11), exceeds the performance of Algorithm 1, quantified in Equation (10) for the constantsc for the functions and the joint density ρ(θ, y lo ) = f lo (y lo | θ)q(θ).
The performance metrics of Algorithm 1, J hi , and of Algorithm 2, J mf [µ], are each the product of the expected simulation time and the variability of the corresponding likelihood-free weighting.
In the case of Algorithm 2, the performance depends explicitly on the free choice of the function µ(θ, y lo ) that determines the conditional mean of the Poisson-distributed number of high-fidelity simulations required at each iteration. We observe from the first factor in Equation (12b) that, when µ is smaller, the total simulation cost is less. However, the second factor of Equation (12b) implies that as µ decreases, the variability of the likelihood-free weighting can increase without bound, which can severely damage the performance. Thus, Equation (12b) illustrates the characteristic multifidelity trade-off between reducing simulation burden while also controlling the increase in sample variance. Using classical results from calculus of variations, it is possible to determine the mean function, µ , that achieves optimal performance of Algorithm 2, in the sense of minimising the functional, J mf [µ].
Lemma 6. The functional J mf [µ] quantifying the performance of Algorithm 2 is optimised by the function µ , where The function µ given by Lemma 6 defines the optimal number of high-fidelity simulations required in any iteration of Algorithm 2, on average, given the parameter value and low-fidelity simulation output. We note that larger values of η(θ, y lo ) = E (ω hi − ω lo ) 2 | θ, y lo , quantifying the expected squared error between ω hi and ω lo , lead to larger values for µ . Intuitively, if the expected error between the likelihood-free weightings at each fidelity is larger, then the requirement to simulate the high-fidelity model should be greater, to reduce the sample variance. Conversely, where c hi (θ, y lo ) is larger, the greater simulation time of the high-fidelity model means that µ should be smaller, and the requirement for the most expensive simulations is reduced. Intuitively, µ acts to balance the trade-off between controlling simulation cost and variance identified in Equation (12b) above.
It follows from Theorem 5 and Lemma 6 that Algorithm 2 can only ever be an improvement over Algorithm 1 if the optimal performance, J mf = J mf [µ ], satisfies J mf < J hi .
Corollary 7. There exists a mean function, µ, such that the performance of Algorithm 2 exceeds the performance of Algorithm 1, if and only if The first term in Equation (14) justifies the key assumption that the average computational cost of the low-fidelity model is as small as possible compared to that of the high-fidelity model,c lo <c hi . The second term is a measure of the total detriment to the performance of Algorithm 2 incurred by the inaccuracy of ω lo versus ω hi as a Monte Carlo estimate of L hi , as quantified by the function η(θ, y lo ) = E((ω hi − ω lo ) 2 | θ, y lo ). This condition justifies two key criteria for the success of the multifidelity method: that low-fidelity simulations are significantly cheaper than high-fidelity simulations, and that the likelihood-free weightings, ω hi and ω lo , agree sufficiently well, on average. A more detailed analysis of Equation (14) is given in Appendix A.3.
The proofs of Theorem 5, Lemma 6 and Corollary 7 are given in Appendix A. However, we note that these analytical results are useful only insofar as the various functions and constants given in Equation (12) are known. In particular, evaluating the optimal mean function, µ (θ, y lo ), in Algorithm 2 requires knowledge of the functions ∆ q (θ), η(θ, y lo ), and c hi (θ, y lo ), and of the constants V mf andc lo . Similarly, certifying whether the multifidelity approach can outperform the high-fidelity importance sampling method relies on knowing and integrating these functions. However, it is typically the case that these functions and constants are unknown a priori, and need to be estimated based on simulations. In the following section, we describe how the analytical results of Section 3.3 can be used to construct a heuristic adaptive multifidelity algorithm that learns a near-optimal mean function, µ(θ, y lo ), as simulations at each fidelity are completed.

Multifidelity implementation
In Section 3.3, we derived the optimal mean function for the Poisson distribution of the number, M , of high-fidelity simulations required in an iteration of Algorithm 2, conditioned on the parameter value, θ, and low-fidelity simulation, y lo . The optimality condition was based on minimising the functional J mf [µ] defined in Equation (12b), with minimiser µ given in Equation (13). While we can derive the analytical form of µ , this cannot generally be determined a priori, but must be learned in parallel with carrying out likelihood-free inference.
In this section, we describe a practical approach to determining a near-optimal mean function for use in multifidelity likelihood-free inference. We rely on two approximations, relative to the analytically optimal mean function µ given in Equation (13). First, we constrain the optimisation of J mf to the space of functions, µ D , that are piecewise constant in an arbitrary, given, finite partition, D, of the global space of (θ, y lo ) values. The resulting optimisation problem is therefore finite-dimensional. However, although this optimisation can be solved analytically, we can observe that its estimation, being based on the ratios of simulation-based Monte Carlo estimates, is numerically unstable. This motivates a second approximation, which is to follow a gradient-descent approach to allow the mean function to adaptively converge towards the optimum.

Piecewise constant assumption
We constrain the space of mean functions, µ, to be piecewise constant. Consider an arbitrary, given collection D = {D k | k = 1, . . . , K} of ρ-integrable sets that partition the global space of (θ, y lo ) values. We denote a D-piecewise constant function, parametrised by the vector ν = (ν 1 , . . . , ν k ), as Substituting this function into Equation (12b), we can quantify the performance of Algorithm 2, using the mean defined by µ D (θ, y lo ; ν), as the parametrised product, with coefficients given by the integrals Similarly to the functional J mf [µ] optimised in Lemma 6 across positive functions, µ, we can optimise the function J D (ν) across positive vectors, ν.
Lemma 8. The optimal function values ν k for µ D (θ, y lo ; ν) that minimise J D (ν) are for c k and V k as defined in Equation (15) and for V mf andc lo as given in Equations (12e) and (12f). The resulting performance of Algorithm 2 with this mean function is Similarly to the analytical results of Section 3.3, to evaluate ν we need to estimate the values of c lo , c k , V mf and V k , which are unknown a priori. Although these values can be estimated based on Monte Carlo simulation, the rational form of ν k means that these estimates can be unstable, particularly for sets D k ∈ D with small volume (measured by the density, ρ). We now consider a conservative approach to determining values for ν that will provide stable estimates of ν .

Adaptive multifidelity likelihood-free inference
Rather than directly targeting ν , based on ratios of highly variable Monte Carlo estimates, we can introduce a gradient-descent approach to updating the vector ν. Taking derivatives of J D with respect to log ν k for k = 1, . . . , K gives the gradient, Thus, if we write ν (r) for the value of ν used in iteration r of Algorithm 2, we intend to update to ν (r+1) in the next iteration using gradient descent, such that log ν Note that we express this updating rule in terms of log ν (r) k to ensure that each ν (r) k is positive, since the updates to ν (r) are multiplicative. As is typical of gradient-descent approaches, Equation (17) requires the specification of the step size hyperparameter, δ. It is straightforward to show that ν is the unique positive stationary point of Equation (17). Furthermore, since each derivative ∂J D /∂ log ν k is quadratic in the variables c j , V j ,c lo and V mf , the numerical instability in estimating Equation (16) as a ratio does not occur when estimating these derivatives. In relatively undersampled regions D k ∈ D with small ρ-volume, the small values of c k and V k ensure that the convergence to the corresponding estimated optimal value, ν k , is more conservative.
We now explicitly set out the Monte Carlo estimates of c j , V j ,c lo and V mf . These estimates can then be substituted into Equation (17) to produce an updating rule for ν (r) k . This adaptive approach is then implemented into multifidelity likelihood-free importance sampling, as described in Algorithm 3.
Lemma 9. Suppose that r iterations of Algorithm 2 have been completed. We denote: • c lo,i for the observed simulation cost of each y lo,i ; • ω lo,i = ω lo (θ i , y lo,i ) for the low-fidelity weighting calculated at each iteration; • µ i = µ (i) (θ i , y lo,i ) for the value of the mean function used to specify the random variable M in iteration i; • m i for the randomly drawn value of M in iteration i, with mean µ i ; • c hi,i,j for the observed simulation cost of each y hi,i,j for j = 1, . . . , m i as the values of the m i high-fidelity weightings calculated in iteration i, noting that m i may be zero; • ω hi,i,j = ω hi (θ i , y hi,i,j ) for j = 1, . . . , m i as the values of the m i high-fidelity weightings calculated in iteration i; •Ḡ (r) hi as the current Monte Carlo estimate ofḠ hi ; hi π(θ i )/q(θ i ) as the importance weighting, centred on the current estimate ofḠ hi .

The simulation-based Monte Carlo quantities
are consistent estimates ofc lo , V mf , c k , and V k , respectively. (18) into the updating rule, Equation (17), we use log ν

Substituting the estimates in Equation
to update ν (r) to ν (r+1) in adaptive multifidelity likelihood-free importance sampling, as outlined in Algorithm 3. In addition to the specification of the step size hyperparameter, δ, Algorithm 3 also requires a burn-in phase, N 0 , to initialise the Monte Carlo estimates in Equation (18). The partition, D = {D 1 , . . . , D K }, is also an input into Algorithm 3. We defer an investigation of how to choose this partition to future work. For the purposes of this paper, however, we can heuristically construct a partition, D, by fitting a decision tree. We use the burn-in phase of Algorithm 3, over iterations i ≤ N 0 , and regress the values of against features (θ i , y lo,i ), using the CART algorithm [Hastie et al., 2009]

as implemented in
DecisionTrees.jl. Note that this regression is motivated by the form of the true optimal mean function, µ , given in Equation (13). The resulting decision tree defines a partition, D = {D 1 , . . . , D K }, used to define the piecewise-constant mean function µ D (θ, y lo ; ν) over i > N 0 .

Example: Biochemical reaction network
The following example considers the stochastic simulation of a biochemical reaction motif. Readers unfamiliar with these techniques are referred to detailed expositions by Warne et al. [2019] and Erban and Chapman [2019]. We model the conversion (over time t ≥ 0) of substrate molecules, labelled S, into molecules of a product, P. The conversion of S into P is catalysed by the presence of enzyme molecules, E, which bind with S to form molecules of complex, labelled C. After nondimensionalising units of time and volume, this network motif is represented by three reactions, parametrised by the vector θ = (k 1 , k −1 , k 2 ) of positive parameters, k 1 , k −1 , and k 2 , which define three propensity functions, where the integer-valued variables S(t), E(t), C(t) and P (t) represent the molecule numbers at time t > 0. At t = 0, we assume there are no complex or product molecules, but set positive integer numbers S 0 = 100 and E 0 = 5 of substrate and enzyme molecules, respectively. Given the fixed initial conditions, the parameters in θ are sufficient to specify the dynamics of the model in Equation (21a). The model is stochastic, and induces a distribution, which we denote f (· | θ), on the space of trajectories x : t → (S(t), E(t), C(t), P (t)) of molecule numbers in N 4 over time.

ABC importance sampling
We assume that we cannot calculate the likelihood function, L(θ) = f (y 0 | θ). Instead, we need to use simulations to perform ABC. Given θ, the model in Equation (21) can be exactly simulated using the Gillespie stochastic simulation algorithm, to produce draws y ∼ f (· | θ) from the exact model [Gillespie, 1977, Erban and Chapman, 2019, Warne et al., 2020. We will use the ABC likelihood-free weighting with threshold value = 5 on the Euclidean distance of the simulation from y 0 , such that ω(θ, y) = 1( y − y 0 2 < 5), to define the likelihood-free approximation to the posterior, L ABC (θ) = E(ω | θ). We combine this likelihood-free weighting in Algorithm 1 with a rejection sampling approach, setting the importance distribution q = π equal to the prior.

Multifidelity ABC
The exact Gillespie stochastic simulation algorithm can incur significant computational burden.
In the specific case of the network in Equation (21), if the reaction rates k ±1 are large relative to k 2 , there are large numbers of binding/unbinding reactions S + E ← − → C that occur in any simulation. In comparison, the reaction C −−→ P + E can only fire exactly 100 times. Michaelis-Menten dynamics exploit this scale separation to approximate the enzyme kinetics network motif. We approximate the conversion of substrate into product as a single reaction step, where the time-varying rate of conversion, k MM (t), given by induces the propensity function v MM (t) = k MM (t)S(t). We assume initial conditions of S(0) = S 0 = 100 and P (0) = 0, and fix the parameter E 0 = 5. Thus, the parameter vector, θ = (k 1 , k −1 , k 2 ), again fully determines the dynamics of the low-fidelity model in Equation (22). We write f lo (y lo | θ) as the conditional probability density for the Gillespie simulation of the approximated model in Equation (22), where y lo is the vector of ten simulated time points y lo,n at which 10n product molecules have been produced.
For a biochemical reaction network consisting of R reactions, the Gillespie simulation algorithm is a deterministic transformation of R independent unit-rate Poisson processes, one for each reaction channel. We can couple the models in Equations (21) and (22) by using the same Poisson process for the single reaction in Equation (22) and for the product formation C −−→ P + E reaction of Equation (21) Baker, 2020, Lester, 2019]. Using this coupling approach, we first simulate y lo ∼ f lo (· | θ) from Equation (22). We then produce the coupled simulation y hi ∼ f hi (· | θ, y lo ) from the model in Equation (21), using the shared Poisson process. We set the corresponding likelihood-free weightings to ω hi (θ, y hi ) = I(|y hi − y 0 | < 5), noting that E(ω hi | θ) = L ABC (θ) is the high-fidelity ABC approximation to the likelihood. Figure 1 illustrates the effect of coupling between low-fidelity and high-fidelity models. The five coupled high-fidelity simulations are significantly less variable than the independent high-fidelity simulations, appearing almost coincident in Figure 1. This ensures a large degree of correlation between the coupled likelihood-free weightings, ω hi and ω lo . Thus, coupling ensures that ω lo is a reliable proxy for ω hi for use in multifidelity likelihood-free inference.
We implement Algorithm 3 by setting a burn-in period of N 0 = 10, 000, for which we generate m i ∼ M = Poi(1) high-fidelity simulations at each iteration, i ≤ N 0 . Once the burn-in period is complete, we define the partition D by learning a decision tree through a simple regression, as described in Section 4. For iterations i > N 0 beyond the burn-in period, we set a step size of δ = 10 3 for the gradient descent update in Equation (17).

Results
Algorithm 1 was run four times, setting the stop condition to i = 10, 000, i = 20, 000, i = 40, 000 and i = 80, 000. Similarly, Algorithm 3 was run five times, setting the stop condition to i = 40, 000, i = 80, 000, i = 160, 000, i = 320, 000 and i = 640, 000. Figure 2a shows how the variance in the estimate,Ĝ, varies with the total simulation cost, C tot , shown for each of the two algorithms. The slope of each curve (on a log-log scale) is approximately −1, corresponding to the dominant behaviour of the MSE being reciprocal with total simulation time, as observed in Equation (11).  (21) and (22) for parameters θ = (k 1 , k −1 , k 2 ) = (50, 50, 1), compared with data used for inference. For one low-fidelity simulation, we generate five uncoupled simulations and five coupled simulations. (b) Ten-dimensional data summarising simulated trajectories in (a). Black represents observed data, y 0 ; the single low-fidelity simulation y lo ∼ f lo (· | θ) is in blue; five uncoupled simulations y hi ∼ f hi (· | θ) are in orange; five coupled simulations y hi ∼ f hi (· | θ, y lo ) are in green (almost coincident).
The offset in the two curves corresponds to the inequality J mf < J hi in the leading order coefficient, thereby demonstrating the improved performance of Algorithm 3 over Algorithm 1.
The values in Figure 2b show the multifidelity weights, w i . We show only those weights not equal to zero or one, corresponding to those iterations where ω lo (θ i , y lo,i ) has been corrected by at least one ω hi (θ i , y hi,i,j ) = ω lo (θ i , y lo,i ). Clearly there is a significant amount of correction applied to the lowfidelity weights. However, as demonstrated by the improved performance statistics, Algorithm 3 has learned the required allocation of computational budget to the high-fidelity simulations that balances the trade-off between achieving reduced overall simulation times and correcting inaccuracies in the low-fidelity simulation.
Each run of Algorithm 3 includes a burn-in period of 10, 000 iterations, at the conclusion of which a partition D is created, based on decision tree regression. In Appendix B, we show how this decision tree is used to define a piecewise-constant mean function, specifically for the partition D used for the final run of Algorithm 3 (i.e. for stopping condition i = 640, 000). In Figure 2c, we show the evolution of the values of ν (i) k used in this mean function, over iterations i. Following the updating rule in Equation (19), the trajectory of ν (i) k converges exponentially towards a Monte Carlo estimate of the optimal value ν k given in Equation (16). However, we can see from Figure 2c that, as more simulations are completed and the Monte Carlo estimates in Equation (18) evolve, the values of each parameter, ν k , track updated estimates. This is illustrated in Figure 2d for ν 1 , where the estimated optimum ν 1 evolves as more simulations are completed. We note that the gradient descent update in Equation (19) at iteration i depends on all ν (i) k values. Thus, the observed convergence of ν (i) 1 to the evolving estimate of ν 1 is not necessarily monotonic. Figure 2d illustrates the motivation for the use of gradient descent rather than simply using the analytically obtained optimum. When very few simulations have been completed, then the estimates in Equation (18) are small and their ratios are numerically unstable, and often far from the true optimum. If ν (i) k values are too small in early iterations, then estimates become more numerically unstable, since fewer high-fidelity simulations are completed for small µ. Instead, using gradient descent ensures that enough high-fidelity simulations are completed for each D k , including those with low volume under the measure ρ, to stabilise the estimates in Equation (18) and thus stabilise the multifidelity algorithm.
Algorithm 1 was run three times, using ω hi (θ, y hi ) dependent on high-fidelity simulations y hi ∼ f (· | θ), alone, and setting the stop condition to i = 2, 500, i = 5, 000 and i = 10, 000. Similarly, Algorithm 3 was run four times using the coupled multifidelity model, setting the stop condition to i = 4, 000, i = 8, 000, i = 16, 000 and i = 32, 000, and initialising with a burn-in of size N 0 = 2, 000. The adaptive step size is set to δ = 10 8 . In both algorithms, we set the number of simulations required for each evaluation of ω hi (θ, (y hi,1 , . . . , y hi,K )) or ω lo (θ, (y lo,1 , . . . , y lo,K )) as K = 100. Figure 3 depicts the performance of multifidelity Bayesian synthetic likelihood (BSL) inference, where Algorithm 3 is applied with BSL likelihood-free weightings, ω lo and ω hi . As with MF-ABC, Figure 3a shows that the MF-BSL generates improved performance over high-fidelity BSL inference, achieving lower variance estimates for a given computational budget. We also note in Figure 3a that the curve corresponding to MF-BSL has slope less than −1. This is due to (a) the overhead cost of the initial burn-in period of Algorithm 3, and also (b) the conservative convergence of ν (i) to the optimum, as shown in Figure 3c-d. Both observations imply that earlier iterations are less efficiently produced than later iterations, meaning that larger samples show greater improvements than expected from the reciprocal relationship in Equation (11).
Comparing Figure 3b to Figure 2b, we note that there are very few negative multifidelity weightings in MF-BSL, in comparison to MF-ABC. We can conclude that the Bayesian synthetic likelihood, constructed using low-fidelity simulations, tends to underestimate the likelihood of the observed data compared to using high-fidelity simulations. We note also in this comparison that the multifidelity likelihood-free weightings are on significantly different scales.

Discussion
The characteristic computational burden of simulation-based, likelihood-free Bayesian inference methods is often a barrier to their successful implementation. Multifidelity simulation techniques have previously been shown to improve the efficiency of likelihood-free inference in the context of ABC. In this work, we have demonstrated that these techniques can be readily applied to general likelihood-free approaches. Furthermore, we have introduced a computational methodology for automating the multifidelity approach, adaptively allocating simulation resources across different fidelities in order to ensure near-optimal efficiency gains from this technique. As parameter space is explored, our methodology, given in Algorithm 3, learns the relationships between simulation accuracy and simulation costs at the different fidelities, and adapts the requirement for high-fidelity simulation accordingly.
The multifidelity approach to likelihood-free inference is one of a number of strategies for speeding up inference, which include MCMC and SMC sampling techniques [Marjoram et al., 2003, Sisson et al., 2007, Toni et al., 2009 and methods for variance reduction such as multilevel estimation [Giles, 2015, Guha and Tan, 2017, Warne et al., 2018, Jasra et al., 2019. A key observation in the previous work of Prescott and Baker [2021] and Warne et al. [2021b] is that applying multifidelity techniques provides 'orthogonal' improvements that combine synergistically with these other established approaches to improving efficiency. Similarly, we envision that Algorithm 3 can be adapted into an SMC or multilevel algorithm with minimal difficulty, following the templates set by Prescott and Baker [2021] and Warne et al. [2021b].
The multifidelity approach discussed in this work is a highly flexible generalisation of existing multifidelity techniques, which can be viewed as special cases of Algorithm 2. In each of MF-ABC Baker, 2020, 2021], LZ-ABC [Prangle, 2016], and DA-ABC [Everitt and Rowińska, 2021], it is assumed that ω hi is an ABC likelihood-free weighting, which we relax in this work. Furthermore, LZ-ABC and DA-ABC both use ω lo ≡ 0, so that parameters are always rejected if no high-fidelity simulation is completed. Clearly, we relax this assumption to allow for any low-fidelity likelihood-free weighting. In all of MF-ABC, LZ-ABC and DA-ABC, the conditional distribution of M , given a parameter value θ and low-fidelity simulation output y lo is Bernoulli distributed, with mean µ(θ, y lo ) ∈ (0, 1]. In this work we change this distribution to Poisson, to ease analytical results, but any conditional distribution for M can be used. These adaptations are explored further in Appendix A. In the case of MF-ABC (as originally formulated by Prescott and Baker [2020]) and DA-ABC [Christen andFox, 2005, Everitt andRowińska, 2021], the mean function, µ(θ, y lo ), depends on a single low-fidelity simulation and is assumed to be piecewise constant in the value of the indicator function 1(d(y lo , y 0 ) < ). LZ-ABC is more generic in its definition of µ = µ(φ(θ, y lo )) to depend on the value of any decision statistic, φ. In this work, we consider more general piecewise constant mean functions, µ D , for heuristically derived partitions D of (θ, y lo )-space. We observe that (θ, y lo ) may be of very high dimension; in the BSL example in Section 5.2, having K = 100 low-fidelity simulations y lo ∈ R 10 means that the input to µ is of dimension 1003. In this situation, it may be tempting to seek a mean function that only depends on θ. However, we recall that the optimal mean function, µ (θ, y lo ), derived in Lemma 6, depends on the conditional expectation E((ω hi − ω lo ) 2 | θ, y lo ). Thus, by ignoring y lo , we would ignore the information about ω hi given by the evaluation of ω lo (θ, y lo ). Furthermore, the high dimension of the inputs to µ suggest that this function is not necessarily well-approximated by a decision tree. Future work may focus on methods to learn the optimal mean function directly without resorting to piecewise constant approximations [Levine and Stuart, 2021]. The key problem is ensuring the conservatism of any alternative estimate of µ , recalling that the variance of w mf is inversely proportional to µ.
In the example explored in Section 5, we considered the use of Algorithm 3 where ω hi and ω lo were first both ABC likelihood-free weightings, and then both BSL likelihood-free weightings.
In principle, this method should also allow for ω lo to be, for example, an ABC likelihood-free weighting based on a single low-fidelity simulation, and ω hi to be a BSL likelihood-free weighting based on K > 1 high-fidelity simulations. However, the success of the multifidelity method depends explicitly on the function η(θ, y lo ) = E((ω hi − ω lo ) 2 | θ, y lo ) being sufficiently small, as quantified in Corollary 7. If ω lo and ω hi are on different scales, as is likely when one is an ABC weighting and one a BSL weighting, then this function is not sufficiently small in general, and so the multifidelity approach fails. We note, however, that we could instead consider the scaled low-fidelity weighting, ω lo = γω lo , in place of ω lo in Algorithms 2 and 3 with no change to the target distribution. Here, γ is an additional parameter that can be tuned with µ when minimising the performance metric, J mf ; the optimal value of this parameter would need to be learned in parallel with the optimal mean function, µ. We defer this adaptation to future work.
Finally, this work follows Baker [2020, 2021] in considering only a single low-fidelity model. There is significant scope for further improvements by applying these approaches to suites of low-fidelity approximations [Gorodetsky et al., 2021]. For example, exact stochastic simulations of biochemical networks, such as that simulated in Section 5, may also be approximated by tauleaping [Gillespie, 2001, Warne et al., 2019, where the time discretisation parameter τ tends to be chosen to trade off computational savings against accuracy: exactly the trade-off explored in this work. Clearly, this parameter therefore has important consequences for the success of a multifidelity inference approach using such an approximation strategy. More generally, a full exploration of the use of multiple low-fidelity model approximations will be vital for the full potential of multifidelity likelihood-free inference to be realised.
A Analytical results: Comparing performance A.1 Theorem 5 Proof. The leading order performance of each of Algorithm 1 and Algorithm 2 is given in terms of increasing computational budget, C tot , in Equation (10) and Equation (11), respectively. For the performance of Algorithm 2 to exceed that of Algorithm 1, we compare the leading order coefficients from Equations (10) and (11), requiring We note that E(w mf | θ) = π(θ)L mf (θ)/q(θ) and E(w hi | θ) = π(θ)L hi (θ)/q(θ). Since L mf = L hi , as shown in Proposition 4, the denominators in Equation (23) are therefore equal. Thus, is the condition for Algorithm 2 to outperform Algorithm 1.
Taking the right-hand side of this inequality first, clearly the expected simulation time is E(C hi ) = c hi , for the constantc hi defined in Equation (12c). Similarly, we can write as given in Equation (12d). Thus, J hi =c hi V hi .
For the left-hand side of the performance inequality, we take each expectation in J mf in turn. We first note that the expected iteration cost of Algorithm 2, E(C mf ), is the sum of the expected cost of a single low-fidelity simulation, and the expected cost of M high-fidelity simulations. By definition, the expected cost of a single low-fidelity simulation y lo ∼ f lo (· | θ) across θ ∼ q(·) is given byc lo . Thus the remaining cost, E(δC mf ) = E(C mf ) −c lo , is the expected cost of M high-fidelity simulations. Conditioning on θ, y lo and M = m, the expected remaining cost is, by definition, E(δC mf | θ, y lo , M = m) = mc hi (θ, y lo ).
Finally, integrating this expression over the density ρ in Equation (12j) gives the first factor of Equation (12b).
It remains to show that We first condition on θ, y lo and M = m, to write where we exploit the conditional independence of the high-fidelity simulations y hi,i and y hi,j , for i = j. On substitution of these conditional expectations, we then rearrange to write where we write the conditional expectation λ hi (θ, y lo ) = E(ω hi | θ, y lo ). At this point we can take expectations over M and rearrange to give Here, we can use the assumption that M conditioned on θ and y lo is Poisson distributed, noting that the statement of Theorem 5 can be adapted for other conditional distributions of M with different conditional variance functions. Under the Poisson assumption, we can substitute Var(M | θ, y lo ) = µ(θ, y lo ) to give E(w 2 mf ∆ 2 | θ, y lo ) = ∆ 2 q λ 2 hi + E (ω hi − ω lo ) 2 | θ, y lo µ(θ, y lo ) .
Finally, we take expectations with respect to the probability density ρ in Equation (12j), and the product in Equation (12b) follows.
A.1.1 Alternative conditional distributions for M The proof above derives the performance measure J mf given in Equation (12b), under the assumption that the conditional distribution of M , given θ and y lo , is Poisson with mean µ(θ, y lo ). The following corollaries adapt the expression for J mf in the case of alternative conditional distributions for M . We first define the MSE, E mf = [λ hi (θ, y lo ) − ω lo (θ, y lo )] 2 ρ(θ, y lo )dθdy lo , between ω lo (θ, y lo ) and λ hi (θ, y lo ) = E(ω hi | θ, y lo ).

A.2 Lemma 6
We return to the assumption that M is conditionally Poisson distributed, given θ and y lo .
Proof. We write the functional J mf [µ] = C[µ]V [µ] in Equation (12b) as the product of functionals, Standard 'product rule' results from calculus of variations allows us to write the functional derivative of J mf with respect to µ as Setting this functional derivative to zero, the optimal function, µ , satisfies ∆ q (θ) 2 η(θ, y lo ) c hi (θ, y lo ) .
The result in Equation (13)
A simple rearrangement of this inequality gives the inequality in Equation (14).
These equalities imply that E(ω hi | θ) 2 ≤ E(E(ω hi | θ, y lo ) 2 | θ) ≤ E(ω 2 hi | θ), where the lower bound is achieved for y hi independent of y lo , while the upper bound would be achieved if y hi were a deterministic function of y lo . In particular, V mf /V hi ≤ 1, and so the first term of Equation (14) is small whenever the low-fidelity model provides significant computational savings versus the high-fidelity model.
The second term in Equation (14) quantifies the detriment to the performance of Algorithm 2 that arises from the inaccuracy of ω lo as an estimate of ω hi . The function η(θ, y lo ) = E((ω hi − ω lo ) 2 | θ, y lo ) is integrated across the density ρ, weighted by the relative computational cost of the highfidelity simulation, c hi (θ, y lo )/c hi , and by the contribution of G(θ) to the variance of the estimated posterior expectation of G. We can conclude that the multifidelity approach requires that the low-fidelity model is accurate in the regions of parameter space where high-fidelity simulations are particularly expensive.
To summarise: if (a) the ratio between average low-fidelity simulation costs and high-fidelity simulation costs is suitably small, and (b) the average disagreement between likelihood-free weightings, as measured by η, is suitably small, then Equation (14) will be satisfied and thus a mean function, µ , exists such that Algorithm 2 is more efficient than Algorithm 1.

B Mean functions
Algorithm 4 Piecewise constant mean function µ D (θ, y lo ; ν) used in MF-ABC Algorithm 3, depicted in Figure 2c, at final iteration.