Implicit Likelihood Inference (ILI)

• Implicit Likelihood Inference is a family of simulation-based methods that enable inference without tractable density evaluation.
• These methods leverage neural density estimation, variational techniques, and test inversion to approximate posteriors and likelihoods from simulated data.
• ILI has broad applications in astrophysics, ecology, and particle physics, highlighting challenges in simulation fidelity, computational front-loading, and uncertainty calibration.

Implicit Likelihood Inference (ILI) denotes a family of inference procedures for settings in which a forward model can generate samples but the likelihood cannot be evaluated in closed form. In the simulator-based literature, ILI is also called likelihood-free inference or simulation-based inference, and methods learn either the posterior, the likelihood, or a likelihood ratio from simulated $(\theta,x)$ pairs. In adjacent work on implicit probabilistic models, the same difficulty appears when a model is defined only by a sampling procedure, motivating estimators that avoid explicit density evaluation while retaining likelihood-based semantics under stated conditions (Ho et al., 2024, Zhao et al., 2022, Li et al., 2018).

1. Conceptual domain and model classes

An implicit probabilistic model is defined only by a procedure for sampling, rather than by an explicit density function. A canonical construction is

$z \sim N(0,I), \qquad x = T_\theta(z),$

where $T_\theta$ is a highly expressive neural-network mapping. The induced density $p_\theta(x)$ exists in principle, but in practice it involves an intractable high-dimensional integral and differentiation, so standard maximum-likelihood inference $\sum_i \log p_\theta(x_i)$ is ill-defined or impossible to evaluate (Li et al., 2018).

The same idea appears in hierarchical form in Hierarchical Implicit Models (HIMs). A HIM extends the usual Bayesian hierarchical factorization

$$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$

but allows the likelihood $p(x_n\mid z_n,\beta)$ to be defined only implicitly by a simulator: $$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$ One can therefore sample from $p(x_n\mid z_n,\beta)$ but cannot evaluate its density in closed form (Tran et al., 2017).

In the simulation-based inference literature, the observational setup is typically written as data generated from an unknown simulator $F_\theta$, with intractable likelihood $z \sim N(0,I), \qquad x = T_\theta(z),$0 but available forward simulation $z \sim N(0,I), \qquad x = T_\theta(z),$1. In this sense, ILI is not restricted to one architecture or one inferential paradigm; it covers posterior estimation, likelihood estimation, ratio estimation, variational approximations, and test inversion, provided the core regime is sample access without tractable density evaluation (Dalmasso et al., 2021, Ho et al., 2024).

2. Statistical objectives and inferential targets

The most common Bayesian target is the posterior

$z \sim N(0,I), \qquad x = T_\theta(z),$2

where the likelihood is unavailable analytically but the simulator can generate $z \sim N(0,I), \qquad x = T_\theta(z),$3. DELFI and related neural-density-estimation methods therefore train a conditional density estimator directly on simulated pairs and use the learned model as an approximation to $z \sim N(0,I), \qquad x = T_\theta(z),$4 or $z \sim N(0,I), \qquad x = T_\theta(z),$5 (Zhao et al., 2022, Wang, 19 Dec 2025).

A distinct Bayesian formulation arises in Likelihood-Free Variational Inference (LFVI) for HIMs. There the posterior

$z \sim N(0,I), \qquad x = T_\theta(z),$6

is approximated by an implicit variational family

$z \sim N(0,I), \qquad x = T_\theta(z),$7

and the ELBO is rewritten in terms of a log-density-ratio

$z \sim N(0,I), \qquad x = T_\theta(z),$8

This reformulation removes the need for closed-form densities for both the model likelihood and the local variational factors (Tran et al., 2017).

A non-Bayesian formulation appears in Implicit Maximum Likelihood Estimation (IMLE). Given data points $z \sim N(0,I), \qquad x = T_\theta(z),$9 and $T_\theta$0 model samples $T_\theta$1, IMLE defines

$T_\theta$2

and optimizes

$T_\theta$3

Under mild regularity and identifiability conditions, the minimizer coincides with that of $T_\theta$4, and if at $T_\theta$5 all $T_\theta$6 coincide then the simple objective $T_\theta$7 has exactly the same minimizer as the likelihood-based one (Li et al., 2018).

ILI also has a frequentist branch. Likelihood-Free Frequentist Inference (LF$T_\theta$8I), described as also called Implicit Likelihood Inference, seeks a confidence set

$T_\theta$9

such that for every $p_\theta(x)$0,

$p_\theta(x)$1

This reframes the problem from posterior approximation to finite-sample or near finite-sample validity of test inversion in settings with intractable likelihoods (Dalmasso et al., 2021).

3. Core algorithmic mechanisms

One major ILI family is neural density estimation. LtU-ILI organizes this space into Neural Posterior Estimation (NPE), Neural Likelihood Estimation (NLE), and Neural Ratio Estimation (NRE). NPE trains $p_\theta(x)$2 directly; NLE trains $p_\theta(x)$3 and combines it with the prior at inference time; NRE trains a classifier $p_\theta(x)$4 to distinguish joint samples from product-of-marginals and uses the standard density-ratio trick

$p_\theta(x)$5

The same ecosystem includes Mixture Density Networks, Masked Autoregressive Flow, Neural Spline Flow, and Rational Quadratic Flows (Ho et al., 2024).

A second mechanism is explicit density-ratio estimation inside a variational objective. In LFVI, a binary classifier is trained to distinguish samples from the model joint $p_\theta(x)$6 and the variational joint $p_\theta(x)$7, with logistic objective

$p_\theta(x)$8

At optimum, $p_\theta(x)$9, allowing Monte Carlo gradients of the surrogate ELBO (Tran et al., 2017).

A third mechanism is odds-based test-statistic construction. LF$\sum_i \log p_\theta(x_i)$0I defines a learned odds function

$\sum_i \log p_\theta(x_i)$1

where $\sum_i \log p_\theta(x_i)$2 labels whether $\sum_i \log p_\theta(x_i)$3 came from $\sum_i \log p_\theta(x_i)$4 or from a reference distribution $\sum_i \log p_\theta(x_i)$5. From this, the framework studies two statistics: ACORE, a max-based likelihood-ratio analogue, and BFF, an average-based statistic that coincides with the Bayes factor under perfect odds estimation (Dalmasso et al., 2021).

A fourth mechanism uses auxiliary information mined from simulators. When a simulator exposes the latent trace $\sum_i \log p_\theta(x_i)$6 and tractable conditional densities along the trace, one can compute the joint score

$\sum_i \log p_\theta(x_i)$7

and the joint likelihood ratio

$\sum_i \log p_\theta(x_i)$8

The resulting losses include Rolr, Rascal, Cascal, and Scandal, which use joint score and joint ratio targets to improve sample efficiency relative to purely sample-based surrogates (Brehmer et al., 2018).

A fifth mechanism concerns summary construction rather than density parameterization. For cosmological fields, Information Maximising Neural Networks (IMNNs) learn summaries $\sum_i \log p_\theta(x_i)$9 by maximizing the determinant of the Fisher matrix under a Gaussian approximation for the summaries. Implicit inference is then carried out on these summaries using ABC or DELFI. This setup is explicitly presented as implicit likelihood inference with maximally informative nonlinear summaries (Makinen et al., 2021).

4. Theoretical guarantees, calibration, and diagnostics

The theoretical claims surrounding ILI are heterogeneous because the inferential goals differ. For IMLE, the central result is equivalence to maximizing likelihood under stated assumptions, including differentiability in $$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$0, continuity in $$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$1, a shift-invariance property, an envelope property guaranteeing a parameter $$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$2 that maximizes local CDFs and densities, uniqueness of the MLE $$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$3, and non-degeneracy of certain gradients. The paper’s Theorem 1 gives an argmin/argmax equivalence between a weighted nearest-sample objective and the log-likelihood objective (Li et al., 2018).

For LFVI, the guarantees are formulated in terms of ratio-estimator consistency, ELBO maximization, and stochastic optimization. If the ratio network is a universal function class and $$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$4 is perfectly minimized, then $$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$5 and the surrogate ELBO recovers the true KL objective. Under standard Robbins–Monro conditions and unbiased mini-batch gradients, the variational parameters converge to a stationary point of the ELBO (Tran et al., 2017).

For LF$$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$6I, the emphasis is coverage. The framework learns a critical-value function $$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$7 by conditional quantile regression and then inverts tests through

$$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$8

Under mild regularity on the quantile regressors and smoothness of the null distribution in $$p(x,z,\beta)=p(\beta)\prod_{n=1}^N p(z_n)\,p(x_n\mid z_n,\beta)$$9, $p(x_n\mid z_n,\beta)$0 converges uniformly to the true $p(x_n\mid z_n,\beta)$1-quantile, the tests have Type-I error $p(x_n\mid z_n,\beta)$2, and the p-value-based procedure also has correct size. LF$p(x_n\mid z_n,\beta)$3I further adds an empirical-coverage diagnostic that regresses the indicator of coverage on $p(x_n\mid z_n,\beta)$4 in order to estimate where coverage falls below or above nominal level (Dalmasso et al., 2021).

Amortized astrophysical pipelines operationalize calibration with simulation-based diagnostics. LtU-ILI includes Probability Integral Transform histograms, P–P plots, the TARP test, test log-likelihood, and classifier two-sample tests, explicitly to assess whether learned posteriors are informative and well calibrated (Ho et al., 2024). Specific applications also report calibration checks. The 21 cm DELFI study states that all posteriors are statistically validated by PIT, copPIT, and HPD tests (Zhao et al., 2022). The set-based galaxy-cluster study evaluates highest-posterior-density credible intervals and Tests of Accuracy with Random Points, reporting coverage within a few percent of perfect calibration across $p(x_n\mid z_n,\beta)$5 (Wang et al., 27 Jul 2025). The neutrino-mass-hierarchy application reports coverage and rank tests on 500 held-out simulations, with rank statistics and P–P plots consistent with uniformity to within statistical errors (Wang, 19 Dec 2025).

A recurrent theme is that valid uncertainty quantification is not automatic. LF$p(x_n\mid z_n,\beta)$6I explicitly states that popular likelihood-free methods do not necessarily lead to valid scientific inference because they do not guarantee confidence sets with nominal coverage in general settings, and LtU-ILI notes that single-network neural density estimators tend to under-cover. This has made posterior validation a central component of contemporary ILI practice (Dalmasso et al., 2021, Ho et al., 2024).

5. Application domains and representative results

ILI has been applied across ecology, text modeling, hidden-state models, particle physics, and a wide range of astrophysical and cosmological problems. The applications differ in observable type, simulator structure, and inferential target, but they share the same operational premise: simulate first, learn an inferential surrogate second.

Domain	Observable or data object	ILI mechanism
Hierarchical simulators	predator-prey populations, discrete data, text generation	LFVI with implicit variational family
Cosmological fields	Gaussian and log-normal fields	IMNN compression with ABC or DELFI
21 cm cosmology	power spectra, light-cones, light-cubes	DELFI, normalizing flows, 3D ScatterNet
Galaxy clusters	projected galaxy dynamics	Deep Sets with conditional normalizing flow
Cosmological parameter inference	CMB power spectra, BAO distances	SNLE in LtU-ILI
Implicit HMMs	latent states and observations	autoregressive-flow smoothing approximation

In the original LFVI paper, the applications include a large-scale physical simulator for predator-prey populations in ecology, a Bayesian generative adversarial network for discrete data, and a deep implicit model for text generation (Tran et al., 2017). In simulator-rich scientific settings, the resulting workflows are often fully amortized. LtU-ILI presents applications to galaxy cluster masses from X-ray photometry, cosmology from matter power spectra and halo point clouds, gravitational-wave progenitor characterization, dust-attenuation parameters from galaxy colors and luminosities, and semi-analytic galaxy-formation models (Ho et al., 2024).

Several cosmology papers provide quantitative benchmarks. For field-level cosmological inference, IMNN compression recovers $p(x_n\mid z_n,\beta)$7 of the analytic Fisher determinant in the Gaussian $p(x_n\mid z_n,\beta)$8 example within $p(x_n\mid z_n,\beta)$9 training epochs, and in the log-normal case a single IMNN trained far from the target saturates to $$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$0 of the Shannon information in $$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$1 epochs, with one retraining at $$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$2 yielding $$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$3 recovery (Makinen et al., 2021). For 21 cm power-spectrum inference, the pure-signal ILI analysis recovers

$$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$4

while the standard MCMC baseline gives wider and shifted intervals; with realistic HERA- or SKA-level noise, ILI and the Gaussian-likelihood MCMC analysis agree on median locations, but ILI still yields slightly smaller $$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$5 widths (Zhao et al., 2022). For 21 cm light-cones, 3D ScatterNet reports $$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$6 for $$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$7 and $$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$8 for $$\beta\sim p(\beta), \qquad z_n\sim p(z_n), \qquad \epsilon_n\sim s(\epsilon), \qquad x_n=g(\epsilon_n;z_n,\beta).$$9 on pure-signal tests, outperforming both power-spectrum summaries and a fine-tuned 3D CNN (Zhao et al., 2023).

In galaxy-cluster inference, the set-based ILI framework learns the posterior over cluster mass from the full projected phase-space of member galaxies rather than a single summary statistic. The baseline $p(x_n\mid z_n,\beta)$0–$p(x_n\mid z_n,\beta)$1 relation has $p(x_n\mid z_n,\beta)$2 dex, whereas the ILI posterior means satisfy $p(x_n\mid z_n,\beta)$3 dex, a factor-of-two improvement; the mean residual is within $p(x_n\mid z_n,\beta)$4 dex across $p(x_n\mid z_n,\beta)$5, and HPD coverage errors are below $p(x_n\mid z_n,\beta)$6 for $p(x_n\mid z_n,\beta)$7 (Wang et al., 27 Jul 2025).

ILI has also been used for calibration and latent-state reconstruction. In the CAMELS-based calibration study, the ILI from the emulated SFRD recovers target observables with relative error $p(x_n\mid z_n,\beta)$8, and the SMF-based inference does so with $p(x_n\mid z_n,\beta)$9; the work reports that stellar-mass functions can break degeneracies present in the cosmic star-formation rate density (Jo et al., 2022). For implicit HMMs, the IDE method learns the posterior over hidden states directly with autoregressive flows and reports hidden-state RMSE comparable to a much more computationally expensive SMC algorithm (Ghosh et al., 2024). In the neutrino-mass-hierarchy application, a two-round SNLE procedure using $F_\theta$0 simulations per round gives

$F_\theta$1

slightly preferring $F_\theta$2, hence the normal hierarchy (Wang, 19 Dec 2025).

6. Limitations, recurrent difficulties, and open directions

A persistent limitation is dependence on simulator fidelity. The LIMFAST multiline-intensity-mapping study states that all inference is conditional on the fidelity of the forward model and that systematic misspecification will generically bias posteriors. The same work also notes that one must densely sample $F_\theta$3 over the entire prior region of interest, since extrapolation beyond the training domain is unsafe, and that foregrounds and interloper-line contamination are not yet included in the proof-of-concept (Sun et al., 8 Sep 2025). LtU-ILI generalizes this concern as model misspecification: if the simulator or training catalog is not fully representative of real data, the posterior can be unreliable (Ho et al., 2024).

A second difficulty is computational front-loading. Although amortized inference is often fast once trained, the training database can be large. The 21 cm DELFI study uses 18,000 simulations for the pure-signal case and 90,000 summaries for the noise-plus-foreground case, with training time of $F_\theta$4 hours on an HPC node; inference is then reduced to approximately $F_\theta$5 minutes on one CPU core (Zhao et al., 2022). LF$F_\theta$6I addresses a related issue by separating simulation for test-statistic estimation, calibration, and diagnostics into amortized modules, so that no per-$F_\theta$7 retraining is required (Dalmasso et al., 2021).

A third limitation is access to simulator internals. Methods based on joint score and joint likelihood ratio require the simulator to expose tractable conditional densities along the latent trace or to be rewritten in an autodiff-friendly framework. The simulator-mining paper notes that not all legacy codes expose these quantities without intrusive rewriting (Brehmer et al., 2018). This distinguishes such methods from black-box density-estimation approaches, which require only simulation.

A fourth issue is the difference between informative posteriors and valid uncertainty statements. A common misconception is that likelihood-free methods automatically deliver calibrated credible or confidence regions. Several papers explicitly reject this: LF$F_\theta$8I argues that popular likelihood-free methods do not guarantee nominal coverage in general settings, and LtU-ILI warns that single-network neural density estimators tend to under-cover, recommending coverage and simulation-based calibration tests as standard practice (Dalmasso et al., 2021, Ho et al., 2024).

The literature also identifies domain-specific remedies and extensions. IMLE emphasizes that its pure minimization objective avoids adversarial-game instability, vanishing gradients, and mode collapse under its nearest-sample formulation (Li et al., 2018). The simulator-mining work points to automatic probabilistic-programming wrappers, active learning, and integration with recent normalizing-flow advances (Brehmer et al., 2018). LtU-ILI highlights deep ensembles, dropout, stochastic weight averaging, and future incorporation of score-based generative models, diffusion, flow matching, and conditional denoising flows for higher-dimensional parameter spaces (Ho et al., 2024). This suggests that ILI is evolving less as a single method than as an increasingly modular inferential infrastructure for implicit models and simulator-defined scientific theories.