Contrastive Neural Ratio Estimation

Likelihood-to-evidence ratio estimation is usually cast as either a binary (NRE-A) or a multiclass (NRE-B) classification task. In contrast to the binary classification framework, the current formulation of the multiclass version has an intrinsic and unknown bias term, making otherwise informative diagnostics unreliable. We propose a multiclass framework free from the bias inherent to NRE-B at optimum, leaving us in the position to run diagnostics that practitioners depend on. It also recovers NRE-A in one corner case and NRE-B in the limiting case. For fair comparison, we benchmark the behavior of all algorithms in both familiar and novel training regimes: when jointly drawn data is unlimited, when data is fixed but prior draws are unlimited, and in the commonplace fixed data and parameters setting. Our investigations reveal that the highest performing models are distant from the competitors (NRE-A, NRE-B) in hyperparameter space. We make a recommendation for hyperparameters distinct from the previous models. We suggest a bound on the mutual information as a performance metric for simulation-based inference methods, without the need for posterior samples, and provide experimental results.


Introduction
: Conceptual, interpolated map from investigated hyperparameters of proposed algorithm NRE-C to a measurement of posterior exactness using the Classifier Two-Sample Test. Best 0.5, worst 1.0. Red dot indicates NRE-A's hyperparameters, γ = 1 and K = 1 [30]. Purple line implies NRE-B [16] with γ = ∞ and K ≥ 1. NRE-C covers the entire plane, generalizing other methods. Best performance occurs with K > 1 and γ ≈ 1, in contrast with the settings of existing algorithms.
We begin with a motivating example: Consider the task of inferring the mass of an exoplanet θ o from the light curve observations x o of a distant star. We design a computer program that maps hypothetical mass θ to a simulated light curve x using relevant physical theory. Our simulator computes x from θ, but the inverse mapping is unspecified and likely intractable. Simulationbased inference (SBI) puts this problem in a probabilistic context [13,64]. Although we cannot analytically evaluate it, we assume that the simulator is sampling from the conditional probability distribution p(x | θ). After specifying a prior p(θ), the inverse amounts to estimating the posterior p(θ | x o ). This problem setting occurs across scientific domains [1,7,10,11,29] where θ generally represents input parameters of the simulator and x the simulated output observation. Our design goal is to produce a surrogate modelp(θ | x) approximating the posterior for any data x while limiting excessive simulation.
Density estimation [5,54,55] can fit the likelihood [2,15,43,56] or posterior [6,21,42,53] directly; however, an appealing alternative for practitioners is estimating a ratio between distributions [12,16,30,34,69]. Specifically, the likelihood-to-evidence ratio p(θ | x) p(x) = p(θ,x) p(θ)p(x) . Unlike the other methods, ratio estimation enables easy aggregation of independent and identically drawn data x. Ratio and posterior estimation can compute bounds on the mutual information and an importance sampling diagnostic. Estimating p(x | θ) p(x) can be formulated as a binary classification task [30], where the classifier σ • f w (θ, x) distinguishes between pairs (θ, x) sampled either from the joint distribution p(θ, x) or the product of its marginals p(θ)p(x). We call it NRE-A. The optimal classifier has Here, σ represents the sigmoid function, • implies function composition, and f w is a neural network with weights w. As a part of an effort to unify different SBI methods and to improve simulation-efficiency, Durkan et al. [16] reformulated the classification task to identify which of K possible θ k was responsible for simulating x. We refer to it as NRE-B. At optimum where an additional bias, c w (x), appears. g w represents another neural network. The c w (x) term nullifies many of the advantages ratio estimation offers. c w (x) can be arbitrarily pathological in x, meaning that the normalizing constant can take on extreme values. This limits the applicability of verification tools like the importance sampling-based diagnostic in Section 2.2.
The c w (x) term also arises in contrastive learning [23,71] with Ma and Collins [45] attempting to estimate it in order to reduce its impact. We will propose a method that discourages this bias instead. Further discussion in Appendix D.
There is a distinction in deep learning-based SBI between amortized and sequential algorithms which produce surrogate models that estimate any posterior p(θ | x) or a specific posterior p(θ | x o ) respectively. Amortized algorithms sample parameters from the prior, while sequential algorithms use an alternative proposal distribution-increasing efficiency at the expense of flexibility. Amortization is usually necessary to compute diagnostics that do not require samples from p(θ | x o ) and amortized estimators are empirically more reliable [31]. Our study therefore focuses on amortized algorithms.
Contribution We design a more general formulation of likelihood-to-evidence ratio estimation as a multiclass problem in which the bias inherent to NRE-B is discouraged by the loss function and it does not appear at optimum. Figure 1 diagrams the interpolated performance as a function of hyperparameters. It shows which settings recover NRE-A and NRE-B, also indicating that highest performance occurs with settings distant from these. Figure 2 shows the relationship of the loss functions. We call our framework NRE-C 1 and expound the details in Section 2.
An existing importance sampling diagnostic [30] tests whether a classifier can distinguish p(x | θ) samples from from samples from p(x) weighted by the estimated ratio. We demonstrate that, when estimating accurate posteriors, our proposed NRE-C passes this diagnostic while NRE-B does not.
Taking inspiration from mutual information estimation [58], we propose applying a variational bound on the mutual information between θ and x in a novel way-as an informative metric measuring a lower bound on the Kullback-Leibler divergence between surrogate posterior estimate p w (θ | x) and p(θ | x), averaged over p(x). Unlike with two-sample testing methods commonly used in machine learning literature [44], our metric samples only from p(θ, x), which is always available in SBI, and does not require samples from the intractable p(θ | x). Our metric is meaningful to scientists working on problems with intractable posteriors. The technique requires estimating the partition function, which can be expensive. We find the metric to be well correlated with results from two-sample tests.
We evaluate NRE-B and NRE-C in a fair comparison in several training regimes in Section 3. We perform a hyperparameter search on three simulators with tractable likelihood by benchmarking the behavior when (a) jointly drawn pairs (θ, x) are unlimited or when jointly drawn pairs (θ, x) are fixed but we (b) can draw from the prior p(θ) without limit or (c) are restricted to the initial pairs. We also perform the SBI benchmark of Lueckmann et al. [44] with our recommended hyperparameters.

Methods
The ratio between probability distributions can be estimated using the "likelihood ratio trick" by training a classifier to distinguish samples [12,19,27,30,50,66,69]. We first summarize the loss functions of NRE-A and NRE-B which approximate the intractable likelihood-to-evidence ratio . We then elaborate on our proposed generalization, NRE-C. Finally, we explain how to recover NRE-A and NRE-B within our framework and comment on the normalization properties.
NRE-A Hermans et al. [30] train a binary classifier to distinguish (θ, x) pairs drawn dependently p(θ, x) from those drawn independently p(θ)p(x). This classifier is parameterized by a neural network f w which approximates log r(x | θ). We seek optimal network weights w * ∈ arg min p(x) given unlimited model flexibility and data. Details can be found in Appendix A.
NRE-B Durkan et al. [16] train a classifier that selects from among K parameters (θ 1 , . . . , θ K ) which could have generated x, in contrast with NRE-A's binary possibilities. One of these parameters θ k is always drawn jointly with x. The classifier is parameterized by a neural network g w which approximates log r(x | θ). Training is done over B samples by finding where θ k ). Given unlimited model flexibility and data NRE-B's ratio estimate converges to g w * (θ, x) = log p(θ | x) p(θ) + c w * (x). Details are in Appendix A.

Contrastive Neural Ratio Estimation
Our proposed algorithm NRE-C trains a classifier to identify which θ among K candidates is responsible for generating a given x, inspired by NRE-B. We added another option that indicates x was drawn independently, inspired by NRE-A. The introduction of the additional class yields a ratio without the specific c w (x) bias at optimum. Define Θ := (θ 1 , ..., θ K ) and conditional probability We set marginal probabilities p(y = k) := p K for all k ≥ 1 and p(y = 0) := p 0 , yielding the relationship p 0 = 1 − Kp K . Let the odds of any pair being drawn dependently to completely independently be γ := Kp K p0 . We now use Bayes' formula to compute the conditional probability We dropped the NRE-C subscript and substituted in γ to replace the p(y) class probabilities. We train a classifier, parameterized by neural network h w (θ, x) with weights w, to approximate (6) by We note that (7) still satisfies K k=0 q w (y = k | Θ, x) = 1, no matter the parameterization.
Optimization We design a loss function that encourages h w (θ, x) = log p(x | θ) p(x) at convergence, and holds at optimum with unlimited flexibility and data. We introduce the cross entropy loss and minimize it towards w * ∈ arg min w (w). We point out that the final term is symmetric up to permutation of Θ, enabling the replacement of the sum by multiplication with K. When γ and K are known, p 0 = 1 1+γ and p K = 1 K γ 1+γ under our constraints. Without loss of generality, we let θ 1 , . . . , θ K ∼ p(θ) and x ∼ p(x | θ K ). An empirical estimate of the loss on B samples is thereforê In the first term, the classifier sees a completely independently drawn sample of x and Θ while θ K is drawn jointly with x in the second term. In both terms, the classifier considers K choices. In practice, we bootstrap both θ K−1 from the same mini-batch and compare them to the same x, similarly to NRE-A and NRE-B. Proof of the above is in Appendix B.
Recovering NRE-A and NRE-B NRE-C is general because specific hyperparameter settings recover NRE-A and NRE-B. To recover NRE-A one should set γ = 1 and K = 1 in (9) yieldinĝ where we dropped the lower index. Recovering NRE-B requires taking the limit γ → ∞ in the loss function. In that case, the first term goes to zero, and second term converges to the softmax function.
is determined by substitution into (9). Both equations are obviously the same as their counterparts.
Estimating a normalized posterior In the limit of infinite data and infinite neural network capacity (width, depth) the optimal classifier trained using NRE-C (with γ ∈ R + ) satisfies the equality: In particular, we have that the following normalizing constant is trivial: This is a result of Lemma 1 in Appendix B. However, practitioners never operate in this setting, rather they use finite sample sizes and neural networks with limited capacity that are optimized locally. The non-optimal function exp(h w (θ, x)) does not have a direct interpretation as a ratio of probability distributions, rather as the function to weigh the prior p(θ) to approximate the unnormalized posterior.
In other words, we find the following approximation for the posterior p(θ | x): where in general the normalizing constant is not trivial, i.e. Z w (x) = 1. As stated above, the NRE-C (and NRE-A) objective encourages Z w (x) to converge to 1. This is in sharp contrast to NRE-B, where even at optimum with an unrestricted function class a non-trivial x-dependent bias term can appear.
There is no restriction on how pathological the NRE-B bias c w (x) can be. Consider a minimizer of (4), the NRE-B loss function, h w * + c w * (x). Adding any function d(x) cancels out in the fraction and is also a minimizer of (4). This freedom complicates any numerical computation of the normalizing constant and renders the importance sampling diagnostic from Section 2.2 generally inapplicable. We report Monte Carlo estimates of Z w (x) on a test problem across hyperparameters in Figure 14.

Measuring performance & ratio estimator diagnostics
SBI is difficult to verify because, for many use cases, the practitioner cannot compare surrogate p w (θ | x) to the intractable ground truth p(θ | x). Incongruous with the practical use case for SBI, much of the literature has focused on measuring the similarity between surrogate and posterior using two-samples tests on tractable problems. For comparison with literature, we first reference a two-sample exactness metric which requires a tractable posterior. We then discuss diagnostics which do not require samples from p(θ | x), commenting on the relevance for each NRE algorithm with empirical results. Further, we find that a known variational bound to the mutual information is tractable to estimate within SBI, that it bounds the average Kullback-Leibler divergence between surrogate and posterior, and propose to use it for model comparison on intractable inference tasks.
Comparing to a tractable posterior with estimates of exactness Assessments of approximate posterior quality are available when samples can be drawn from both the posterior θ ∼ p(θ | x) and the approximation θ ∼ q(θ | x). In the deep learning-based SBI literature, exactness is measured as a function of computational cost, usually simulator calls. We investigate this with NRE-C in Section 3.3.
Based on the recommendations of Lueckmann et al.
[44] our experimental results are measured using the Classifier Two-Sample Test (C2ST) [17,40,41]. A classifier is trained to distinguish samples from either the surrogate or the ground truth posterior. An average classification probability on holdout data of 1.0 implies that samples from each distribution are easily identified; 0.5 implies either the distributions are the same or the classifier does not have the capacity to distinguish them.  Importance sampling diagnostic An accurate likelihood-to-evidence weight transforms the data distribution into the likelihood by p(x | θ) = p(x)r(x | θ). Since NRE necessitates simulator access, we can test the ratio estimator by training a classifier to distinguish unweighted p(x | θ) samples from weighted p(x)r(x | θ) samples, wherer implies an estimate. Indistinguishability between samples implies either that the approximate ratio is accurate for parameter θ or that the classifier does not have sufficient power to find predictive features. Issues with classification power can be detected by assessing the classifier's ability to distinguish p(x) from p(x | θ). The performance can be visualized in a receiver operating curve (ROC) or measured by the area under the curve (ROC AUC). This diagnostic has been used for ratio estimators before [12,30] but it comes from training models under covariate shift [62]. It is particularly appealing because it does not require samples from p(θ | x).
Durkan et al. [16] do not mention this diagnostic in their paper, but due to its intrinsic bias NRE-B does not fulfill the identity necessary for this diagnostic to hold at optimum. The unknown factor that depends on x implies p(x | θ) = p(x) exp •g w (x | θ). We provide empirical evidence of this issue in Figure 3. Although NRE-B accurately approximates the true posterior, it demonstrably fails the diagnostic. Given the limited options for verification of SBI results, this presents a major problem by significantly limiting the trustworthiness of NRE-B on any problem with an intractable posterior. In Appendix B, we show that the unrestricted NRE-B-specific c w (x) bias means approximating p(x | θ) with normalized importance weights will not solve the issue.
Mutual information bound Selecting the surrogate model most-similar to the target posterior remains intractable without access to p(θ | x o ). Nevertheless, practitioners must decide which surrogate should approximate the posterior across training and hyperparameter search. Unfortunately, the validation losses between different versions of NRE and different K and γ settings are not comparable.
A good heuristic is to choose the model which minimizes the Kullback-Leibler divergence on average over possible data p(x). In Appendix D, we prove the relationship between I(θ; x), the mutual information with respect to p(θ, x), our models' variational bound I The non-negativity of all terms in (15) implies I(θ; x) ≥ I  w (θ; x) best satisfies our heuristic. We propose to approximate −I (0) w (θ; x) with Monte Carlo using held-out data as a metric for model selection during training and across hyperparameters. The expectation values average over p(θ, x), p(θ), and p(x). We can sample from all of these distributions in the SBI context. Since the second term in (16) computes the average log partition function, our metric can compare NRE-B-based surrogates to NRE-C-based ones. However, the metric comes with the normal challenges of estimating the log partition function which can be very expensive or high variance in some cases. We go into more depth in Appendix D, including mentioning the relevance to Neural Posterior Estimation [16,21,42,53]. While the application to SBI is novel, bounds on the mutual information have been broadly investigated for contrastive representation learning and mutual information estimation [4,22,23,25,58,71].
Empirical expected coverage probability For a candidate distribution to qualify as the posterior, integrating over data must return the prior. A measurement that follows from calibration to the prior is called expected coverage probability. Expected coverage probability can be estimated with samples from p(θ, x) and any amortized SBI method. Although important, ability to compute this metric does not distinguish NRE-C. We refer the interested reader to Hermans et al. [31]. We note that popular sequential techniques generally render this diagnostic inapplicable, with exceptions [10,47,48].

Experiments
We perform experiments in three settings to measure the exactness of surrogate models under various hyperparameters and training regimes. Section 3.1 aims to identify whether data or architecture is the bottleneck in accuracy by drawing from the joint with every mini-batch. The next experiments aim to determine how to optimally extract inference information given a limited amount of simulation data. Section 3.2 leverages a cheap prior by drawing new contrastive parameters with every mini-batch. Finally, Section 3.3 applies the commonplace training regime for deep learning-based SBI literature of fixed data and bootstrapped contrastive parameters. In this setting, we also benchmark NRE-C on ten inference tasks from Lueckmann et al.
On all hyperparameter searches we consider three simulators from the simulation-based inference benchmark, namely SLCP, Two Moons, and Gaussian Mixture [44]. SLCP is a difficult inference task which has a simple likelihood and complex posterior [21,56]. Parameters are five dimensional and the data are four samples from a two dimensional random variable. Two Moons is two dimensional with a crescent-shaped bimodal posterior. Finally, the Gaussian Mixture data draws from a mixture of two, two-dimensional normal distributions with extremely different covariances [3,63,64,70].
Our surrogate models are parameterized by one of these architectures: Small NN is like the benchmark with 50 hidden units and two residual blocks. Large NN has 128 hidden units and three residual blocks. We use batch normalization, unlike the benchmark. We compare their performance on a grid of γ and K values. We report post-training results using the C2ST, and mid-training validation loss for insight into convergence rate. We generally use residual networks [28] with batch normalization [33] and train them using adam [37]. We also run NRE-B with the same architecture for comparison. To compare with NRE-B we set the number of total contrastive parameters equal.
What does fair comparison mean in our experimentation? We compare models across fixed number of gradient steps. This implies that models with more classes, i.e., greater K, evaluate the classifierin-training on more pairs at a given training step than a model with fewer classes. An alternative which we do not apply in this paper: Vary the number of gradient steps, holding the number of pair evaluations fixed, i.e., a model with higher K sees the same number of pairs as a model with lower K but the model with lower K has been updated more times. We leave this analysis for future work.

Behavior with unlimited data
What is responsible for inaccuracies of the surrogate model when training has saturated? (a) the amount of training data (b) the flexibility of the model? In this section we provide new simulation and parameter data with every mini-batch and train until saturation, thereby eliminating the possibility of (a). The newly drawn mini-batch parameters Θ are bootstrapped for the contrastive pairs.
The setting is similar to REJ-ABC where simulations are drawn until the posterior has converged. The results of this study will provide a baseline to compare with limited-data results and help us understand how the deep learning architecture's limitations are affected by our introduced hyperparameters. The results are reported in the top row of Figure 5.
The trend is that increasing the number of contrastive examples helps NRE-B and NRE-C. There is not a clear trend with respect to γ in the C2ST. Study of the detailed validation losses in Appendix C reveals that high γ is associated with higher variance validation loss during training, but here it does not seem to strongly affect convergence rate. Saturation is reached before the maximum number of epochs have elapsed and Large NN converges to a better result. Although obfuscated in the averaged results in Figure 5, the result is obvious on SLCP, the most difficult inference task that provides an opportunity for a more flexible architecture to improve.
We argue that generally the performance bottleneck has to do with the number of contrastive parameters, network flexibility, or training details rather than the amount of data. Despite saturating validation losses, the C2ST improves based on these factors. Intuitively a more flexible model could continue to extract information from this unlimited set of jointly drawn (θ, x) pairs. This has some consequence because network flexibility may limit performance in the benchmark case. Appendix C contains more detailed results.

Leveraging fast priors (drawing theta)
In practice, the simulator is often slow but one can often draw from the prior at the same pace as training a neural network for one mini-batch. Our goal is to understand how our hyperparameter are affected by this setting and whether it is valuable to use this technique in practice. To our knowledge, this training regime has not been explored in the deep learning-based SBI literature.
Initially we draw a fixed set of around 20,000 samples (θ, x) ∼ p(θ, x). For every mini-batch during training, we sample all necessary contrastive parameters from the prior. For each term inˆ γ,K we take the same batch of contrastive parameters and reshuffle them, thereby equalizing the number of samples seen by NRE-B and NRE-C, up to bootstrap. The averaged C2ST results from the three inference tasks are reported in the middle row of Figure 5.
The resulting estimators are markedly less sensitive to the number of contrastive parameters than in Section 3.1. The number of network parameters has a positive effect, although again this is best seen on the hardest task, SLCP, in the Appendix C. Empirically, lower values of γ improve performance with Large NN; however, there is still high variance and this result may be from noise. In this setting, γ = 1 has a very small negative effect on convergence rate compared to unity, as seen in the validation loss plots in Appendix C. Since training has saturated, we claim this result implies drawing contrastive parameters from the prior helps extract the maximum amount of information from fixed simulation data x-without being strongly affected by the number of contrastive parameters K. w (θ; x), for the SLCP task estimated over the validation set versus training epochs using (a) NRE-B and (b) NRE-C with various values of γ and K, a Large NN architecture, and fixed training data. The bound permits visualization of the convergence rates and pairwise comparison across models. When K is fixed, small γ slows convergences but leads to better optima than large γ. When γ is fixed, small K slows convergence but leads to better optima than high K. Unlike when computing the validation loss directly, this metric does not exhibit a γ and K dependent bias as seen in Figures 7, 8, and 9.

Simulation-based inference benchmark
In this section, we assume the traditional literature setting of limited simulation and prior budget. Once we've selected hyperparameters based on a grid search of a subset of the SBI benchmark, we perform the entire benchmark with those hyperparameters. (mid) A fixed number of dependent (θ, x) pairs were drawn, but p(θ) was sampled at every mini-batch during training. In this regime, K has a smaller effect. (bot) The training data is completely fixed. Contrastive parameters are drawn in a bootstrap from the mini-batch. On the problems with fixed simulation data x, higher K improves accuracy and small γ with larger architectures slightly improves performance. The effects of the architecture are more clearly seen on difficult problems like SLCP in Appendix C.
Hyperparameter search To compare to the previous two sections, the amount of training data was fixed such that each epoch was comparable to Section 3.1, this amounts to about 20,000 samples. We first inspect the C2ST results in the last row of Figure 5. The sensitivity to K appears less than in Section 3.1 but more than in Section 3.2. Smaller γ slightly improves Large NN's performance but the noise makes this result uncertain. The larger network performs better, see Appendix C.
We do a Monte Carlo estimate of our proposed metric, denoted −Î , to show performance on the SLCP task as a function of epochs in Figure 4. For both NRE-B and NRE-C, increasing K tends to positively affect the convergence rate. Generally, with fixed γ and high K, peak performance is negatively impacted, in contrast to results from Figure 5 using the C2ST. Generally, with fixed K, small γ slows convergences but leads to better optima than large γ. When γ is fixed, small K slows convergence but leads to better optima than high K. See Figure 12.
Our take-away is that γ ∈ [0.1, 10] has a small effect, otherwise learning can become unstable or very slow. γ = 1 has a good compromise on high convergence rate without sacrificing too much performance. Generally, we saw improved performance on the C2ST by increasing contrastive parameters K. However, larger K also increases the magnitude of the normalizing constant, and extreme values reduced the performance on −Î (0) w (θ; x), see Appendix D. Since the benchmark compares C2ST, we optimized our architecture based on that metric. Due to bootstrapping parameters from the batched θ, the maximum K is B/2 without reusing any θ to compare with x. This is considering both terms in our loss function.
Benchmark We performed the benchmark tasks from Lueckmann et al. [44] with NRE-C. Our architecture and hyperparameters deviate based on the results of the last paragraph. We applied the largest number of computationally practical contrastive parameters, namely K = 99, and set γ = 1. Since we identified potential architecture bottlenecks in Sections 3.1 and 3.3, we used Large NN for our results. The averaged results are visible in Table 1. Sequential estimates are prepended with an S, all other estimates are amortized. We trained a five amortized estimators per task with five seeds, then we predicted the posterior for each of the ten observations x 1 , . . . , x 10 . Lueckmann et al. [44] trained an estimator for every x with a new seed. The reported C2ST for our method is averaged across ten posteriors per seed with five seeds, i.e., fifty C2ST computations per task. The other methods are averaged across one posterior per seed with ten seeds, i.e., ten C2ST computations per task. NRE-C performed better than NRE-B at all budgets and generally well among the amortized algorithms. The high performance at 10 5 simulation budget implies that our design has high capacity and may scale well. Sequential methods naturally lead to high exactness since they are tailored to the specific observation x o , but they come with drawbacks since they cannot practically perform empirical expected coverage tests [31,48] nor bound the mutual information, as discussed in Section 2.2. Detailed results per task can be found in Appendix C. We introduced a generalization of NRE-A and NRE-B which discourages the NRE-B specific c w (x) bias term and produces an unbiased estimate of the likelihood-to-evidence ratio at optimum. This property has implications for the importance sampling diagnostic, which is critical to practitioners. We suggested using a variational bound on the mutual information as a model-selection metric that could replace the C2ST averaged over several observations. It is significantly more practical since it does not require access to the ground truth posterior; however, it remains a lower bound so we don't know the overall quality of our estimator. Further discussions and derivations related to the mutual information are available in Appendix D.

Conclusion
In the context of NRE-C, we found that increasing the number of contrastive parameters K improved the C2ST in most cases, at the price of increasing the normalizing constant. Setting γ ≈ 1 was generally the fastest to saturate and lower values seemed to slightly improve performance with the Large NN and fixed simulation data. These results indicate that better performance can be achieved by using NRE-C with these hyperparameters than the other version of NRE. When both training parameters θ and simulation data x are unlimited, the architecture size plays an important role in the quality of surrogate model. We tried drawing contrastive parameters from the prior, a commonly available practical use-case, and found that it damped the effect of K on the C2ST. According to the C2ST, for highest accuracy with a fixed budget set K large, up to the mini-batch size divided by 2, and γ ≈ 1.0 for the best convergence rate. When the prior is sampled for every mini-batch, applying a higher K may be less valuable. In situations where the normalizing constant should be very close to unity, i.e., for practitioners who want to run diagnostics, drawing θ from the prior with every mini-batch and using a smaller K and γ is best.
Broader Impact The societal implications of SBI are similar to other inference methods. The methods are primarily scientific and generally lead to interesting discoveries; however, one must be careful to use an accurate generative model and to carefully test empirical results. Model mismatch and untested inference can lead to incorrect conclusions. That is why we emphasize the importance of the diagnostics in our paper. This nuance can be missed by practitioners doing inference in any field; however, special care should be taken when producing inference results that may be used for making decisions in areas like predicting hidden variables that describe human behavior or determining what factors are responsible for climate change. This list is non-comprehensive and not specific to SBI.
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] The code will be made available via a link in the supplemental material. It is a python package with example calls and self-explanatory installation instructions. The details will be shown in Section C. Results will be released upon acceptance. (b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] Additionally to the hyperparameter search which is the primary experimental result, seen in Section 3, we included more details in Appendix C. (c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [Yes] The hyperparameter searches indeed included error bars in Section 3. We did not include error bars in the computation of the diagnostic since we emphasize that it is made possible by our method, not its statistical properties. The details of the benchmark, including some uncertainty in the last experimental section are shown in Appendix C. (d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] This is discussed in Appendix C.
(a) If your work uses existing assets, did you cite the creators? [Yes] We cited the benchmark [44] several times in the paper. Also the sbi package [68] (b) Did you mention the license of the assets? [Yes] The license is mentioned in Appendix A (c) Did you include any new assets either in the supplemental material or as a URL? [Yes] I will link to my code in the supplemental material as discussed above. Upon acceptance results will be made available. (d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] The publicly available data was used in accordance with the aforementioned license in [44]. There was no need to get consent or mention it in the paper. (e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] The data is generated by toy mathematical models which abstractly relate parameters to data. It is highly unlikely any culture would find them offensive. 5. If you used crowdsourcing or conducted research with human subjects... A Relationship to other SBI methods Sampling in SBI All NRE (and SBI) algorithms require samples from the joint distribution p(θ, x). NRE additionally requires samples from the product of marginals p(θ)p(x). Sampling the joint with a simulator requires drawing from the prior θ ∼ p(θ) then passing that parameter into the simulator to produce x ∼ p(x | θ). Sampling the product of marginals is simple, just take another sample from the prior θ ∼ p(θ) and pair it with our simulation from before. Then we have (θ, x) ∼ p(θ, x) and (θ , x) ∼ p(θ)p(x). In practice, we refer to this operation as a bootstrap within a mini-batch where we take θ from other parameter-simulation pairs and reuse them to create samples drawn from the product of marginals p(θ)p(x). NRE-C sometimes requires more θ samples, often represented as Θ. These can be generated by merely concatenating several samples from p(θ). Some SBI methods, e.g., sequential methods, sometimes replace the prior with a proposal distributionp(θ).
Amortized and sequential SBI Recently, significant progress has been made in SBI, especially with so-called sequential methods that use active learning to draw samples from the posterior for a fixed observation- A core aspect of our paper focus on the effects of the NRE-B ratio estimate biased by c w (x) at optimum. Ma and Collins [45] also investigate the effects on the partition function when applying a binary or multi-class loss variant for estimating conditional energy-based models. Due to the similarity in the loss functions, they exhibit a similar bias in their partition function. In order to correct this bias they estimate c w (x) directly. We did not attempt to do this for our problem, so we do not have intuition about the effectiveness of such an approach to the likelihood-based diagnostic in NRE. Although, since c w (x) is completely unconstrained, it may be quite difficult to estimate. We believe that this would be an alternative direction for future work.
Durkan et al. [16] emphasize the connection between NRE and contrastive learning in their paper which created a framework such that Neural Posterior Estimation (NPE) and NRE-B can be trained using the same loss function. The addition of an independently drawn (Θ, x) set means that the NRE-C framework is not trivially applicable for computing the normalizing constant on atoms necessary for their sequential version of NPE.
Software Our experiments used sbi [68], which is released under an AGPL-3.0 license, and sbibm [44], which is released under an MIT license. They implement various neural SBI algorithms and benchmark problems respectively.

A.1 Full derivation of other NRE methods using our framework
Here we present a derivation of the previous works NRE-A [30] and NRE-B [30] in our framwork.
NRE-A To estimate r(x | θ), Hermans et al. [30] introduce an indicator variable y which switches between dependently and independently drawn samples. We have conditional probability Each class' marginal probability is set equally, p(y = 0) = p(y = 1). Dropping the NRE-A subscript, the probability that (θ, x) was drawn jointly is encoded in the another conditional probability NRE-A estimates logr(x | θ) with neural network f w with weights w. Training is done by minimizing the binary cross-entropy of q w (y = 1 | θ, x) := σ • f w (θ, x) relative to p(y, θ, x). For B samples, where is bootstrapped from within the mini-batch to produce θ (b) . NRE-A's ratio estimate converges to f w = log p(x | θ) p(x) = log r(x | θ) given unlimited model flexibility and data NRE-B Durkan et al. [16] estimate r(x | θ) by training a classifier that selects from among K parameters Θ := (θ 1 , . . . , θ K ) which could have generated x. In contrast with NRE-A, one of these parameters θ k is always drawn jointly with x. Let y be a random variable which indicates which one of K parameters simulated x. The marginal probability p(y = k) := 1/K is uniform. That means defines our conditional probability for parameters and data. Bayes' rule reveals a conditional distribution over y, dropping the NRE-B subscript, therefore . (21) NRE-B estimates logr(x | θ) with a neural network g w . Training is done by minimizing the cross entropy of q w (y = k | Θ, where θ k ) over B samples. In our parameterization and given unlimited flexibility and data, g w (θ, x) = log p(θ | x) p(θ) + c w (x) at convergence. The extra term enters because the optimal classifier for (22) need not be normalized.

B Theoretical Arguments
We present first the proof of convergence for NRE-C. Afterwards, we discuss the properties of estimated importance weights in NRE-B.
We have: = arg min f E p(z,y) log which is minimized, when KLD = 0, thus: which implies that for p(z)-almost-all z: So we get with the definition of q(y | f * , z): Dividing the latter by the former gives for k = 1, . . . , K: implying for k = 1, . . . , K and p(z)-almost-all z: This shows the claim.
We used the symbol KLD to imply the Kullback-Leibler Divergence. The proof in Lemma 1 is slightly more general than necessary for our typical case. All our functions f k are typically the same, namely the evaluation of a neural network with weights w. Rather than searching for the function f which minimizes the objective, we search for the weights, but these are equivalent. Finally to make everything fit, we set z := (θ, x).

B.2 Properties of the importance sampling diagnostic on biased ratio estimates
In Section 2.2 we discuss the importance sampling diagnostic, in which the estimated ratio is tested by comparing weighted samples from x ∼ p(x) to unweighted samples from x ∼ p(x | θ). The estimated ratio from NRE-C is merely exp(h w (θ, x)) and is therefore not restricted to have the properties of a ratio of probability distributions, except at optimum with unlimited flexibility and data. The importance sampling diagnostic is designed to test whether the estimated ratio is close to having these properties for a fixed θ. One way to improve ratios estimated by NRE-C is to compute the normalizing constant Z w (x), which should be close to one, and replace the ratio with this "normalized" version.
The unrestricted nature of the NRE-B-specific c w (x) bias means that the normalization constant Z w (x) does not have to be close to one. Further, we show that the unrestricted bias means that the normalizing constant for two NRE-B ratio estimators, which are equivalent in terms of their loss function (22), is not unique, i.e., the estimate is ill-posed. This property means that normalizing NRE-B will not, in general, produce an estimator which passes the diagnostic.
Consider two ratio estimates with the relationshipr 1 (x | θ) =r 2 (x | θ)/C(x) where C is an positive function of x resulting from the aggregation of the exponentiated bias. Given N samples x n ∼ p(x) with weights w 1 (x n ) =r 1 (x n | θ) and w 2 (x n ) =r 2 (x n | θ), we can compute the importancenormalized weights bȳ However, if we substitute this constant back into our expression, we find that the weights do not agreē Therefore, normalization does not "protect" against scaling bias introduced by a function in x. NRE-B does not penalize functional biases like C(x), so even ratios considered optimal by NRE-B can easily fail the diagnostic. Meanwhile, NRE-C encourages terms like C(x) towards one. That has the effect of making the C(x) ≈ 1 drop out of the normalized importance weights, thus making performance on the diagnostic indicative of a better ratio estimate, given enough classifier flexibility.

C Experimental Details
Computational costs Our experiments were performed on a cluster of Nvidia Titan V graphics processing units. The primary expense was the hyperparameter search within Sections 3.1, 3.2, and the first part of Section 3.3. The run of those experiments took about 50,000 gpu-minutes total considering all of our current data. Since there were several iterations, we multiply this by three to estimate total compute. The next expense was the computation of sbibm in Section 3.3 that took about 24,000 gpu-minutes. Together these equal about 2880 gpu-hours. An estimate for the total carbon contributions corresponds to 311.04 kgCO 2 . Luckily, our clusters run completely on wind power offsetting the contribution. Estimations were conducted using the MachineLearning Impact calculator presented in [39].
Architecture and training The centerpiece of our method is the neural network h w which is trained as a classifier. The hyperparameter choices here were fairly constant throughout the experiments. Any hyperparameters for training or architecture which were consistent across experiments are listed in Table 2. Hidden features and number of RESNET blocks depend on the architecture. Large NN uses three resnet blocks with 128 features, while Small NN uses two resnet blocks and 50 features. In this section, we trained many neural networks with different hyperparameter settings on three different tasks from Lueckmann et al. [44]. We trained both architectures, Large and Small NN, on a grid of γ and K values with NRE-B and NRE-C. No matter whether we were training with unlimited draws from the joint, prior, or fixed data we designed an epoch such that it has 20 training mini-batches and 2 validation mini-batches. For fixed initial data and prior this corresponds to a simulation budget of 22,528. The mean validation losses per-epoch for these networks is visible in Figures 7, 8, and 9.
For one specific problem, we also computed validation loss at a fixed γ and K no matter the training regime in an effort to produce comparable validation losses, specifically γ = 1 and K = 1 aka the NRE-A regime. This training setting reverses some of the trends we've seen when validated on the same loss as training data. We did not apply it further since this could bias results. See Figure 11. To compare, we also plot the validation loss on the sample problem but using the γ and K each model was trained with. See Figure 10, note that the bias depending on γ and K is clearly visible (just as in the other validation loss plots).
Once the networks were trained, we drew samples from ten per-task posteriors based on ten predefined observations x i with i ∈ 1, 2, . . . , 10. This leveraged the amortized property of NRE-C. Samples were drawn on these problems using rejection sampling. Once the samples were drawn, they were compared to ground truth posterior samples from the benchmark with the C2ST. The detailed plot which shows per-task behavior is available in Figure 6.    where colors indicate different contrastive parameter counts. The plot shows the convergence rates of each model. Since validation loss is a function of K and γ, the relative performance of models is not comparable. A grid search of K and γ indicates that increasing K leads to earlier convergence at fixed γ. With K fixed, γ < 1 has a negative effect on convergence rate and γ > 1 is ambiguous. where colors indicate different contrastive parameter counts, i.e, K + 1. A fixed validation loss was used, namely 1,1 . Although the validation loss is now comparable so that we can see the performance of the different classifiers on the same task, some of the convergence rate trends disagreed with SLCP in Figure 10. Naturally, the classifier trained to distinguish two samples often performed best on this task.

C.2 Simulation-based inference benchmark
We trained NRE-C with γ = 1.0 and K = 100 using Large NN for 1000 epochs on all of the tasks from the simulation-based inference benchmark [44]. In contrast with our previous hyperparameter search experiments, we used their Markov-chain Monte Carlo sampling scheme which applies slice sampling [51]. The details results are presented in Table 3 and Table 4. There were three tasks which did not succeed with this sampling procedure, there we used rejection sampling.
Our training slightly diverged from the benchmark. We enumerate the ways below.
• Large NN is a bigger architecture than the one which produced their reported NRE-B values.
In fact, the NRE-B version has the same settings as Small NN; however, we found that larger networks performed better in our search and chose the bigger one for that reason. • In the benchmark, networks are trained with early stopping but we always trained for 1000 epochs. We selected the network which performed best on the validation loss, using the same values for K and γ in the validation loss as in the training. • Even though NRE-B is an amortized method, the benchmark trained a network for every observation. Instead, we leveraged the amortization properties of NRE-C and trained a single network which drew samples from each approximate posterior given by each of the ten observations.
Task details We provide a short summary of all of the inference tasks in the SBI benchmark by Lueckmann et al. [44].
Bernoulli GLM This task is a generalized linear model. The likelihood is Bernoulli distributed. The data is a 10-dimensional summary statistic from an 100-dimensional raw vector. The posterior is 10-dimensional and it only has one mode. Bernoulli GLM Raw This is the same task as above, but instead the entire 100-dimensional observation is shown to the inference method rather than the summary statistic. Gaussian Linear A simple task with a Gaussian distributed prior and a Gaussian likelihood over the mean. Both have a Σ = 0.1 · I covariance matrix. The posterior is also Gaussian. It is performed in 10-dimensions for the observations and parameters. Gaussian Linear Uniform This is the same as the task above, but instead the prior over the mean is a 10-dimensional uniform distribution from -1 to 1 in every dimension. Gaussian Mixture This task occurs in the ABC literature often. Infer the common mean of a mixture of Gaussians where one has covariance matrix Σ = 1.0·I and the other Σ = 0.01·I. It occurs in two dimensions. Lotka Volterra This is an ecological predator-prey model where the simulations are generated from randomly drawn initial conditions by solving a parameterized differential equation.
There are four parameters that control the coupling between the generation and destruction of both prey and predators. The priors are log normal. The data is a twenty dimensional summary statistic. SIR An epidemiological model simulating the progress of an contagious disease outbreak through a population. Simulations are generated from randomly drawn initial conditions with a parameterized differential equation defining the dynamics. There are two parameters with a log normal prior. The data is a ten dimensional summary statistic. SLCP A task which has a very simple non-spherical Gaussian likelihood, but a complex posterior over the five parameters which, via a non-linear function, define the mean and covariance of the likelihood. There are five parameters each with a uniform prior from -3 to 3. The data is four-dimensional but we take two samples from it. It was introduced in [56]. SLCP with Distractors This is the same task as above but instead the data is concatenated with 92 dimensions of Gaussian noise. Two Moons This task exhibits a crescent shape posterior with bi-modality-two of the attributes often used to stump MCMC samplers. Both the data and parameters are two dimensional. The prior is uniform from -1 to 1.

D Mutual Information
Estimating the mutual information is closely related to estimating the likelihood-to-evidence ratio. We reference various bounds on the mutual information and show how NRE-C can estimate them numerically. These bounds obey the variational principle and might be a practical candidate for validating the performance of SBI methods for scientific purposes, i.e., when the ground truth posterior is intractable. The approximation of the mutual information could synergize with other diagnostics like empirical, expected coverage testing and the importance sampling diagnostic. See section 2.2.
Mutual information and expected Kullback-Leibler divergence If p(θ | x) denotes the true posterior and p w (θ | x) an approximate posterior then the quality of that approximation can be measured via the (forward) Kullback-Leibler divergence: However, in the SBI-setting we have only access to the likelihood p(x | θ) via samples. Since we have also access to the prior p(θ) (analytically and) via samples we can sample from the joint p(θ, x). So using the expected Kullback-Leibler divergence is more tractable in the SBI setting to measure the discrepancy between true and approximated prior: In all considered cases in this paper the approximate posterior p w (θ | x) is given by: wherer w (x | θ) comes from a (trained) neural network and Z w (x) denotes the normalization constant. With the above notations we get for the expected Kullback-Leibler divergence the expression: where I(θ; x) is the mutual information w.r.t. p(θ, x). Since we aim at minimizing the expected Kullback-Leibler divergence we implicitely aim to maximize our mutual information approximation: which can be estimated via Monte-Carlo by sampling i.i.d. θ n , θ n,m ∼ p(θ), x n ∼ p(x | θ n ), n = 1, . . . , N , m = 1, . . . , M and then compute: exp(h w (θ n,m , x n )) .
Since in all mentioned methodsr w (x | θ) is meant to approximate the ratio p(θ | x) p(θ) and estimating the normalizing constant is expensive a naive alternative to approximate the expected Kullback-Leibler divergence is by naively plugging the unnormalized distributionq w (θ | x) :=r w (x | θ)p(θ) into the above formula and using the estimate: While this is justified for the training objectives for NRE-A and NRE-C, which encourage a trivial normalizing constant Z w (x) ≈ 1 at optimum, the same is not true for NRE-B, which leads to an additional non-vanishing, possibly arbitrarily big, bias term: Another way to address the normalizing constant is the use of the Kullback-Leibler divergence that also works for unnormalized distributions p(z), q(z): KLD(p(z) q(z)) := log p(z) q(z) + q(z) p(z) − 1 p(z) dz, which is always ≥ 0 with equality if p(z) = q(z) for p(z)-almost-all z.
This gives for the expected Kullback-Leibler divergence between the posterior p(θ | x) and the unnormalized approximate posteriorq w (θ | x) :=r w (x | θ)p(θ): We see that the normalizing constant Z w (x) re-appears, but with a different term. Similar to before the above can be used for another mutual information approximation given by: (exp(h w (θ n,m , x n )) − 1) .
Note that since log(r) ≤ r − 1 we always have the inequalities: showing that I (0) w (θ; x) leads to a tighter approximation to the mutual information I(θ; x) than I (1) w (θ; x). The procedures above require estimating the normalizing constant or the partition function. Generally this can be quite expensive and may require techniques like Nested Sampling [65]; however, it is tractable with Monte Carlo on problems with parameters within low dimensional compact regions.
Bounds on the mutual information There is a connection between the training objective of NRE-B and a multi-sample lower bound on the mutual information [58], as noted in Durkan et al. [16]. Contrastive learning has been explored for estimating the mutual information by Van den Oord et al. [71] also discussed by Belghazi et al. [4]. The bounds we define above are also discussed in detail by Poole et al. [58], although we find computingÎ  w (θ; x) to be tractable within SBI-although potentially expensive.
We attempted to train a ratio estimator by minimizingÎ (1) w (θ; x) directly on fixed data. Our preliminary experiments found that the C2ST was near unity on the SLCP task and other works claim that this mutual information bound has extremely high variance. We therefore ended our investigation. w (θ; x) versus the C2ST on the SLCP task with every point corresponding to a different set of values for γ and K. C2ST scales from 0.5 to 1.0 with 0.5 implying that the classifier could not distinguish the approximate posterior from the ground truth. Just like in the fixed data regime in the main experiment, see Figure 6, we found that on this task a lower γ improved the C2ST. We also find that the mutual information is correlated with the average C2ST across 10 pieces of data, but the mutual information has the practical advantage that we can bound it without knowing the ground truth posterior. The C2ST requires being able to sample from the ground truth posterior. This data represents the same set of experiments as in Figure 4, Figure 12, and Figure 14.