The frontier of simulation-based inference

Simulation-based inference

Statistical inference is performed within the context of a statistical model, and in simulation-based inference the simulator itself defines the statistical model. For the purpose of this paper, a simulator is a computer program that takes as input a vector of parameters $\theta$, samples a series of internal states or latent variables $z_{i}\sim p_{i}(z_{i}|\theta,z_{<i})$, and finally produces a data vector $x\sim p(x|\theta,z)$ as output. Programs that involve random samplings and are interpreted as statistical models are known as probabilistic programs, and simulators are an example. Within this general formulation, real-life simulators can vary substantially:

The parameters $\theta$ describe the underlying mechanistic model and thus affect the transition probabilities $p_{i}(z_{i}|\theta,z_{<i})$. Typically the mechanistic model is interpretable by a domain scientist and $\theta$ has relatively few components and a fixed dimensionality. Examples include coefficients found in the Hamiltonian of a physical system, the virulence and incubation rate of a pathogen, or fundamental constants of Nature.

The latent variables $z$ that appear in the data-generating process may directly or indirectly correspond to a physically meaningful state of a system, but typically this state is unobservable in practice. The structure of the latent space varies substantially between simulators. The latent variables may be continuous or discrete and the dimensionality of the latent space may be fixed or may vary depending on the control flow of the simulator. The simulation can freely combine deterministic and stochastic steps. The deterministic components of the simulator may be differentiable or may involve discontinuous control flow elements. In practice, some simulators may provide convenient access to the latent variables, while others are effectively black boxes. Any given simulator may combine these different aspects in almost any way.

Finally, the output data $x$ correspond to the observations. They can range from a few unstructured numbers to high-dimensional and highly structured data, such as images or geospatial information.

Consider for instance the systems shown in Fig. 1. Particle physics processes often only depend on a small number of parameters of interest such as particle masses or coupling strengths. The latent process combines a high-energy interaction, rigorously described by a quantum field theory, with the passage of the resulting particles through an incredibly complex detector, most accurately modeled with stochastic simulations with billions of latent variables; this second part often does not depend on the parameters of interest. The output data consist, in their raw form, of millions of sensor read-outs, though there is an established pipeline that compresses this raw data to tens to hundreds of observables. Epidemiological simulations can be based on a network structure with geospatial properties, and the latent process consists of many repeated structurally identical stochastic time steps. In contrast, cosmological simulations of the evolution of the Universe may consist of a highly structured stochastic initial state followed by a smooth, deterministic time evolution.

These differences mean that there is no one-size-fits-all inference method. In this review we aim to clarify the considerations needed to choose the most appropriate approach for a given problem.

2 Inference

Scientific inference tasks differ by what is being inferred: given observed data $x$, is the goal to infer the input parameters $\theta$, or the latent variables $z$, or both? Sometimes only a subset of the parameters (or latent variables) are of interest, while the rest are nuisance parameters (i. e. parameters that we are not directly interested in but must account for because they influence the distributions of the data). We will focus on the common problem of inferring $\theta$ in a parametric setting, we will comment on methods that allow inference on $z$, and we will not focus on non-parametric inverse problems.

Inference may be performed either in a frequentist or a Bayesian approach and may be limited to point estimates $\hat{\theta}(x)$ or extended to include a probabilistic notion of uncertainty. In the frequentist case, confidence sets are often formed from inverting hypothesis tests, based on the likelihood ratio test statistic. In Bayesian inference, the goal is typically to calculate the posterior

for observed data $x$ and a given prior $p(\theta)$. In both cases the likelihood function $p(x|\theta)$ is a key ingredient.

The fundamental challenge for simulation-based inference problems is that the likelihood function $p(x|\theta)$ implicitly defined by the simulator is typically not tractable, as it corresponds to an integral over all possible trajectories through the latent space, i. e. all possible execution traces of the simulator. That is,

where $p(x,z|\theta)$ is the joint probability density of data $x$ and latent variables $z$. For a simple sequential data generation procedure, the joint likelihood can be written as $p(x,z|\theta)=p(x|\theta,z)\prod_{i}p_{i}(z_{i}|\theta,z_{<i})$. For real-life simulators with large latent spaces, it is clearly impossible to compute this integral explicitly. Since the likelihood function is the central ingredient to both frequentist and Bayesian inference, this is a major challenge for inference in many fields. This paper reviews simulation-based or likelihood-free inference techniques that enable frequentist or Bayesian inference despite this intractability. These methods can be seen as a specialization of inverse Uncertainty Quantification (UQ) on the model parameters in situations with accurate, stochastic simulators.

There is a second, more widely appreciated source of intractability. In the case of Bayesian inference, the evidence—the denominator of (1)—involves an integral over the parameters $\theta$. In problems with high-dimensional parameters this becomes intractable, independently of the intractability of the likelihood function. This challenge is commonly addressed with Markov Chain Monte Carlo (MCMC) methods (10, 11) or variational inference (VI) (12).

3 Traditional methods

The problem of inference without tractable likelihoods is not a new one, and two major approaches have been developed to address it. Arguably the most well-known is Approximate Bayesian Computation (ABC) (2, 3). Until recently, it was so established that the terms “likelihood-free inference” and “ABC” were often used interchangeably. In the simplest form of rejection ABC, the parameters $\theta$ are drawn from the prior, the simulator is run with those values to sample $x_{\text{sim}}\sim p(\cdot|\theta)$, and $\theta$ is retained as posterior sample if the simulated data is sufficiently close to the observed data. In essence, the likelihood is approximated by the probability that the condition $\rho(x_{\text{sim}},x_{\text{obs}})<\epsilon$ is satisfied, where $\rho$ is some distance measure and $\epsilon$ is a tolerance. The accepted samples then follow an approximate version of the posterior. We show a schematic workflow of this algorithm in Fig. 3a (for a more elaborate Markov Chain Monte Carlo algorithm with a proposal function).

In the limit $\epsilon\to 0$, inference with ABC becomes exact, but for continuous data the acceptance probability vanishes. In practice, small values of $\epsilon$ require unfeasibly many simulations. For large $\epsilon$, sample efficiency is increased at the expense of inference quality. Similarly, the sample efficiency of ABC scales poorly to high-dimensional data $x$. Since the data immediately affect the rejection process (and in more advanced ABC algorithms the proposal distribution), inference for new observations requires repeating the entire inference algorithm. ABC is thus best-suited for the case of a single observation or at most a few i.i.d. data points.

Lacking space to do the vast ABC literature justice, we refer the reader to a review of ABC methods, see Ref. (13), and highlight the combination with MCMC (14) and Sequential Monte Carlo (SMC) (15, 16).

The second classical approach to simulation-based inference is based on creating a model for the likelihood by estimating the distribution of simulated data with histograms or kernel density estimation (1). Frequentist and Bayesian inference then proceeds as if the likelihood were tractable. We sketch this algorithm in Fig. 3e (replacing the green learning step with a classical density estimation method). This approach has enough similarities to ABC to be dubbed “Approximate Frequentist Computation” by the authors of Ref. (17). One advantage over ABC is that it is amortized: after an upfront computational cost at the simulation and density estimation stage, new data points can be evaluated efficiently. In Fig. 3e this manifests itself as the blue “data” box only entering at the inference stage and not affecting the expensive simulation step. This property makes density estimation–based inference particularly well-suited for problems with many i.i.d. observations, a key reason for its wide-spread use in particle physics measurements.

Both of the traditional approaches suffer from the curse of dimensionality: in the worst case, the required number of simulations increases exponentially with the dimension of the data $x$. Therefore both approaches rely on low-dimensional summary statistics $y(x)$ and the quality of inference is tied to how well those summaries retain information about the parameters $\theta$. Traditionally, the development of powerful summary statistics has been the task of a domain expert, and the summary statistics have been prescribed prior to inference.

Frontiers of simulation-based inference

These traditional simulation-based inference techniques have played a key role in several fields for years. However, they suffer from shortcomings in three crucial aspects:

Sample efficiency: Both ABC and classical density estimation techniques suffer from the curse of dimensionality. The poor scaling means that the number of simulated samples needed to provide a good estimate of the likelihood or posterior can be prohibitively expensive.

Quality of inference: The reduction of the data to low-dimensional summary statistics invariably discards some of the information in the data about $\theta$, which results in a loss in statistical power. Large values of the $\epsilon$ parameter in ABC or bandwidth parameter for kernel density estimation lead to poor approximations of the true likelihood. Both reduce the overall quality of inference.

Amortization: Performing inference with ABC for a new set of observed data requires repeating most steps of the inference chain, in particular if the proposal distribution depends on the observed data. The method scales poorly when applied to large numbers of observations. On the other hand, inference based on density estimation is amortized: the computationally expensive steps do not have to repeated for new observations. This is particularly desirable for the case with i.i.d. observations.

In recent years, new capabilities have become available that let us improve all three of these aspects. We loosely group them into three main directions of progress:

The revolution in machine learning allows us to work with higher-dimensional data, which can improve the quality of inference. Inference methods based on neural network surrogates are directly benefitting from the impressive rate of progress in deep learning.

Active learning methods can systematically improve sample efficiency, letting us tackle more computationally expensive simulators.

The deep integration of automatic differentiation and probabilistic programming into the simulation code, as well as the augmentation of training data with additional information that can be extracted from the simulator, are changing the way the simulator is treated in inference: it is no longer a black box, but exposed to the inference workflow.

We sketch these trends in Fig. 2, broadly categorizing the landscape of inference tasks in a two-dimensional plane of the dimensionality of the data (vertical axis) and the complexity of the simulator (horizontal axis).

Over the last decade, machine learning techniques, in particular deep neural networks, have turned into versatile, powerful, and popular tools for a variety of problems (18). Neural networks initially demonstrated breakthroughs in supervised learning tasks such as classification and regression. They can easily be composed to solve higher-level tasks, lending themselves to problems with a hierarchical or compositional structure. Architectures tailored to various data structures have been developed, including dense or fully connected networks aimed at unstructured data, convolutional neural networks that leverage spatial structure for instance in image data, recurrent neural networks for variable-length sequences, and graph neural networks for graph-structured data. Choosing an architecture well suited for a specific data structure is an example of inductive bias, which more generally refers to the assumptions inherent in a learning algorithm independent of the data. Inductive bias is one of the key ingredients behind most successful applications of deep learning, though it is difficult to characterize its role precisely.

Another class of approaches use autoregressive models, in which the probability distribution of a high-dimensional variable is factorized into successive conditional densities of the individual components (22, 23, 24, 25, 26, 27, 28, 29, 30). These models are expressive, have a tractable (conditional) density, and can be used to generate synthetic data. While autoregressive models are somewhat disfavored in industrial applications because generating samples from them can be slow, the sequential nature is more closely aligned with the way simulators are written and offers an opportunity to align the neural networks latent variables with the semantically meaningful latent variables of simulators.

Generative Adversarial Networks (GANs) are an alternative type of generative model based on neural networks. Unlike in normalizing flows and autoregressive models, the transformation implemented by the generator is not restricted to be invertible. While this allows for more expressiveness, the density defined by the generator is intractable. Since maximum likelihood is not a possible training objective, the generator is pitted against an adversary, whose role is to distinguish the generated data from the target distribution. We will later discuss how the same idea can be used for simulation-based inference, using an idea known as the “likelihood ratio trick”.

2 Active learning

A simple, but very impactful idea is to run the simulator at parameter points that are expected to increase our knowledge the most. This can be done iteratively such that after each simulation the knowledge resulting from all previous runs is used to guide which parameter point should be used next. There are multiple technical realizations of this idea of active learning. It is commonly applied in a Bayesian setting, where the posterior can be continuously updated and used to steer the proposal distribution of simulator parameters (35, 36, 37, 38, 39, 40, 41). But it applies equally well to efficiently calculating frequentist confidence sets (42, 43, 44). Even simple implementations can lead to a substantial improvement in sample efficiency.

Similar ideas are discussed in the context of decision making, experimental design, and reinforcement learning, and we expect further improvements in inference algorithms from the cross-pollination between these fields. For instance, a question that is occasionally discussed in the context of reinforcement learning (45, 46) or Bayesian optimization (47), but has not yet been applied to the likelihood-free setting, is how to make use of multi-fidelity simulators offering multiple levels of precision or approximations.

3 Integration and augmentation

Both machine learning and active learning can substantially improve quality of inference and sample efficiency compared to classical methods. However, overall they do not change the basic approach to simulation-based inference dramatically: they still treat the simulator as a generative black box that takes parameters as input and provides data as output, with a clear separation between the simulator and the inference engine. A third direction of research is changing this perspective, opening this black box to access more information and integrating inference and simulation more tightly.

One example of this shift is the probabilistic programming paradigm. Gordon et al. (48) describe probabilistic programs as the usual functional or imperative programs with two added constructs: (1) the ability to draw values at random from distributions, and (2) the ability to condition values of variables in a program via observations. We have already described simulators as probabilistic programs focusing on the first construct, which does not require opening the black box. However, conditioning on the observations requires a deeper integration as it involves controlling the randomness in the generative process. This approach abstracts the capabilities needed to implement Particle Filters and SMC (49). Previously, this required writing the program in a special purpose language; however, recent work allows these capabilities to be added to existing simulators with minimal changes to their codebase (50). Ultimately, probabilistic programming aims at providing the tools to infer the incredibly complex space of all execution traces of the simulator conditioned on the observation.

A complementary development is the observation that additional information that characterizes the latent data-generating process can be extracted from the simulator and used to augment the data used to train surrogates. This augmented training data can be exploited in supervised learning objectives, and can dramatically increase the sample efficiency for surrogate training. Those developing inference algorithms and those familiar with the details of the simulator should consider whether, in addition to the sole ability to sample $x\sim p(x|\theta)$, the following properties in which the corresponding quantities in the simulator are well-defined and tractable.

$p(x|z,\theta)$: the probability density of the output data given the latent variables.

$t(x,z|\theta)\equiv\nabla_{\theta}\log p(x,z|\theta)$: the joint score is the gradient of the joint log probability density of output data and latent variables with respect to the parameters.

$\nabla_{z}\log p(x,z|\theta)$: the gradient of the same quantity, but with respect to the latent variables.

$r(x,z|\theta,\theta^{\prime})\equiv p(x,z|\theta)/p(x,z|\theta\textquoteright)$: the ratio of the joint probability density of output data and latent variables for two parameter points $\theta$ and $\theta^{\prime}$.

$\nabla_{\theta}(x,z)$: the derivative of output data and latent variables with respect to the parameters.

$\nabla_{z}x$: the gradient of the output data with respect to the latent variables.

These quantities can then be used to augment the usual output $x$ from the simulator and can be exploited in supervised learning objectives, and can dramatically increase the sample efficiency for surrogate training (51, 52, 53), as we will detail later.

Many of the quantities above involve derivatives, which can now be efficiently calculated using automatic differentiation (often referred to simply as autodiff) (54). Autodiff is a family of techniques similar to but more general than the backpropagation algorithm that is ubiquitous in deep learning. Automatic differentiation, like probabilistic programming, involves non-standard interpretations of the simulation code and has been developed by a small but established field of computer science. In recent years several have advocated that deep learning would be better described as differential programming (55, 56). With this view, incorporating autodiff into existing simulation codes is a more direct way to exploit the advances in deep learning than trying to incorporate domain knowledge into an entirely foreign substrate such as a deep neural network.

Extracting the necessary information from the simulator again requires integration deep in the code. While technologies to incorporate probabilistic programming paradigm into existing code bases are just emerging, the development of tools to enable autodiff in the most commonly used scientific programming languages is well advanced. We highlight that two of the quantities above (II and III) involve both autodiff and probabilistic programming. The integration of inference and simulation as well as the idea of augmenting the training data with additional quantities have the potential to change the way we think about simulation-based inference. In particular, this perspective can influence the way simulation codes are developed in order to provide these new capabilities.

Workflows for simulation-based Inference

This wide array of capabilities can be combined in different inference workflows. Some of these are structurally identical to the traditional ABC and density estimation–based methods, while others are fundamentally different. As a guideline through this array of different workflows, let us first discuss common building blocks, and the different approaches that can be taken in each of these components. In Fig. 3 and the following sections we will then piece these blocks together into different inference algorithms.

An integral part of all inference methods is running the simulator, visualized as a yellow pentagon in Fig. 3. The parameters at which the simulator is run are drawn from some proposal distribution, which may or may not depend on the prior in a Bayesian setting, and can either be chosen statically, or iteratively with an active learning method. Next, the potentially high-dimensional output from the simulator may be used directly as input to the inference method, or reduced to low-dimensional summary statistics, which may be prescribed or learned from data.

The inference techniques can be broadly separated into those which, like ABC, use the simulator itself during inference, and methods which construct a surrogate model and use that for inference. In the first case, the output of the simulator is directly compared to data, see the top panels of Fig. 3. In the latter case, the output of the simulator is used as training data for an estimation or machine learning stage, shown as the green boxes in the bottom panels of Fig. 3. The resulting surrogate models, shown as red hexagons, are then used for inference.

The algorithms address the intractability of the true likelihood in different ways: some methods construct a tractable surrogate for the likelihood function, others for the likelihood ratio function, both of which make frequentist inference straightforward. In other methods, the likelihood function never appears explicitly, for instance when it is implicitly replaced by rejection probability (an approach that does not lend itself to frequentist inference).

The final target for Bayesian inference is the posterior. Methods differ in whether they provide access to samples of parameter points sampled from the posterior, for instance from MCMC or ABC, or a tractable function that approximates the posterior function. Similarly, some methods require specifying which quantities are to be inferred early on in the workflow, while others allow this decision to be postponed.

Let us now discuss how these blocks and computational capabilities can be combined into inference techniques, beginning with those which, like ABC, use the simulator directly during inference. We sketch some of these algorithms in the top panels of Fig. 3.

One of the major shortcomings of ABC is its reliance on low-dimensional summary statistics. Classifier ABC (57) removes the requirement of compressing the data into summary statistics by instead training a classifier to estimate the discrepancy between observed and simulated data.

A reason for the poor sample efficiency of the original rejection ABC algorithm is that the simulator is run at parameter points drawn from the prior, which may have a large mass in regions that are in strong disagreement with the observed data. Different algorithms have been proposed that instead run the simulator at parameter points that are expected to improve the knowledge on the posterior the most (35, 36, 37, 38, 39). Compared to vanilla ABC, these techniques improve sample efficiency, though they still require the choice of summary statistics, distance measure $\rho$, and tolerance $\epsilon$.

In the case where the final stage of the simulator is tractable or the simulator is differentiable (respectively, properties I and VI from the list in Sec. 2.2.3), asymptotically exact Bayesian inference is possible (51) without relying on a distance tolerance or summary statistics, removing ABC’s main limitations in terms of quality of inference.

The probabilistic programming paradigm presents a more fundamental change to how inference is performed. First, it requires the simulator to be written in a probabilistic programming language, though recent work allows these capabilities to be added to existing simulators with minimal changes to their codebase (50). In addition, probabilistic programming either requires a tractable likelihood for the final step $p(x|z,\theta)$ (quantity I) or the introduction of an ABC-like comparison. When these criteria are satisfied, several inference algorithms exist that can draw samples from the posterior $p(\theta,z|x)$ of the input parameters $\theta$ and the latent variables $z$ given some observed data $x$. These techniques are either based on MCMC, see Fig. 3c, or on training a neural network to provide proposal distributions (58) as shown in Fig. 3d. The key difference to ABC is that the inference engine controls all steps in the program execution and can bias each draw of random latent variables to make the simulation more likely to match the observed data, improving sample efficiency.

A strength of these algorithms is that they allow to infer not only the input parameters into the simulator, but the entire latent process leading to a particular observation. This allows us to answer entirely different questions about scientific processes, adding a particular kind of physical interpretability that methods based on surrogates do not possess. While standard ABC algorithms in principle allow for inference on $z$, probabilistic programming solves this task more efficiently.

2 Surrogate models

A key disadvantage of using the simulator directly during inference is the lack of amortization. When new observed data becomes available, the whole inference chain has to be repeated. By training a tractable surrogate or emulator for the simulator, inference is amortized: after a (computationally expensive) upfront simulation and training phase, new data can be evaluated very efficiently. This approach scales particularly well to data consisting of many i.i.d. observations. As discussed in Sec. 11.3, this is not a new idea, and well-established methods use classical density estimation techniques to create a surrogate model for the likelihood function. But the new computational capabilities discussed in Sec. 2.1 have given new momentum to this class of inference techniques. They can be incorporated in various ways, providing approximations of the true parameters given data, defining suitable summary statistics, learning the likelihood, the likelihood ratio, or the posterior; we will briefly go through these different options to organize the inference. Selected algorithms are visualized in the bottom panels of Fig. 3.

Perhaps the most obvious of these approaches is to directly invert the parameter-to-data process implemented by the simulator and to train a model to estimate the true parameters $\hat{\theta}(x)$ as a function of observed data $x$ (59, 60, 61, 62). However, point estimates are not always useful and a probabilistic interpretation in terms of the likelihood or posterior of these methods is not obvious, so will not focus on this approach.

A powerful probabilistic approach is to train neural conditional density estimators such as normalizing flows as a surrogate for the simulator. The conditional density can be defined in two directions: the network can either learn the posterior $p(\theta|x)$ (63, 64, 40, 65, 66, 67) or the likelihood $p(x|\theta)$ (41, 68, 69, 70). We show these three techniques in Figs. 3e and 3f; note that the likelihood surrogate algorithm is structurally identical to the classical density estimation–based approach, but uses more powerful density estimation techniques.

Relatedly, neural networks can be trained to learn the likelihood ratio function $p(x|\theta_{0})/p(x|\theta_{1})$ or $p(x|\theta_{0})/p(x)$, where in the latter case the denominator is given by a marginal model integrated over a proposal or the prior (71, 72, 73, 4, 74, 57, 75, 76, 77). We sketch this approach in Fig. 3g. The key idea is closely related to the discriminator network in GANs mentioned above: a classifier is trained using supervised learning to discriminate two sets of data, though in this case both sets come from the simulator and are generated for different parameter points $\theta_{0}$ and $\theta_{1}$. The classifier output function can be converted into an approximation of the likelihood ratio between $\theta_{0}$ and $\theta_{1}$! This This manifestation of the Neyman-Pearson lemma in a machine learning setting is often called the likelihood ratio trick.

These three surrogate-based approaches are all amortized: after an upfront simulation and training phase, the surrogates can be evaluated efficiently for arbitrary data and parameter points. They require an upfront specification of the parameters of interest, the network then implicitly marginalizes over all other (latent) variables in the simulator. All three classes of algorithms can employ active learning elements such as an iteratively updated proposal distribution to guide the simulator parameters $\theta$ towards relevant parameter region, improving sample efficiency. Using neural networks eliminates the requirement of low-dimensional summary statistics, leaving it to the employed model to learn the structures in high-dimensional data and potentially improving quality of inference.

Despite these fundamental similarities, there are some differences between emulating the likelihood, the likelihood ratio, and the posterior. Learning the posterior directly provides the main target quantity in Bayesian inference, but induces a prior dependence at every stage of the inference method. Learning the likelihood or the likelihood ratio enables frequentist inference or model comparisons, though for Bayesian inference an additional MCMC or VI step is necessary to generate samples from the posterior. The prior independence of likelihood or likelihood ratio estimators also leads to extra flexibility to change the prior during inference. An advantage of training a generative model to approximate the likelihood or posterior over learning the likelihood ratio function is the added functionality of being able to sample from the surrogate model. On the other hand, learning the likelihood or posterior is an unsupervised learning problem, whereas estimating the likelihood ratio through a classifier is an example of supervised learning and often a simpler task. Since for the higher-level inference goal the likelihood and the likelihood ratio can be used interchangeably, learning a surrogate for the likelihood ratio function may often be more efficient.

Another strategy that allows us to leverage supervised learning is based on extracting additional quantities from the simulator that characterize the likelihood of the latent process (e. g. II and IV from the list in Sec. 2.2.3). This additional information can be used to augment the training data for surrogate models. The resulting supervised learning task can often be solved more efficiently, ultimately improving the sample efficiency in the inference task (52, 78, 17, 53).

Surrogate-based approaches benefit from imposing suitable inductive bias for a given problem. It is widely acknowledged that the network architecture of a neural surrogate should be chosen according to the data structure (e. g. images, sequences, or graphs). Another, potentially more consequential, way of imposing inductive bias is to have the surrogate model reflect the causal structure of the simulator. Manually identifying the relevant structures and designing appropriate surrogate architectures is very domain-specific, though has been shown to improve the performance on some problems (79, 80, 81). Recently attempts are being made to automate the process of creating surrogates that mimic the simulation (82). Looking further ahead, one would like to learn surrogates that reflect the causal structure of a coarse grained system. If this is possible, it would allow the surrogate to model only the relevant degrees of freedom for the phenomena that emerge from the underlying mechanistic model.

3 Preprocessing and postprocessing

There are a number of additional steps that can surround these core inference algorithms, either in the form of preprocessing steps that precede the main inference stage, or as an afterburner following the main inference step.

One preprocessing step is to learn powerful summaries $y(x)$. Because of the curse of dimensionality, both ABC and classical density estimation–based inference methods require a compression of the data into low-dimensional summary statistics. They are usually prescribed, i. e. hand-chosen by domain scientists based on their intuition and knowledge of the problem at hand, but the resulting summaries will generally lose some information compared to the original data. A minimally invasive extension of these algorithms is to first learn summary statistics that have certain optimality properties, before running a standard inference algorithm such as ABC. We sketch this approach in Fig. 3b for ABC, but it applies equally to inference based on density estimation.

The score $t(x|\theta)\equiv\nabla_{\theta}p(x|\theta)$, the gradient of the log (marginal) likelihood with respect to the parameters of interest, defines such a vector of optimal summary statistics: in a neighborhood of $\theta$, the score components are sufficient statistics, and they can be used for inference without loss of information. Just like the likelihood, the score itself is generally intractable, but it can be estimated based on quantity V and an exponential family approximation (83, 84). If quantity II is available, augmented data extracted from the simulator can instead be used to train a neural network to estimate the score (52) without requiring such an approximation. Learned summaries can also be made robust with respect to nuisance parameters (85, 86).

Even if it is not necessary to reduce the data to low-dimensional summary statistics, in some fields the measured raw or “low-level” data can be very high-dimensional. It is then common practice to compress them to a more managable set of “high-level” features of moderate dimensionality and to use these compressed data as input to the inference workflow.

Inference compilation (58) is a preprocessing step for probabilistic programming algorithms, shown in Fig. 3d. Initial runs of the simulator are used to train a neural network used for sequential importance sampling of both the parameters $\theta$ and the latent variables $z$.

After the completion of the core inference workflow, an important question is whether the results are reliable: can the outcome be trusted in the presence of imperfections such as limited sample size, insufficient network capacity, or inefficient optimization?

One solution is to calibrate the inference results. Using the ability of the simulator to generate data for any parameter point, we can use a parametric bootstrap approach to calculate the distribution of any quantity involved in the inference workflow. These distributions can be used to calibrate the inference procedure to provide confidence sets and posteriors with proper coverage and credibility (73, 17). While possible in principle, such procedures may require a large number of simulations.

Other diagnostic tools that can be applied at the end of the inference stage involve training classifiers to distinguish data from the surrogate model and the true simulator (73), checking known expectation values of estimators of the likelihood, likelihood ratio, or score (17); varying reference parameters that should leave the inference result invariant (73); ensemble methods; and comparing distributions of network output against known asymptotic properties (87, 88, 89). Passing these sanity checks does not guarantee that an estimator is correct, but failing them is an indication of problems. Some of these methods may be used for uncertainty estimation, though the statistical interpretation of such error bars is not always obvious.

None of these diagnostics address the issues encountered if the model is misspecified and the simulator is not an accurate description of the system being studied. Model misspecification is a problem that plagues inference with both prescribed and implicit models equally. Usually this is addressed by expanding the model to have more flexibility and introducing additional nuisance parameters.

4 Recommendations

The considerations needed to choose which of the approaches described above is best for a given problem will include the inference goals, the dimensionality of model parameters, the latent variables, and the data; whether good summary statistics are available; the internal structure of the simulator; the computational cost of the simulator; the level of control over how the simulator is run; and on whether the simulator is a black box or whether any of the quantities discussed in Sec. 2.2.3 can be extracted from it. Nevertheless, we believe that the existing body of research lets us provide a few general guidelines.

First, if any of the quantities discussed in Sec. 2 are available, they should be leveraged. Powerful algorithms are available for the case with differentiable simulators (51), for simulators for which the joint likelihood of data and latent variables is accessible (52), and for simulators explicitly written as a probabilistic model in a probabilistic programming framework (90). Probabilistic programming is also the most versatile approach when the goal is not only inference on the parameters $\theta$, but also the latent variables $z$.

If powerful low-dimensional summary statistics are established, traditional techniques can still offer a reasonable performance. In most cases, however, we recommend trying methods based on training a neural network surrogate for the likelihood (41, 69) or the likelihood ratio (73, 57, 76). If generating synthetic data from the surrogate is not important, learning the likelihood ratio rather than the likelihood allows us to leverage powerful supervised learning methods.

Finally, active learning techniques can improve the sample efficiency for all inference techniques. There is a tradeoff between active learning, which tailors the efficiency to a particular observed data set, and amortization, which benefits from surrogates that are agnostic about the observed data. A good compromise here will depend on the number of observations and the sharpness of the posterior compared to the prior.

Discussion

Until recently, scientists confronted with inverse problems and a complex simulator as a forward model had little recourse than to choose ABC or methods based on classical density estimation techniques. While these approaches have served some domains of science quite well, they have relied heavily on experts providing powerful summary statistics. As a result, the techniques are labor intensive and do not lend themselves well to high-dimensional data where powerful summary statistics are not obvious. While not explicit, there is a frontier where these traditional methods are no longer useful, beyond which scientists must resort to other heuristics not framed as a statistical statement tied to the underlying mechanistic model.

The term likelihood-free inference has served as a point of convergence for what were previously disparate communities, and a new lingua franca has emerged. This has catalyzed significant cross-pollination and led to a renaissance in simulation-based inference. The advent of powerful machine learning methods is enabling practitioners to work directly with high dimensional data and to reduce the reliance on expert-crafted summary statistics. New programming paradigms such as probabilistic programming and differentiable programming provide new capabilities that enable entirely new approaches to simulation-based inference. Finally, taking a more systems-level view of simulation-based inference that brings together the statistical and computational considerations has taken root. Here active learning is leading the way, but we expect more advances like this as simulation-based inference matures.

The rapidly advancing frontier means that several domains of science should expect either a significant improvement in inference quality or the transition from heuristic approaches to those grounded in statistical terms tied to the underlying mechanistic model. It is not unreasonable to expect that this transition may have a profound impact on science.

KC and JB are supported by the National Science Foundation under the awards ACI-1450310, OAC-1836650, and OAC-1841471 and by the Moore-Sloan data science environment at NYU. GL is recipient of the ULiège-NRB Chair on Big Data and is thankful for the support of NRB. \showacknow

References