Benchmarking Simulation-Based Inference

Motivation

• Existing SBI comparisons are difficult because studies use inconsistent tasks and metrics, commonly used metrics may be inappropriate or biased, and reimplementation is wasteful.These problems hinder identifying promising approaches and choosing algorithms, especially across closely related methods developed in different disciplines.Source: §1 · pp. 1–2
• An extensible benchmark can support community-driven algorithm development, help practitioners choose methods, and improve transparent and reproducible computational research.The framework is intended to address future challenges including high-dimensional data, active simulation selection, and gray-box settings.Source: §1 · p. 2

Benchmark contribution

• The benchmark provides initial SBI algorithms, ten tasks, analyzed performance metrics, reference posteriors, open code, and an interactive results website.Tasks with evaluable likelihoods enable ground-truth posterior computation for rapid algorithm evaluation.Source: §1 · p. 2

Initial findings

• Metric choice is critical, algorithms retain substantial room for improvement, and sequential estimation generally improves sample efficiency.These findings are among the benchmark’s initial insights.Source: §1 · p. 2
• For small and moderate simulation budgets, neural-network-based approaches outperform classical ABC algorithms, but no algorithm dominates across tasks.Algorithm rankings are task-dependent, motivating better guidance or automated procedures for algorithm selection.Source: §1 · p. 2

2 Benchmark

• The benchmark tests approximate posteriors q(θ|xo) against reference posteriors p(θ|xo) using defined performance metrics.Each task specifies a prior, simulator, observation, and inference output, which is evaluated against the reference posterior.Source: §2 · p. 2

2.1 Algorithms

• The algorithm suite covers classical ABC and neural approaches to likelihood, posterior, and density-ratio estimation, with sequential and non-sequential variants.Sequential methods adaptively acquire simulations to improve sample efficiency, but introduce additional algorithmic and hyperparameter choices.Source: §2.1 · pp. 2–3, Figure 1 · p. 3

2.2 Performance metrics

• Classifier two-sample-test accuracy is preferred for comparison because commonly used alternatives can be unreliable, including probability-of-true-parameters, PPCs, and MMD.Probability-of-true-parameters requires averaging over many prior-sampled observations, PPCs can accept poor posterior estimates, and MMD is sensitive to kernel hyperparameters, especially for multimodal posteriors.Source: §2.2 · p. 4

2.3 Tasks

• The tasks combine eight statistical problems and two applied-domain problems with varied parameter and data dimensionalities, including multimodal, truncated, distractor-rich, epidemiological, and ecological settings.Reference posterior samples enable two-sample testing; examples include SLCP, Two Moons, SIR, and Lotka-Volterra.Source: §2.3 · pp. 4–5

2.4 Experimental Setup

• Each task uses 10 prior-sampled parameter-observation pairs, 10k reference-posterior samples per observation, and algorithm budgets from 1k to 100k simulations.Performance is reported across simulation budgets because simulation time is assumed to dominate total computation cost.Source: §2.4 · p. 5

2.5 Software

• The released framework provides extensible tasks, reference posteriors, metrics, plotting, and tooling, with tasks implemented as probabilistic programs in Pyro.It supports external algorithm toolboxes and ODE-based tasks through an interface between DifferentialEquations.jl and PyTorch.Source: §2.5 · p. 5, Code · p. 5

Reproducibility.

• A public website enables interactive exploration of all benchmark results.The code and reproducibility instructions are publicly available, and the web interface accompanies the released framework.Source: Reproducibility · p. 5, §2.5 · p. 5

Metric analysis

• Performance conclusions depended critically on the metric: C2ST reflected budget and sequential-estimation effects that MMD and MEDDIST did not, while KSD encountered numerical problems.MMD was sensitive to hyperparameters and struggled with multimodal posteriors under the median heuristic; median distance could be gamed by point estimates.Source: §3 · p. 6
• The authors used C2ST for reported performance because ground-truth log-probability correlated weakly with C2ST when averaged over few observations.The reported correlations were 0.3 on Two Moons and 0.6 on SLCP.Source: §3 · p. 6

Main results

• C2ST showed that several benchmark tasks remained unsolved within the specified simulation budgets, demonstrating substantial room for more powerful algorithms.C2ST uses an interpretable scale from 1 to 0.5; SLCP and Lotka-Volterra were examples where no algorithm solved the task.Source: §3 · p. 6
• Sequential algorithms generally outperformed non-sequential methods, especially on harder tasks, but their gains diminished as simulation budgets increased.The difference was small for linear Gaussian tasks and pronounced for most other tasks.Source: §3 · p. 6, Figure 3 · p. 7
• Neural density- or ratio-estimation methods generally outperformed REJ-ABC and SMC-ABC, while classical rejection methods remained competitive for low-dimensional problems with cheap simulators.Classical methods require orders-of-magnitude less computation because they avoid network training.Source: §3 · p. 6, Figure 3 · p. 7
• No single sequential method—SNLE, SNRE, or SNPE—was consistently preferable because their distinct strengths interacted with problem structure.The benchmark therefore supports practical, problem-dependent algorithm selection rather than a one-size-fits-all recommendation.Source: §3 · pp. 7–8

Algorithm diagnosis

• Reference posteriors helped diagnose MCMC chains getting stuck in single modes, and revised sampling strategies and unbounded parameter transformations improved performance and speed.These issues particularly affected (S)NLE and (S)NRE, including on Lotka-Volterra.Source: §3 · pp. 7, 9

Practitioners’ advice

• Method choice should reflect simulator access, dimensionality, simulation cost, and reuse: classical ABC suits low-dimensional settings with abundant cheap simulations, whereas model-based or sequential methods suit harder or costly problems.NPE is preferable for repeated inference on different datasets, while NLE and NRE suit many i.i.d. observations; sufficient summaries are required to preserve the full posterior.Source: §3 · p. 8

Algorithm coverage

• The benchmark omits several SBI families, including active-learning or Bayesian-Optimization sequential methods, gray-box methods, Gaussian-Process alternatives, and many other algorithms.Its deliberately small initial selection enabled careful hyperparameter investigation and emphasized less sophisticated sequential acquisition schemes in black-box settings.Source: §4 · p. 9

Task coverage

• The tasks exclude high-dimensional spatial structure such as images, which would require algorithms to learn summary statistics while exploring data structure.Such methods remain an active research area.Source: §4 · p. 9

Evaluation bias

• Results may reflect the authors’ areas of expertise despite extensive tuning-choice investigations and implementation comparisons.Source: §4 · p. 9

Computational cost

• Simulation-call performance can obscure substantial runtime differences: (S)ABC was much faster than network-training approaches, while sequential neural algorithms had the longest runtimes.Source: §4 · p. 9

Generalisability

• Because the benchmark focuses on problems with computable reference posteriors, its insights may not generalize directly to real-world simulators, although it remains a necessary stepping stone toward practical selection strategies.Direct comparison on real-world simulators is difficult because appropriate metrics are unavailable.Source: §4 · p. 9

Practical selection

• Algorithm choice depends on difficult-to-quantify factors including available problem information, inference goals, and simulator speed.The paper provides practical considerations and recommendations in Box 1.Source: §4 · p. 9

Benchmark risks

• Benchmarks can encourage overfitting and implementation-focused illusions of progress, while novel ideas may initially appear uncompetitive; therefore, hyperparameters were selected for performance across tasks rather than task-specific tuning.Because simulators are available, excessive task-specific tuning could use samples or likelihoods that would not transfer to practice.Source: §4 · p. 9

• The initial benchmark provides a basis for quantitatively evaluating, comparing, and improving simulation-based inference algorithms.Used sensibly, it is intended to clarify and expedite progress in the field.Source: §5 · p. 9
• The benchmark’s results and recommendations are intended to inform algorithm development and help practitioners select appropriate algorithms.Source: §5 · p. 9
• Precomputed reference posteriors, an interactive website, and open-source code support participation and rapid comparison of algorithms.The stated resources include sbi-benchmark.github.io and github.com/sbi-benchmark/sbibm.Source: §5 · p. 9

A.1 Rejection Approximate Bayesian Computation (REJ-ABC)

• REJ-ABC accepts prior-sampled parameters whose simulated data lie within a distance threshold of the observation, but accuracy requires balancing tolerance and simulation budget.As data dimensionality increases linearly, maintaining accuracy requires an exponential increase in simulation budget.Source: §A.1 · p. 17

A.2 Sequential Monte Carlo Approximate Bayesian Computation (SMC-ABC)

• SMC-ABC constructs successive perturbed particle populations approaching the posterior, improving sampling efficiency over REJ-ABC while requiring careful tuning and retaining dimensionality problems.Its tunable components include the distance measure, perturbation kernel, and acceptance-threshold schedule.Source: §A.2 · p. 18

A.3 Neural Likelihood Estimation (NLE)

• NLE trains a conditional neural density estimator to approximate p(x|θ), then uses the surrogate likelihood with MCMC to obtain approximate posterior samples.With enough simulations and sufficient flexibility, the estimator can approximate the likelihood within the prior support.Source: §A.3 · p. 19

A.4 Sequential Neural Likelihood Estimation (SNLE)

• SNLE extends NLE across rounds by sampling parameters from the previous posterior estimate, concentrating simulations near the observation to improve sample efficiency.This proposal-focused training can be more efficient than drawing all training parameters from the prior.Source: §A.4 · p. 20

A.5 Neural Posterior Estimation (NPE)

• NPE directly learns the posterior qψ(θ|x) with conditional density estimation, recovering the ground-truth posterior in the infinite-data, sufficiently flexible-model limit.The estimator is trained by minimizing the negative log probability of parameters conditioned on simulated data.Source: §A.5 · p. 21

A.6 Sequential Neural Posterior Estimation (SNPE)

• SNPE sequentially targets posterior regions to increase sample efficiency, but requires correcting for proposal distributions that differ from the prior.The benchmark uses the atomic-proposal approach of Greenberg et al. (2019), referred to as SNPE.Source: §A.6 · p. 22

A.7 Neural Ratio Estimation (NRE)

• NRE trains classifiers to estimate likelihood ratios, enabling posterior evaluation or MCMC sampling without directly evaluating the intractable likelihood.The benchmark uses the formulation of Durkan et al. (2020) with a residual-network classifier and MCMC sampling.Source: §A.7 · p. 23

A.8 Sequential Neural Ratio Estimation (SNRE)

• SNRE sequentially refines a ratio estimator using posterior-focused proposals, increasing sample efficiency but requiring new classifiers for different observations.The proposal for each round is the previous posterior estimate, and posterior samples are obtained through MCMC.Source: §A.8 · p. 24

A.9 Random Forest Approximate Bayesian Computation (RF-ABC)

• RF-ABC uses random-forest regression of summary statistics on parameters to assign weights to simulated parameters for approximate posterior sampling.The benchmark used 500 trees, a minimum node size of 5, and 10k posterior samples generated from the resulting weights.Source: §A.9 · p. 25
• RF-ABC is applicable only to one-dimensional parameters or multidimensional spaces where the posterior factorizes, so it is reported only for compatible benchmark tasks in Supplementary Figure 2.The compatible tasks were Gaussian Linear, Gaussian Linear Uniform, and Gaussian Mixture.Source: §A.9 · p. 25

A.10 Synthetic Likelihood (SL)

• Synthetic Likelihood estimates a Gaussian likelihood from simulated summary statistics and uses it within MCMC, capturing noise that distance-based ABC may treat as problematic.The benchmark implementation used Slice Sampling and 100 simulations to estimate the Gaussian likelihood at each sampling step.Source: §A.10 · p. 26
• SL required approximately 10^8 to 10^9 simulations per run to generate 10k posterior samples, making direct comparison with other algorithms problematic and motivating its supplementary reporting.The implementation used covariance regularization for selected tasks and reported SL in Supplementary Figure 3.Source: §A.10 · p. 26

B.1 Reference posteriors

• Reference posteriors used analytic, custom, or likelihood-based sampling schemes across tasks, with repeated runs producing indistinguishable C2ST results.For multimodal tasks, the authors avoided direct MCMC because standard methods struggled with multimodal posteriors and potential chain-correlation bias.Source: §B.1 · p. 27

B.2 Code

• sbibm implements benchmark tasks, reference posteriors, metrics, and scalable tooling in an open, modular framework that integrates with SBI toolboxes.The framework was evaluated with pyABC, pyabcranger, and sbi, while remaining toolbox-agnostic.Source: §B.2 · p. 27

F Figures

• Supplementary ABC analyses found that linear regression adjustment sometimes improved results but decreased performance on the bimodal Two Moons task, while SASS matched abcpy performance.These results used post-hoc LRA and a third-order polynomial feature expansion for SASS.Source: Figure 1 · p. 28
• RF-ABC performance is limited on tasks with non-factorized posteriors because it predicts parameters independently, despite comparisons against REJ-ABC and SMC-ABC using C2ST.The factorized assumption is appropriate only for Gaussian Linear, Gaussian Linear Uniform, and Gaussian Mixture.Source: Figure 2 · p. 29
• SL required approximately 10^8 to 10^9 simulations and was omitted from the main paper because its budget greatly exceeded those of other algorithms.Stable covariance estimation was also unavailable for SLCP Distractors and Lotka-Volterra, so those tasks were omitted from the SL analysis.Source: Figure 3 · p. 30
• On Two Moons, MMD’s median heuristic ranked REJ-ABC slightly ahead of SNLE despite C2ST and visual evidence favoring SNLE, whereas a manually chosen length scale reversed the comparison.The reported median-heuristic MMD values were 0.00729 for REJ-ABC versus 0.00772 for NLE; with the manually chosen length scale they were 0.00738 versus 0.00035, respectively.Source: Figure 4 · p. 31

H Hyperparameter Choices

• Hyperparameters were selected across multiple or all tasks rather than per task, prioritizing settings that generalize and avoid task-specific overfitting.More than 10 000 runs explored previously unreported choices and revealed substantial improvements, enabling systematic design of more robust SBI algorithms.Source: H Hyperparameter Choices · p. 33

H.1 REJ-ABC

• REJ-ABC used an l2 distance metric, task-specific acceptance tolerances, quantile-based selection, and resampling of 10k retained samples for posterior comparison.A KDE alternative was evaluated to test whether resampling impaired performance.Source: §H.1 · p. 33
• REJ-ABC performed best using the 0.1, 0.01, 0.001 quantile schedule, retaining the top 100 samples and fitting a KDE.The schedule performed best across tasks and simulation budgets, and KDE particularly improved Gaussian-task performance.Source: §H.1 · p. 33, Figure 6 · p. 34

H.2 SMC-ABC

• SMC-ABC tuning swept population size, epsilon-selection quantile, and Gaussian-kernel covariance scale while fixing the Gaussian perturbation kernel and l2 distance metric.The implementation was also compared with pyABC, with optional KDE used to generate samples for two-sample metrics.Source: §H.2 · p. 35
• SMC-ABC performed best with population size 100, covariance scale 0.5, epsilon quantile 0.2, implementation A, and KDE, with population size 1000 at 100k simulations.Version A was slightly better on Gaussian tasks, while both implementations were qualitatively similar; KDE especially improved Gaussian-task performance.Source: §H.2 · p. 35, Figure 7 · p. 36, Figure 8 · p. 37

H.3 MCMC for (S)NLE and (S)NRE

• Using 100 MCMC chains with resampling-based initialization was essential for multimodal posteriors, while unbounded transformations and z-scoring prevented sampling failures on some tasks.A single chain often missed modes on Two Moons, and untransformed sampling could enter endless loops on Lotka-Volterra.Source: §H.3 · p. 38
• Vectorizing Slice Sampling enabled parallel likelihood evaluations and reduced SNLE runtime on Gaussian Linear from over 36 hours to under 2 hours for 100k simulations.Batching flow-based density-estimator evaluations enabled the speed-up.Source: §H.3 · p. 38

H.4 Density estimator for (S)NLE

• NSFs generally performed worse than MAFs for (S)NLE, with significant reductions on SIR and Lotka-Volterra, so MAFs were reported in the main paper.The added flexibility of NSFs was not needed for (S)NLE across the evaluated tasks.Source: §H.4 · p. 38, Figure 9 · p. 39

H.5 Density estimator for (S)NPE

• NSF flexibility improved direct posterior estimation on tasks with complex multimodal structure such as Two Moons and SLCP.Artifacts from density transformation visible in Two Moons posteriors disappeared when using NSFs.Source: §H.5 · p. 40
• To the authors’ knowledge, results on (S)NPE with NSFs had not previously been published.Source: §H.5 · p. 40
• NSFs generally outperformed MAFs for (S)NPE across tasks, especially on Gaussian Mixture, Two Moons, and SIR.Performance was measured using C2ST across tasks, with means and 95% confidence intervals over 10 observations.Source: Figure 10 · p. 40, §H.5 · p. 40
• The authors reported NSF-based (S)NPE results in the main paper because NSFs generally performed better across tasks.The comparison used C2ST performance across tasks.Source: Figure 10 · p. 40

H.6 Classifier choice for (S)NRE

• MLP and ResNet classifiers produced similar (S)NRE results on most tasks, but ResNets performed better on Two Moons and SIR.The comparison evaluated C2ST performance across tasks, with means and 95% confidence intervals over 10 observations.Source: §H.6 · p. 41, Figure 11 · p. 41
• ResNets were used for (S)NRE in the main paper because their advantage appeared at low to medium simulation budgets on Two Moons and SIR.Source: §H.6 · p. 41

M.1 Negative log probability of θo (NLTP)

• NLTP needs many parameter-observation pairs for reliable evaluation, but becomes prohibitively expensive for non-amortized algorithms across hundreds of observations.With few pairs, an incorrect posterior can assign higher probability to the true parameter than the true posterior; asymptotically, NLTP relates to posterior KL divergence.Source: §M.1 · p. 42

M.2 Simulation-based calibration (SBC)

• SBC checks whether posterior samples follow the prior under exact inference, but cannot detect uninformative approximations such as the prior itself and is therefore omitted from the benchmark.The benchmark excludes SBC because it requires inference for many prior-generated observations, making hundreds of runs computationally prohibitive for most algorithms.Source: §M.2 · p. 42

M.3 Median distance (MEDDIST)

• MEDDIST measures the median L2 distance between observed data and posterior-predictive samples, using 10,000 samples, but is only a check and may not test posterior structure.It is computed in data space and requires additional potentially expensive simulator calls; the median reduces sensitivity to outliers relative to the mean.Source: §M.3 · p. 43

M.4 Maximum Mean Discrepancy (MMD)

• MMD compares reference and approximate posteriors with Gaussian kernels and 10,000 samples each, but simple translation-invariant kernels can have low power for complex multimodal structure.On SLCP and Two Moons, kernel limitations were reflected in low correlation between MMD and C2ST; learned kernels are suggested as a remedy.Source: §M.4 · p. 43

M.5 Classifier-based tests (C2ST)

• C2ST trains a classifier to distinguish true from estimated-posterior samples, with indistinguishable samples corresponding to chance accuracy of 0.5.The benchmark used a two-layer ReLU network, 10,000 z-scored samples from each distribution, and five-fold cross-validation.Source: §M.5 · p. 43

M.6 Kernelized Stein Discrepancy (KSD)

• KSD is a one-sample goodness-of-fit test comparing algorithm samples with the gradient of the unnormalized true posterior density.The benchmark used Gaussian kernels with median-heuristic length scales and 10,000 samples from each algorithm.Source: §M.6 · p. 43

Runtime setup

• Runtimes were compared on single CPU cores without GPUs, because the shallow, narrow networks used in the benchmark gained little from GPU training and could incur data-transfer overhead.SNPE on SLCP was slightly slower on a GeForce GTX 1080 than on CPU.Source: R Runtimes · p. 44
• Neural-network algorithms ran on single 3.6 GHz AWS C5 CPU cores, while ABC algorithms used single 2.4 GHz cluster CPUs, with cross-hardware ABC differences below 100ms.The benchmark recorded runtimes for every algorithm and task; Figure 12 reports means and 95% confidence intervals.Source: R Runtimes · p. 44, Figure 12 · p. 45

Runtime results

• Sequential neural-network algorithms took longer than their non-sequential counterparts, while SNPE was the slowest sequential method; REJ-ABC and SMC-ABC were short because they require neither neural-network training nor MCMC.Sequential methods involve 10 rounds of neural-network training.Source: R Runtimes · p. 44, Figure 12 · p. 45
• MAF-based runtimes tended to be faster than NSF-based runtimes, with SNPE on SLCP at 100k simulations differing by about 50 minutes on average.The comparison concerns the density-estimator choices within the neural-network algorithms.Source: R Runtimes · p. 44
• Vectorizing MCMC reduced SNLE runtime on Gaussian Linear from more than 36 hours to less than 2 hours.The reported (S)NLE runtimes depended on vectorized MCMC sampling.Source: R Runtimes · p. 44

T.1 Gaussian Linear; T.2 Gaussian Linear Uniform

• The passages describe inference of the mean of a 10-dimensional Gaussian model with fixed covariance, using a conjugate Gaussian prior in one case and a uniform prior in another.Source: §T.1 · p. 46, §T.2 · p. 46

T.3 SLCP

• T.3 uses five parameters and four two-dimensional observations from a Gaussian likelihood with nonlinear mean and variance, producing a complex posterior despite a simple likelihood.Source: §T.3 · p. 46

T.4 SLCP with Distractors

• T.4 extends T.3 by adding uninformative, randomly reordered dimensions to the observations as distractors.Source: §T.4 · p. 46

T.5 Bernoulli GLM; T.6 Bernoulli GLM Raw

• T.5 infers a smoothness-regularized 10-parameter Bernoulli GLM from sufficient statistics, while T.6 uses the same process with raw observations; both have duration T = 100.Source: §T.5 · p. 47, §T.6 · p. 47

T.7 Gaussian Mixture; T.8 Two Moons

• One task infers the common mean of a two-dimensional mixture of two Gaussian distributions with equal 0.5 weights and covariances I and 0.01 ⊙ I; another illustrates global bimodality and local crescent-shaped posterior structure.Source: §T.7 · p. 47, §T.8 · p. 48

T.9 SIR

• T.9 infers SIR contact and recovery rates from ten infectious-count observations, using LogNormal priors and a simulator with population size N = 1000000 and duration T = 160.Source: §T.9 · p. 48

T.10 Lotka-Volterra

• The task infers four species-interaction parameters from 20 summary statistics comprising both population sizes at 10 evenly-spaced time points, with LogNormal noise of 0.1, duration T = 20, and initial conditions (X(0), Y(0)) = (30, 1).Source: §T.10 · p. 48

Companion website

• The companion website supports interactive metric comparisons and posterior-sample inspection for all runs, which the authors found useful for selecting hyperparameters and diagnosing implementation issues.Figure 13 illustrates C2ST comparisons for sequential algorithms on SLCP and an SNLE posterior for x(1) at 100k simulations.Source: W Website · p. 49, Figure 13 · p. 49