Calibrating Neural Simulation-Based Inference with Differentiable Coverage Probability
Authors: Maciej Falkiewicz
Naoya Takeishi
Imahn Shekhzadeh
Antoine Wehenkel
Arnaud Delaunoy
Gilles Louppe
Alexandros Kalousis
Venue: NeurIPS 2023 poster (2023)
Keywords: simulation-based inference
inverse problem
bayesian inference
uncertainty quantification
generative modeling
TL;DR: Regularization method for neural simulation-based inference that uses a differentiable computation of posterior's expected coverage probability.

Abstract: Bayesian inference allows expressing the uncertainty of posterior belief under a probabilistic model given prior information and the likelihood of the evidence. Predominantly, the likelihood function is only implicitly established by a simulator posing the need for simulation-based inference (SBI). However, the existing algorithms can yield overconfident posteriors (Hermans *et al.*, 2022) defeating the whole purpose of credibility if the uncertainty quantification is inaccurate. We propose to include a calibration term directly into the training objective of the neural model in selected amortized SBI techniques. By introducing a relaxation of the classical formulation of calibration error we enable end-to-end backpropagation. The proposed method is not tied to any particular neural model and brings moderate computational overhead compared to the profits it introduces. It is directly applicable to existing computational pipelines allowing reliable black-box posterior inference. We empirically show on six benchmark problems that the proposed method achieves competitive or better results in terms of coverage and expected posterior density than the previously existing approaches.

Calibrating Neural Simulation-Based Inference with Differentiable Coverage Probability
Maciej Falkiewicz!?2 Naoya Takeishi Imahn Shekhzadeh Antoine Wehenkel
Arnaud Delaunoy Gilles Louppe" Alexandros Kalousis
CCmputer Science Department, University of Geneva HHES-SO/HEG Genève
Thhe University of Tokyo "RIKEN UUiiverstty Liège
{maciej falkiewicz, imahn. shekhzadeh, alexandros .at.anrrrsaallussss Saeesehheeeeeccc ch
ntakeōg. - ecc. uutokyo..ac.jpp
Vantoine.wehenkee a. d.leaunny,, g. g.Louppe/@uliege.
Abstract
Bayesian inference allows expressing the uncertainty of posterior belief under probabilistic model given prior information and the likelihood the evidence. Predominantly, the likelihood function is only implicitly established by simula- tor posing the need for simulation- based inference (SBI). However, the existing algorithms yield overconfident posteriors (Hermans et al., 2022) defeating the inaccurate. We
whole purpose of credibility if the uncertainty quantification propose to include calibration term directly into the training objective of the neural model in selected amortized SBI techniques. By introducing a relaxation of the classical formulation of calibration error we enable end-to-end backpropagation. The proposed method is not tied to any particular neural model and brings mod- erate computational overhead compared the profits introduces. It is directly applicable to existing computational pipelines allowing reliable black-box posterior inference. We empirically show on six benchmark problems that the proposed method achieves competitive or better results in terms of coverage and expected posterior density than the previously existing approaches.
Introduction
Inverse problems [46] appear in variety of fields of human activity, including engineering (2]
geology [8) medical imaging [29], economics (1I. and particle physics M. In some domains,
uncertainty quantification is required for results to be used in the real-world [5I]. A particularly
well-suited tool for this task Bayesian inference, where the procedure of updating prior beliefs (knowledge) with observed evidence is formally described in the framework of probability theory. In order to compute the posterior probability according to Bayes' rule, the likelihood function has to be evaluated. Often, it is only implicitly established by simulator which given input parameters generates (possibly stochastically) outputs. Such a situation requires the use of likelihood- techniques, also referred to simulation-based inference (SBI) (
There exist two major classical approaches that address the problem Bayesian inference in the likelihood-free setting. In the first one, simulated data are used to estimate the distribution low-dimensional summary statistics with a nonparametric method like kernel density estimation [14]. Next, the approximate density is used the problem into likelihood-bard based inference. In the second approach, termed as Approximate Bayesian Computation (ABC) 40 38 48
of
likelihood surrogate turning
47 data is simulated upfront, but instead for every newly appearing observation (instance the inverse problem), an independent sampling- rejection procedure with simulator in-the-loop is run. Given some proposal distribution, candidate
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

parameters are accepted if the simulator's output is in sufficient agreement with the observation. While ABC effectively sorts out the lack of explicit likelihood, it generally costly, especially when the underlying simulator is computationally expensive.
Deep Neural Networks have been increasingly utilized in SBI to alleviate these costs, by acting as universal and powerful density estimators or surrogating other functions. The following approaches of neural SBI can be identified: Neural Likelihood Estimation (NLE) [35], where the missing likelihood function is approximated based on the simulator; Neural Ratio Estimation (NRE)
19
15 32), where
the likelihood- ceroooddtt--viieeee to-evidence ratio (other ratios have also been considered) approximated; Neural Posterior Estimation (NPE) [4 30 I8), where the learned density estimator directly solves the inverse problem for a given observation. In the last category the score-based methods [4] [6] are emerging, where the posterior is implicitly expressed by the learned score function. All of the approaches mentioned above amortize the use of simulator, however, inference with NLE and NRE still requires the use of Markov chain Monte Carlo 77,, while NPE can be fully amortized with the use of appropriate generative models 47
The question of reliability is raised with the use of approximate Bayesian inference methods. Many works have been devoted to inference evaluation techniques in the absence of the ground-truth posterior reference [45] 53 27 (28]. Lueckmann et al. [3I propose an SBI benchmarking protocol taking into account the effect of the simulation budget, with the observation that the assessment's conclusions will strongly depend on the metric chosen. Hermans et al. [20] suggest relying on coverage analysis of posterior's multi-dimensional credible regions, showing that numerous SBI techniques yield overconfident posteriors. In follow-up work [1]] they introduce a balancing regularizer for training NRE, and show empirically that it leads to the avoidance of overconfident solutions. Very recently, an interesting study by Delaunoy et al. [12] extends the applicability of balancing regularizer to the NPE approach. An alternative perspective is studied in Dey et al. 13 where an additional regression model for posthoc re-calibration based on Probability Integral Transform is learned for the outputs of a pre-trained inference model.
In this work, we introduce and empirically evaluate a new regularizer that directly targets the coverage of the approximate posterior.
applicable to NRE and NPE based on neural density estimators such as Normalizing Flows [26 36). As an additive term, it is easy for practitioners to use as a plug-in to existing pipelines. The additional computational overhead introduced by the regularizzr is justified in view of the advantages it brings.
The rest of the work is structured follows. Section 2 introduces the formal description of the studied problem and presents the coverage analysis that is the motivation for the presented method. Section 3 describes the proposed method and differentiable formulation of the regularizer that allows brasee learning. Section shows the empirical results of applying the method on a wide range of benchmark problems, as well as providing apper-parameter sensitivity analysis. In we provide summarize our work with conclusions and potential future directions.
section 5
2 Background
In this work, we consider the problem of estimating posterior distribution p(0]zz) of model parameters 8, given the observation x. We approach the problem from a frree [43 perspective, where the model's likelihood p(x[00) cannot be evaluated. Instead, there exists simulator that generates observations x based on the input parameters 0. We can prepare a dataset ((0,, x; x;)); by sampling parameters from a prior p(0) and running them through the simulator. In that way, the problem of finding the approximate posterior marked as p(0]z) turns into a supervised learning problem. However, it is important to note that usually there is not only single generating parameter 0 for given æ possible. Instead, for each observation multiple parameter values are plausible [23]. These values form a posterior distribution [0]zz) under prior distribution p()).
Assessing the performance of an SBI technique is challenging due to the aon-availability of the ground-truth posterior p(H]æ)) in the general case. Theoretically, methods such as the rejection- sampling ABC [44) are proven to converge to p(0(zr)) under an infinite simulation budget but often to achieve good approximation. Avoiding this costly procedure is exactly what motivates the introduction of different techniques. Therefore, it common to resort to performance metrics evaluated on test datasets sampled from the simulator.
require large number of simulation runs

---

If a Bayesian inference method allows p(0)xz) evaluation, its performance can be assessed by looking the expected value of the approximate log posterior density of the data-generating parameters (nominal parameters) Epce.a) llgg [0[==)] estimated over a test dataset. In domains such as natural sciences, it is also accurate uncertainty quantification that matters, i.e. reliable posterior density estimation given the available evidence in data. It has been shown in Hermans et al. [20] that approaches yielding high expected posterior density values can give overconfident posteriors, and therefore, verifying the quality of credible regions returned by inference engines as well as designing more reliable ones essential. In this work, we define a regularizer suitable for training NRE and NPE based the Simulation-Baate paper self-contained, in the rest of Section 2 we will introduce SBC formally, adapting the definitions and lemmas from Lemos et al. [27]. and making the link with the coverage analysis from Hermans et al. 22
Calibration (SBC) diagnostic of Talts al. [45]. To make the
2.1 Coverage analysis
Coverage analysis offers way to evaluate the quality of an approximate posterior. Given an observation and an estimated distribution, one can assess the probability that a certain credible region [22] of the parameter space contains the nominal parameter value found in the dataset. For a conditional probability distribution, a credible region at level - a is a region of the support to which the parameter used to generate the conditioning observation belongs with 1 - probability, called the credibility level. There are infinitely many credible regions at any given level, and coverage analysis requires particular way of choosing one. In this work, will consider the Highest Posterior Density Region (HPDR).
Definition (Highest Posterior Density Region).
θ HPDR (1 a) == arg
mig min O s.t.
p(0)r)d0 =
(1)
p(0
Joee
where [+| : Se du is the volume of a subset O of the parameters' space, under some measure μ.
Next, we introduce the central concept of coverage analysis.
Definition 2 (Expected Coverage Probability). The Expected Coverage Probability (ECP) over joint distribution p(0, x) for a given O (1 - a) is
ECP HAAE p(o)z) (1 (1-a)) = Ep(..)) [1[0 M(l()) (18)) (1-
HPDR
(2)
where [[:] is the indicator function.
While the expectation in eq. can be easily estimated by using a set of i..dd samples from the joint distribution p(0, æ), the challenge of finding p(0]=) (1 - a) remains non-trivial In general, directly following eq. ① requires solving an optimization problem. This would be prohibitive to do for every instance during training, thus it crucial to benefit from the two lemmas we introduce below.
Lemma 1. A pair (0", xx) and distribution p(0)z)) uniquely define an HPDR:
p(0)=) (1 ourDo((p 0" ,,xx))))))))""))))"))))) = {0 (00)==)) > p001)=))) |x²)},
(3)
where
CHPpR( xx)00,x))== = |(000)1(((0))) < $(00)>"))00
(4)
See proof in Appendix A Lemma 2.
ECP
e
HPDR p(0)z)
(1 (1-a)) = Ep(o.x) [1[AHPDR((2_ 0) > a]]
(5)
See proof of Lemma in Lemos et al. 27
If ECP HPDR (1 a)) is equal to 1 - a the inference engine yields posterior distributions that are (()
calibrated at level a, or simply calibrated the condition is satisfied for all a. The ground-truth posterior p(0 z) is calibrated, but another such distribution is the prior p(0) discussed in Hermans
3

---

et al. (20] and Lemos et al. 127 Therefore, coverage analysis should be applied in conjunction with the standard posterior predictive check. Formally, we aim for:
p" (0 arg max pp(e..) [lgg((0)>)) s.t. Va ECP p(8)z)
HPDR p(0)) (1 a)) =1-a
(6)
In Hermans et al. the equality sign in the condition in eq. (6 can be relaxed and replaced with ECP @HHPRR p(o)x)
20 it has been pointed out that for practical applications often the requirement of
(1 a)
above equal the level. Such posterior distribution is called conservative. We define the calibration error as
M
M (1 0:) - ECP ( (1 az))
Σ
1=1
and the conservativeness error as
M
max
(1 (1-os) -ECP (
p(oto) (1-))()))))))))) p(Olz)
HPDR
(8)
M
where {a is a set of credibility levels of arbitrary cardinality such that 4, < a, < 1; we exclude a and construction. For motivating graphical example of coverage analysis of an overconfident distribution see appendix D
because calibvatiio/ccossrraaiveee
these levels is achieved by
3 Method
We propose a method supporting the training of neural models underlying SBI engines such that the resulting distributions are calibrated or conservative against the real p(0,z) the data. In addition to the standard loss term in the optimization objective function, we include a regularizer corresponding differentiable relaxations of the calibration (conservativeness) error as defined in eq. (eq. (8).
In the remainder of this section, we will show how to efficiently estimate the expected coverage probability during training and how alleviate the arising non con-continuities allow gradient-based optimization. After a few simplifications, the regularizer will become a mean square eтror between the sorted list of estimated oHpdr (p, 0, .TSs and a step function on the (0, 1) interval.
The proposed method is based on the Monte Carlo (MC) estimate of the integral in eq.
atPOR (p. 00,x*))
Σ
< p(0*)=)))
L
where 0.,.... 0g are samples drawn from p(0 (r). The entity defined in eq. 9 is a special case of what is called rank statistic in Talts et al. 45 (with the 1/L normalization constant skipped there). The non-normalized variant is proven to be uniformly distributed over the integers [0. L. for p(0)x:) (f) V(, ~ p(0,z)) x), which is the accurate posterior defined in Lemos et al. [27]. Taking -> oo and applying the normalization (as done in eq. (9) one shows that distribution of dRPPR (p, 8", x²) converges to standard continuous uniform distribution (proof in appendix A.2)
We return to coverage analysis with the following remark:
Remark 1. An inference engine yields calibrated posterior distributions and only if it is character- ized by uniformly distributed ourPR (f. 0, 2..)s.
See proof in Appendix A
Following remark we will enforce the uniformity of the normalized rank statistics which coincides with necessary, though not sufficient, condition for the approximate posterior to be accurate. The one-sample Kolmogorov-Sniron test provides a way check whether a particular sample set follows some reference distribution. We will use its test statistic to check, and eventually enforce, the condition. An alternative learning objective to directly minimize the calibration (conservativeness) error at arbitrary levels with ECP computed following eq. (

---

3.1 One-sample Kolmogoroy-Smirnox test
Let B be random variable taking values 3 E B. Given sample set B = [B,]i = I...., N). the
one-sample Kolmogorov-S Smimnov test statistic [5] for r.v. B is defined as:
DN sup (x (x(3)- BEB
(10)
where Fy is the empirical cumulative distribution function (ECDF) of the r.v. B, with regard to the sample set B [T]]. F is the reference cumulative distribution function (CDF) against which we compare the ECDF.
In the proposed method, we instantiate the KS test statistic for r.v. CHIPDR (p. 0, xx), where (0,z)~ ~ p(0,x), thus Fx is its ECDF, which we estimate from a sample set Á == (åfIPDR (p. 0g (2,, 0.) p(x, f), = N.}. The reference CDF, F, is that of the stan- dard uniform distribution because we want to check and enforce the uniformity of the elements of A.
By the Glivenko-Cantelll theorem 44 we know that DN 0 almost surely with N - if samples follow the F distribution. The claim the opposite direction, that is if Dy + 0 then samples follow F, is also true and we show that in appendix AA3] As result, minimizing DN will enforce uniformity and thus posterior's calibration. We can upper-boume the supremum in eq. (IO) with the sum of absolute errors (SAE), but instead, stabilize gradient descent, of squared errors (MSE), which shares its optimum with SAE. Minimizing the one-sample KS test statistic is just one of the many ways of enforcing uniformity. We can deploy any divergence measure between the sample set and the uniform distribution.
the mean
3.2 Approximating DN
In practice cannot evaluate the squared errors the complete support of the r.v,, but only finite subset. There are two ways of evaluating this approximation which differ by how we instantiate (a) and F(o.).
Direct computation In the first, we compute (a = E, [aa; S al/N summing over a; e A and use the fact that F(a) a for the target distribution. To approximate Dx we aggregate the squared error over an arbitrary set of ag values: (Fx(()) regularizer means that we need to backpropagate through Fx (ox). Such an approach can be seen as directly aiming for calibration at arbitrary levels an i.e. applying eq. ( to evaluate the calibration error.
ax)². Including this term as a
Sorting-bast based computation Alternatively, we can sort the values in A to get the ordered sequence (a.li := ....N). Given the sorted values, ag and their indices, the Fx(a,) values are given directly by the indices normalized by the sample size (i/N), and by the corresponding target distribution F(a,)= ag which leads to the (,i/ ax)² approximation. Under this approach, the gradient will now flow through F(a.). This is the approach that we use in our experiments.
3.3 Importance Sampling integration
In the definition of CHPDR (p, 0* x²) in eq. as well in its MC estimate åfiPDR (p, ("'x") defined in eq. (g) we sample 0,s from the approximate posterior 6(0)z). Optimizing the model constituting the approximate posterior requires backpropagating through the sampling operation which can be done with the reparametrization trick (24). However, the reparametrization trick is not always applicable; for example, in NRE an MC procedure is required sample from the posterior [99 Here
instead self-normalized Importance Sampling (IS) to estimate ÖHPRR (6, 0, x) and lift the requirement
of sampling from the approximate posterior. We will be sampling from known proposal distribution and only evaluate the approximate posterior density. We thus propose to use IS integration [e.g., 17 estimate the integral in eq. (4) in the following way:
L.IS a (p. 0* x²) Σt= p(0, [T00) < 800 [²))
(11)
Σラ 1(0
5

---

Algorithm Computing the regularizer loss with calibration objective.
Require: Data batch {(0t, model p(0)))), number of samples I, proposal distribution (())
Ensure: Regularizer's loss R
for I N do
2:
p: - x(,,
3: for j 1 to L do
4:
8: I(0)
sampling from proposal distribution
5:
p + p00 (x,)
6:
end for
@HPPR (f. 0; vi) &- Σi Ei / (2:1lp//:0 <p<
7:
D eq.
8:
end for
(arli ,,..., N) t sort(( (a (p. 0, z;)li 1, N.)
9:
10:
R N E (i/N - as)²
D- the second term in eq. 12
The d,s are sampled from the proposal distribution I(0). For uniform proposal distribution, the I(0,) terms are constant and thus cancel out. By default, we suggest using the prior as the proposal distribution. While IS is not needed in the case of NPE, we use it also for this method to keep the regularizer's implementation consistent with the NRE setting.
3.4 Learning objective
So far we have shown how to derive regularizer that enforces a uniform distribution of the normalized rank statistics that characterizes calibrated model as shown in remark T This also necessary condition for the accuracy of the posterior as shown in Talts et al.
(55 but not sufficient one, as exemplified by the approximate posterior equal the prior for which the regularizer trivially obtains its minimum value. Thus the final objective also contains the standard loss term and is a weighted sum of the two terms:
N
Σ log (0g((()))s)) +
1 N Σ (index (LLIS (f, 0c, ar. /N IN- dHHDR (ò. 0,, mi))
LLIS
(12)
N
which we will minimize with Adam W optimizer. The first term in which we use for the NF-based NPE (or its corresponding cross-entropy loss for the NRE) and the second is the uniformity enforcing regularizer. The index(-) operator gives the index of its argument in the sequence of sorted elements of L.IS (6. 0,,,,) conservativeness we accept aHPDR (p, 0*, .xx) ) to first-order stochastically dominate [39] the target uniform distribution, thus the differences in the uniformity regularizer are additionally activated with a rectifier max(:; 0) analogously to eq. 8). In appendix A.1 we sketch a proof of how minimizing such regularizer affects the learned model. The two main syper-parameters of the objective function are the number of samples L we use to evaluate eq. (TT), and the weight of the regularizer À in D ww provide the pseudo-code compute the uniformity regularizer.
is the likelihood loss
NV. When regularizing for
eq. (12). In algorithm
To allow backpropagation through indicator functions, we use a Straight- Thrrohgh Estimator 4 with the Hard Tanh (37) backward relaxation. In order to backpropagate through sorting operation its differentiable implementation from Blondel et al. [6].
Experiments
In our experiments]] we basically follow the experimental protocol introduced in Hermans al. 20 for evaluating SBI methods. We focus on two prevailing amortized neural inference methods, i.e. NRE approximating the iikelihood-to-evideme model. In the case of NRE, there has been already regularization method proposed [II] and we include it for comparison with our results, and mark BNRE. In the main text of the paper, focus on the conservativeness objective. For the calibration objective, the results can be found in
ratio and NPE using conditional NF as the underlying
'The code is available //ttttss://giiiii

---

SLOP
Weinberg
Spatial SIR
M/G/1
Lotion Volterra
Gravitational Waves
1.0
NRE
0.5
1.0
DNRE
0.5
0.0
1.0
CANRE
0.5
0.0
NPC
1.0
Expected coveaae
CaNPE
0.5
0.0
Crediddity
1024
4095
16384
63036
Calibrated astimater
2048
8192
32708
131022
Figure 1: Evaluation on six benchmark problems (for each column) in the form of coverage curves estimated on set of 10k test instances. the curve is on (above) the diagonal line, then p(0) [x)) is calibrated (conservative). Expected coverage of HPDRs at 19 evenly spaced levels on the [0.05; 0.95] interval is estimated following eq. (2) by solving an optimization problem to find HPPRR p(f) (1 a) for
every test instance. The average of retraining five times from random initialization is shown. Top there marked as NRE, BNRE, and CaINRE (ours) show the NRE trained without regularization, with balancing regularizer, and the proposed regularizer with conservativeness objective, respectively. The bottom rows marked NPE and CaINPE show the NPE trained without regularization and with the proposed regularizer with conservativeness objective, respectively. Best viewed color.
appendix C where report positive outcomes only for part of the studied SBI problems. As the proposal distribution for IS, we use the prior p(f). In the main experiments, we set the weight of the regularizer À to 5, and the number of samples L to 16 for all benchmarks. In section 4.3 we provide insights into the sensitivity of the method with respect to its hyper-parameters.
4.1 Principal experiments
In fig. [ present cross-sectional analysis of the results of the proposed method, based six
benchmark SBI problems described in detail in Hermans et al. (20), for each considering eight differ- ent simulation budgets. Expected coverage of HPDRs at 19 evenly spaced levels on the [0.05; 0.95] interval is estimated following eq. (2) by solving an optimization problem to find p(0 (1 a) for every test instance. We intentionally use a different ECP estimation method than in the regularizer to avoid accidental overfitting to the evaluation metric. The top three rows show empirical expected coverage curves the test set obtained by training the same model with three different objectives NRE: classification loss; BNRE: classification loss with balancing regularizer 11: CaINRE: classifi- cation loss with the proposed regularizer. The results for CaINRE show that the proposed method
7

---

SLCP
Weinberg
Spatial SIR
M/G/1
Lotta Volterza
Gravitational Waves
0.22
22.5
0.4
+30
10'
10
10'
Simulation budget
NRE
BNRE
CANRE
NPE
CaINPE
Priucc
Expected value of the approximate log posterior density of the nominal parameters Bp(0..)) [log p(0[z)] of the SBI approaches in fig. [ estimated over 10k test instances. The solid line the median over five random initializations, and the boundaries of the shaded area are defined by the minimum and maximum values. The horizontal black line indicates the performance of the prior distribution. Best viewed in color.
Figure 2:
Computation on CPU
Computaaton GPU
Figure 3: The computational overhead of the proposed method applied to NRE the SLCP problem. The Y-axis indicates how many times training with the regularizer lengthens relative to training without it. For every number of samples experiments were conducted on both CPU (left chart) and GPU (right chart). Averaged over five runs, error bars show a double standard deviation range. Best viewed color.
Doathnter cmn ootboel
Typer printiter
systematically results in conservative posteriors - with SLCP simulation budget 1024 standing out, due to too weak regularization weight A. Moreover, in fig. 2 it can be seen that CaINRE typically outperforms BNRE in terms of (e))).)) [log ((0[x)), being close to the performance NRE (which is not conservative), and demonstrating less variance across initializations.
The bottom two rows fig. [ show the same type of results for NPE-based models, with CaINPE being trained with the proposed regularizer. For all budgets in the Gravitational Waves benchmark, CaINPE does not meet the expected target, but this can easily fixed by adjusting the regularization strength X. In fig. 2 one finds that CaINPE sometimes yields higher Ep(e,a) [log pa [f]]z) than the non-regularized version, which a counterintuitive result. We hypothesize that by regularizing the model we mitigate the classical problem of over-fitting on the training dataset.
4.2 Computation time
In fig. 3 we show how applying the proposed regularizer affects the training time in the case of NRE for the SLCP benchmark, with the distinction to CPU and GPU. As a reference, we take the training time of the non-regularized model. While on CPU, the training time grows about linearly with L, on GPU due to its parallel computation property the overhead is around x2, with some variation depending on the simulation budget and the particular run. However, this is only true long as the GPU memory is not exceeded, while the space complexity is of order ((Ndim(0)A)
where N the batch size. Comparing with fig. 2 we conclude that applying the regularizer can result in similar (p(8,=) [logr((0(z>) (00(=)) as in the case of increasing the simulation budget (what extends training time) while keeping comparable computational time, and avoiding over-confidence. Moreover, in situations where the acquisition of additional training instances very expensive even impossible, the proposed method allows for more effective use of the available data at a cost sub-linear in L.
4.3 Sensitivity analysis
In fig. we show the sensitivity of the method with respect to its hyper- pppammppraeeters on the example of NRE applied for the Weinberg benchmark. In fig. 4a] the analysis of regularization strength indicates tradcoff between the level of conservativeness and Epce.)) [log ()] Therefore, choosing the
8

---

Fiffet
Fllert experted
Filre: Chompon MIC
Elleet 1o
Foimnoel purtenee Granto
0094
Chavsige 3DC
0.07
0.00
500
0.10
(Libmz))
Ohdiidiig lnve
Hippr pecsmeter
Hyper genancteetterr..
(a) Sensitivity analysis with respect to regularizer's 's weight A.
Fleers on portetioos
Bleert on expprtted coveray
Bllect Cinerage ATC
Bliet Bes.. lon, AlMeel
6.00
1.3
0.04
CU33
1.0
11.00
-0.0
0.00
-0.15
-1100
0.5
2.0
0.0
10
XRR (158572)
Hypeeppanntte
(b) Sensitivity analysis with respect the number of samples in the regularizer.
Figure 4: Sensitivity analysis study of CaINRE on the Weinberg benchmark with the smallest simulation budget of 1024 training instances with respect to in (a) and L in (b). The yellow color in each chart shows the results of non-regularized NRE on the highest simulation budget of 131072 training instances. Results in the first two charts show the performance of ensemble model built from five randomly initialized members. The vertical line in the first chart marks the nominal parameter value. Points in the third chart indicate the average AUC over five randomly initialized models, with two times the standard deviation marked with the shaded area. The shaded area and the error bars in the fourth chart span between minimum and maximum results over the five randomly initialized models with the median indicated by point (line). Best viewed in color.
right configuration for practical applications is multi-objective problem. The lowest possible value that yields a positive AUC (A = the figure) should be used to maximize predictive power as measured by Ep(e,2) [log p(0|==)] From fig. 4b
read that increasing the number of samples in the regularizer brings the coverage curve closer to the diagonal while increasing Ex(e.2) [log p(0)=)) This suggests that more accurate estimation of conservativeness errors during training allows better focus on the main optimization objective.
Conclusions and future work
In this work, we introduced a regularizer designed avoid overconfident posterior distributions, which is an undesirable property in SBI. It is applicable in both the NRE and NPE settings, whenever the underlying model allows for amortized conditional density evaluation. the evaluation metric based on Expected Coverage Probability, used to assess whether the posterior distribution is calibrated or conservative, depending on the target. We formally show that minimizing the regularizer to 0 error under infinite sampling guarantees the inference model to be calibrated (conservative).
derived directly from
We empirically show that applying the regularizer in parallel to the main optimization objective systematically helps to avoid overconfident inference for the studied benchmark problems. We admit that calculating the additional term adds computation cost at training time, but this can be mitigated with the use of modern hardware accelerators well suited for parallel computation.
The presented method is introduced in the context of simulation based inference, but is suitable for application in general conditional generative modeling, where calibration (conservativeness) can for example help mitigate mode collapse [50].
As future direction, we identify the need evaluate the proposed method in problems where the parameters' space has a high dimension. Moreover, currently, the method restricted to neural models acting density estimators, but not general generative models that only allow for sampling from the

---

learned distribution. A sampling- based relaxation of the regularizer would extend its applicability to a variety of models including Generative Adversarial Networks [33. score/diffusion-bene -based models S2 and Variational Autoencoders 24
Acknowledgments and Disclosure of Funding
We acknowledge the financial support of the Swiss National Science Foundation within the MIGRATE project (grant no. 209434). Antoine Wehenkel and Arnaud Delaunoy recipients of an F.R.S.- FNRS fellowship and acknowledge the financial support of the FNRS (Belgium). Naoya Takeishi was supported by JSPS KAKENHI Grant Number JP20K 19869. The computations were performed at the University of Geneva on "Baobab" and "Yggdrasil" HPC clusters.
References
[1] G. Aad, T. Abajyan, B. Abbott, J. Abdallah, S. A. Khalck, A. A. Abdolalim, R. Aben, B. Abi, M. Abolins, O. AbouZeid, et al. Observation of particle in the scarch for the Standard Model Higgs boson with the ATLAS detector at the LHC. Physics Letters B, 716(1):1-29 2012.
[2] V. Akbarzadeh, J.-C. Lévesque, C. Gagné, and M. Parizeau. Efficient sensor placement optimization using gradient descent and prohabilistic coverage. Sensors, 1455511555--15555, 2014.
[3] M. A. Beaumont, J.MM Cornuct, J.-M. Marin, and C. P. Robert. Adaptive approximate Bayesian computa- tion. Biometrika, 96(4):983 990, 2009.
Y. Bengio, N. Léonard, and A. Courville. Estimating or propagating gradients through stochastic neurons conditional computation. arXiv preprint arXiv: 2013.
V. W. Berger and Y. Zhou. Kolmogoro--mirrin Test: Overview. John Wiley & Sons, Ltd, 2014.
[6] M. Blondel, O. Teboul, Q Berthet, and J. Djolonga. Fast differentiable sorting and ranking. In International Conference on Machine Learning, pages 950-959. PMLR, 2020.
[7] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng. Handbook of Markov chain Monte Carlo. CRC press, 2011.
[8] A. Cividini, L. Jurina, and G. Gioda. Some aspects of 'characterization' problems in gcomechanics. In International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, volume 18, pages 487-503 Elsevier, 1981.
[9] Cranmer. Brehmer, and Louppe. frontier simulation-baat based inference. Proceedings of the National Academy of Sciences, 117(48):30055-3006,, 2020.
[10] F. Dekking, C. Kraaikamp, H. P. Lopuhaä, and L. E. Meester. A Modern Introduction to Probability and Statistics: Understanding why and how, volume 488, page 219. Springer, 2005.
[11] A. Delaunoy, Hermans. Rozet, A. Wehenkel, and G. Louppe. Towards reliable simulation- based inference with balanced neural ratio estimation. In A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, editors, Advances in Neural Information Processing Systems, 2022.
[12] A. Delaunoy, B. K. Miller, P. Forré, C. Weniger, and G. Louppe. Balancing simulation cause inference for conservative posteriors. arXiv preprint arXiv:2304. 10978, 2023.
[13] B. Dey, D. Zhao, J. A. Newman, B. H. Andrews, R. Izbicki, and A. B. Lee. Calibrated predictive distributions via diagnostics for conditional coverage. arXiv preprint arXiv:2205. 2022.
[14] P. Diggle and R. Gratton. Monte Carlo methods of inference for implicit statistical models. Journal of the Royal Statistical Society, Series B (Methodological), 16(2))193-22) 227, 1984.
[15] C. Durkan, I. Murray, and G. Papamakarios. On contrastive learning for likelihood-free inference. In International Conference on Machine Learning, pages 2771-2781. PMLR, 2020.
[16] T. Geffner, G. Papamakarios, and A. Mnih. Score modeling for simulation- inference. In NeurIPS 2022 Workshop on Score-Based Methods, 2022.
[17] J. Geweke, Bayesian inference in econometric models using Monte Carlo integration. Econometrica: Journal of the Econometric Society, pages 1317-1339, 1989.
10

---

[18] D. Greenberg, M. Nonnenmacher. and J. Macke. Automatic posterior transformation for likelihood-fiees inference. In International Conferencee on Machine Learning, pages 2404-2414. PMLR, 2019.
[19] J. Hermans, V. Begy, and G. Louppe. Likelihond-free MCMC with amortized approximate ratio estimators. In International conference on machine learning, pages 4239-4248. PMLR, 2020.
[20] J. Hermans, A. Delaunoy, F. Rozet, A. Wehenkel, V. Begy, and G. Louppe. A crisis in simulation-based inference? Beware, your posterior approximations can be unfaithful. Transactions on Machine Learning Research, 2022.
[21] J. L.. Horowitz. Ill-posed inverse problems in economics. Annu. Rev. Econ., 6(1):21-51, 2014.
[22] R. Hyndman. Computing and graphing highest density regions. The American Statistician, 50(2): 120-126, 1996.
[23] E. T. Jaynes. Prior information and ambiguity in inverse problems. Inverse problems, 14:151-166, 1984.
[24] D. P. Kingma and M. Welling. Auto-encoding variational Bayes. arXiv preprint arXiv: 14, 2013.
[25] G. Klambauer, Unterthiner, A. Mayr, and S. Hochreiter. Self-normalizing neural networks. Advances in neural information processing systems, 30, 2017.
[26] Kobyzev, S. J. Prince, and M. A. Brubaker. Normalizing flows: An introduction and review of current methods. IEEE transactions on pattern analysis and machine intelligence, 43(11):3964 3979, 2020.
[27] P. Lemos, A. Coogan, Y. Hezaveh, and L. Perreault-Levasseus Sampling- baaeed accuracy testing of posterior estimators for general inference. A. Krause, E. Brunskill, K. Cho, Engelhardt, S. Sabato, and J. Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 19256 19273. PMLR, 23 23-29 Jul 2023.
[28] Linhart, A. Gramfort, and P. L.. Coelho Rodrigues. Validation Diagnostics for SBI algorithms based on Normalizing Flows. In NeurlPS 2022 Machine Learning and the Physical Sciences Workshop, pages 1-7, New Orleans, United States, Nov. 2022.
[29] T. Liu, P. Spincemaille, L. De Rochefort, B. Kressler, and Y. Wang. Calculation of susceptibility through multiple orientation sampling (COSMOS): method for conditioning the inverse problem from measured magnetic field map to susceptibility source image in MRI. Magnetic Resonance in Medicine: An Official Journal of the International Sociery for Magnetic Resonance in Medicine, 61(1):196-200 2009.
[30] J.-M. Lueckmann, P. J. Gonçalves, G. Bassetto, K. Öcal, M. Nonnenmacher, and J. H. Macke. Flexible statistical inference for mechanistic models of neural dynamies. In Proceedings of the 31st International Conference on Neural Information Processing Systems, pages 1289-1299. 2017.
[31] J.-M. Lueckmann, J. Boelts, D. Greenberg, P. Goncalves, and J. Macke. Benchmarking simulation-based inference. In International Conference on Artificial Intelligence and Statistics, pages 343-351. PMLR, 2021.
[32] B. K. Miller, C. Weniger, and P. Forré. Contrastive neural ratio estimation. Advances in Neural Information Processing Systems, 35:3262-33 0278,, 2022.
[33] Z. Pan, W. Yu, X. Yi. A. Khan, F. Yuan, and Y. Zheng. Recent progress on generative adversarial networks (GANs): A survey. IEEE access, 7:36322-36333, 2019.
[34] G. Papamakarios and I. Murray. Fast s-free inference of simulation models with Bayesian conditional density estimation. Advances in neural information processing systems, 29, 2016.
[35] G. Papamakarios, D. Sterratt, and Murray. Sequential neural likelihood: Fast likelihood free inference with autoregressive flows. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 837-848. PMLR, 2019.
[36] G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. The Journal of Machine Learning Research, 2211):2617- 2021.
[37] J. Pemnington, S. Schocnholz, and S. Ganguli. The emergence of spectral universality in deep networks. In International Conference on Artificial Intelligence and Statistics, pages 1924-1932. PMLR, 2018.
[38] J. K. Pritchard, M. T. Seielstad, A. Perez-Lezaun and M. W. Feldman. Population growth of human Y chromosomes: a study of Y chromosome microsatellites. Molecular biology and evolution, 16(12):
1791-1798, 1999.
11

---

[39] J. P. Quirk and R. Saposnik. Admissibility and measurable utility functions. The Review of Economic Studies, 29(2):140-146, 1962.
[40] D. B. Rubin. Bayesianly justifiable and relevant frequency calculations for the applied statistician. The Annals of Statistics, 12(4):1151 1172, 1984.
[41] L. Sharrock, J. Simons, S. Liu, and M. Beaumont. Sequential neural score estimation: Likelihood- free inference with conditional score based diffusion models. arXiv preprint arX1v:2210.04872, 2022.
[42] S. A. Sisson, Y. Fan, and M. M. Tanaka. Sequential Monte Carlo without likelihoods. Proceedings of the National Academy of Sciences, 10666)1760-1765, 2007.
[43] S. A. Sisson, Y. Fan, and M. Beaumont. Handbook of approximate Bayesion computation. CRC Press, 2018.
[44] M. Sunnáker, A. G. Busetto, E. Numminen, J. Corander, M. Foll, and C. Dessimoz. Approximate Bayesian computation. PLOS Computational Biology, 9(1):1-10, 01 2013.
[45] S. Talts, M. Betancourt, D. Simpson, A. Vehtari, and A. Gelman. Validating Bayesian inference algorithms with simulation-based calibration. arXiv preprint arXiv: 1804.06788, 2018.
[46] A. Tarantola. Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial and Applied Mathematics, 2005.
[47] J. M. Tomezak. Deep generative modeling. Springer, 2022.
[48] T. Toni and M. P. H. Stumpf. Simulation- based model selection for dynamical systems in systems and population biology. Bioinformatics, 26(1):104-11) 110, 10 2009.
[49] F. Topsse. On the Glivenko-Cantelli theorem. Zeitschrift für Wahrechheimiichheisshhccii und verwandte Gebiete, (4(3):239-250 1970.
[50] L. Vaitl, K. A. Nicoli, S. Nakajima, and P. Kessel. Gradients should stay on path: better estimators of the reverse- and forward KL divergence for normalizing flows. Machine Learning: Science and Technology. (4):045006, oct 2022.
[51] V. Volodina and P. Challenor. The importance of uncertainty quantification in model reproducibility. Philosophical Transactions of the Royal Society A, 379(2197):20000711 2021.
[52] L. Yang. Z. Zhang, Y. Song, S. Hong. R. Xu, Y. Zhao, W. Zhang, B. Cui. and M.-H. Yang. Diffusion models: comprehensive survey of methods and applications. ACM Comput, Surv, sep 2023. Just Accepted.
[53] D. Zhao, N. Dalmasso, R. Izbicki, and Lee. Diagnosties conditional density models and Bayesian inference algorithms, Uncertainty Artificial Intelligence, pages 1830-1840. PMI., 2021.
12

---

A Proofs
Lemma 1. pair (0",x") and distribution p(0)z)) uniquely define an HPDR:
UPDR p(0 (1 C.HPDR(P, 0- ,x")) = {0 1(0)=²) > p }}
(3)
where
(H,DR((,000x")) p(01=)111(01=") < 1(0"12")]<t
(4)
Proof Lemma
aHPDR (0, x²))
@NDOR p())))) (1- (2.00 x)d0
1010001=1)20000000 p(0)x)d0-
/ ((01)11p((107 > p(0 (13)
while
/ ;(01=)11((0"") > p(0*==))] / < -10')))
(14)
because p(0)z) ) probability density function.
Remark An inference engine yields calibrated posterior distributions if and only if it is character- ized by uniformly distributed currra(p, (p, 0,
Proof of Remark We start by observing that RHS of eq. 5
be reformulated as Ep(0,x) [[ OHPPR (Ø, 0, x) < ai]] which is estimated using MC in the following way:
N Σ
1[1 OHPPR 01, < 11-a].
(15)
and for a calibrated model it is equal to - a. The empirical ECP expressed above, takes the form of the ECDF of 1 @HPDR (p, z;, 1 Thus, following Glivenko-Canteli theorem conclude aHpDr (p, 0,, 2i) has a standard uniform distribution because is bounded in [o, 1] interval by construction. Thus, the same holds for @HPDR (f, 2;)
that for a calibrated model,
The opposite direction is trivial show.
A.1 Stochastic dominance rank statistics over (0 1| as over-confidenee aversion
Let us consider the following minimization objective
DN sup max (ax - F(a),0)
(16)
which leads first-order stochastic dominance of (x(a)) (ECDF of (OHPDR (ý, 0,, N over F(a) (target standard uniform distribution). At the minimum Dx = O we know that
Va [ 0, < al/N.
(17)
But the LHS can reformulated as ECDF of infinitely many samples generated with the ground-trule distribution. Without the loss of generality, assume the number of samples is equal to N both sides.
V2
L|HHFOR (p. 0,, & aj [OHPDR( 2,)) & a).
(18)
13

---

With equality between LHS and RHS, this means calibrated p(0(z) which is special case of conservativeness in our formulation. Let us focus only on the inequality. Then, the condition above translates into
a 3 OHPPR (p, 3g) & a & OHPPR (z. > a,
(19)
thus
Ya 3 aHPDR (p. xg) < OHPDR (p.
(20)
Now, insert the definition of CHPPR (: xz);
3 1 ((00snn))))))) 1(0,
(((0)))22)) p(0,l,))]
(21)
where Un) are sampled from the ground truth posterior p(0|z), and O() from the approximate posterior p(0]z). The following two steps show equivalent conditions:
Va 3j 3 (n) (r2)) p(0,la,)) & ((0))))))) ((0,1,)),
(22)
Va 7, 7 (()))))))))) ((0
(23)
The condition in eq. (23 says that ground-truth posterior generates at least sample with higher density measured by the approximate posterior p(0]z) than that generated by the approximate posterior. Thus, have shown that model generating OHPDR 0, 24)5 stochastically dominating standard uniform distribution has an aversion to over-confidence.
A.2 Convergence uniform
Let us consider discrete uniform distribution over 7. + 1 evenly spaced values on the [0. 1] interval, which has the following CDF:
0,
x<0 0 x <1
F2(2)=
= |
1,
x>1
The CDF of the standard continuous uniform distribution has the following form:
x F(x)
I,
05 x <1
1, x x
For we have
[F(x) - (x)] = - = 1 [xb]]] 1/LL,
where the inequality a property of the floor
|.] operator.
Taking the limit of - oo we have
lim FL()) lim 1/L 0, L
what proofs convergence the distribution of Fz to F(x).
A.3 GC theorem inverse direction
We are given a set of iid samples a o....., ON from G distribution implicitly given by the trained model. With regularizer optimized to loss, we assume that
DN sup - F(a))
We follow the derivation from httppp//aaattt..aaaaacccaaaaaaee
(version: 2018-10-23).
We follow the derivation from
(version: 2020-08-29)
14

---

SLCP
Weinberg
Spatial SIR
M/G/1
Lotia Volterta
Gravitational Waves
NRE
0.0
0.11
0.2
0.2
0.1
BVRE
0.0
CANRE
0.1
0.0
0.1
0.2
0.2
0.1
NPC
0.0
0.2
CMNPE
Coverage AUe
0.1
0.00
0.1
0.2
130
budges
Figure 5: AUC scores of the coverage curves presented in fig. 1 Points indicate the average, while
the shaded area spans over the standard deviation range. Best viewed in color.
From the Glivenko-Canter theorem, we have that
DN = sup (Fx()) G(a)) $ 0.
The triangle inequality satisfied for any norm, so:
sup Fx(a) F(a)) + sup [N(()) G(a)) sup GG(a) - F(a)),
thus
sup [G(a) F(a)) 3 0,
where the LHS does not depend on N, and we conclude that G' = F", what to be shown.
Additional conservative posterior results
In this section, we supplement the results shown in section 4.1 by presenting the coverage AUC scores in fig. 5 We show explicitly, that applying the proposed method with conservativeness objective leads to conservative posteriors. For the two non-conservative cases CaINRE for SL.P,, and CaINPE for Gravitational Waves we expect a positive result when choosing the appropriate value of the regularization weight, already stated in the main text.
15

---

SLOP
Weinberg
Spatial SIR
M/G/1
Lotion Volterra
Gravitational Waves
1.0
NRE
0.5
1.0
CAINRE
0.5
0.0
1.0
NPC
0.5
Expeetoo contang
CANPE
Crviibilly level
10384
Calibrated
2013 S192
32768
131072
Figure 6: Evaluation on six benchmark problems (for each column) in the form coverage curves estimated on set of 10k test instances. If the curve is on (above) the diagonal line, then p [+) is calibrated (conservative). Expected coverage of HPDRs at evenly spaced levels on the [0.05; 0.95] interval is estimated following eq. (2) by solving optimization problem to find HPDR -(()a) (1 a) for
every test instance. The average of retraining five times from random initialization is shown. Top two rows marked as NRE, and CaINRE (ours) show the NRE model trained without regularization, and the proposed regularizer with calibration objective, respectively. The bottom rows marked as NPE and CaINPE show the NPE model trained without regularization and with the proposed regularizer with calibration objective, respectively. Best viewed in color.
C Regularization for calibrated posterior
In this section, we present the results experiments analogous those in section 4.1 but with the proposed regularizer configured as calibration eтrоr (see eq. (7).. Obtaining calibrated posterior task way more difficult than the case conservative one, and this what we empirically confirm.
In fig. [6] we find characteristic S-shaped coverage emerging for the regularized models. For the low credibility levels (small credible regions) posteriors become overconfident, while with the high levels (big credible regions) they become conservative, by a certain margin. Overall, this results in an approximately 0 AUC as shown in fig. S However, it is only by chance that it becomes positive (on average conservative), and negative (on average overconfident). From fig, 7 we also see that the predictive power of the regularized models measured by Emjo,)) [log p((]z)) estimated on the test degrades compared to the conservativeness setting fig.
16

---

SLCP
Weinberg
Spatial SIR
M/G/1
Lotla Voltera
Gravitational Warre
0.0
0.2
2.5
-300
Simulation budget
NRE
CaIRRE
NPE
CaINPE
Pricr
Figure Expected value of the approximate log posterior density of the nominal parameters Ep(8,a) [log p ([][=)) of the models in fig. 6 estimated over 10k test instances. The solid line is the median over five random initializations, and the boundaries of the shaded area defined by the minimum and maximum values. The horizontal black line indicates the performance of the prior distribution. Best viewed in color.
SLCP
Weinbery
Sputial SIR
M/G/11
Lot'a Voltema
Gruvitational Waws
0.7
0.1
NRC
0.0
0.1
11.2
0.2
CAUNRE
0.2
0.1
NPE
0.0
011
0.2
Contrage AUC
CaNNE
10
201
10
100
107
10'
10.
10
10
10
100
117
107
102
207
10
Simulation besigr
Figure 8: AUC scores of the coverage curves presented in fig. 6 Points indicate the average, while the shaded area spans over the two standard deviation range. Best viewed in color.
17

---

D Intuition about coverage analysis
In this section, build toy example intuitively introduce the concept of Expected Cover- age Probability and the wver-/under-conffde
that results from it. We rely on non-conditional one-dimensional distributions for clarity. In fig. | we introduce a black multimodal ground-truth distribution and a red one intentionally under-dispersed, but with the same modes, acting as the approximate distribution. We find the which consists of two segments concentrated around the modes. To mimic the real-world scenario we sample N 1000 points from the ground -ttuth distribution and use them to estimate the expected coverage probability, which is trivial in one dimension knowing the segments' boundaries. According the plan, the coverage analysis showed that the red distribution overconfident for evidence from the data.
0.9 Highest Density Region for the red distribution,
E Implementation details
For the details of the architecture of Neural Ratio Estimation models please refer to Delaunoy et al. M.. The SELU function [25] is used as an activation function in the feed-forward networks.
For the details of the architecture of Neural Posterior Estimation models please refer to Hermans et al. 20
All the models with the proposed regularizer are trained for 500 epochs with 128 batch size with AdamW optimizer with 0.001 learning rate.
Additional remarks:
We use the following PyTorch implementation of Blondel et al. [o]: https://giit gitbbb http://ddyyokeeetorheceeee We use with the default parameters.
In order to stabilize the training of regularized NRE (CaINRE) we introduced Gradient
Norm Clipping.
Training of CaINPE for the M/G/1 in the conservativeness setting with the last four simula- tion budgets was challenging due to extreme values in observations. It is reflected in the decrease of Ep(0,a) [log p [0|z)] in fig. 2 as the simulation budget increases.
When computing the regularizer the same pC ((="==") is repeatedly evaluated, therefore, implementation evaluates x= s embedding only once (per regularizer evaluation) and reuses it for all samples.
$ A I00, (xo((0:9) 0.793
Figure 9: Plots of two probability density functions of 2-mixture of Gaussians (0.7/0.3 ratio). with the red one (01 = 0.7, 22 = 0.2) being under-dispersed compared to the black one (a1 = 0.9, 02 = 0.4). The blue region marks 0.9 Highest Posterior Density Region of the red distribution. The black points on the x-axis are random samples from the black distribution. Below the graph, the empirical Expected Coverage Probability of the red distribution against N 1000 samples of the black distribution is computed. Best viewed in color.
18


Reviews:
[{"review_id": "8bn7e0vgAT", "reviewer": "Reviewer_T8qr", "summary": "The paper introduces a training algorithm for posterior distribution learning in the likelihood free setting that combats overconfident models.  The authors focus on the expected coverage probability (ECP) from Hermans et al., 2022 to measure if a model posterior is conservative.    Specifically, when its value is equal to the credibility level used in its calculation, then the posterior is calibrated.  An equivalent condition is that the distribution of credibility levels constructed from samples (lemma 1) is uniformly distributed.  With this in mind, the authors propose a regularizer that penalizes how far the distribution of these credibility levels are from a uniform distribution.   The main contribution of the paper is the algorithm used to train with this regularizer.  The authors use a variety of techniques in their algorithm including a one-sample Kolmogorov-Smirnov test to work with finite samples, differentiable sorting to implement the test, and importance sampling to choose useful samples to use.  The experimental results demonstrate that the learned models are more conservative than their unregularized counterparts while not sacrificing performance in terms of likelihood.  Furthermore, the authors showed the effect of the hyperparameters on performance and also showed the that the training algorithm can be expensive if not run on GPUs.", "questions": "- Why choose the one-sample Kolmogorov-Smirnov test to measure the divergence over others?", "limitations": "Yes, the authors acknowledge that their method is computationally expensive and may not scale to high dimensions as is.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Directly addresses the \"trust crisis in simulation-based inference\" described in Hermans et al., 2022.\n- Strong empirical results to demonstrate that the training algorithm works.\n- The paper progressed smoothly from the background to the proposed algorithm.\n- The layout of the text and plots are visually easy to digest.\n", "weaknesses": "- It wasn't too clear to me after reading the main text why having uniformly distributed credibility levels was the right thing to want.\n- Similarly, it took a bit of drawing to see how all of the values in section 2 were related to each other.  Adding something like figure 9 to the main text could greatly help introduce the background material.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer,\n\nThank you very much for taking the time to review our manuscript and for your comments. Below we would like to address the questions:\n\n- *Why choose the one-sample Kolmogorov-Smirnov test to measure the divergence over others?*\n\n   The one-sample Kolmogorov-Smirnov test was our first choice. Its test statistic is straightforward to use as a minimization objective. The differentiable relaxation was also easy to find. We did not investigate alternative approaches.\n\n   As common alternatives for the KS test, the Anderson–Darling and Cramér–von Mises tests are typically listed. Both of them require ordered samples, thus only the sorting-based (alternatively, ranking-based) computation is available. Moreover, both of the tests rely on tabulated critical values which necessitates the introduction of significance level and is less convenient to use as a minimization objective. \n"}, {"review_id": "enuysDCt7f", "reviewer": "Reviewer_Lpwd", "summary": "The authors suggest a new objective function for simulation-based inference that adds a penalty term to the “expected score” objective that is used by many other works. This penalty term encourages the resulting posterior approximations to be well-calibrated in the sense that the $1-\\alpha$ highest probability density regions contain the ground-truth parameter $100(1-\\alpha)$% of the time. The penalty term is an extension of the Kolmogorov-Smirnov test statistic for goodness of fit between a $U(0,1)$ distribution and the rank statistics of eq. 7, which are asymptotically uniformly distributed when $p(\\theta \\mid x)$ is calibrated. Minimizing this penalty term alongside the usual objective function should yield posteriors that are calibrated while substantially different from the degenerate case of the prior (which is trivially calibrated). \n", "questions": "* The method called neural posterior estimation (NPE) has been shown in other work to be equivalent to minimizing a forward KL divergence between exact and approximate posteriors (Reweighted Wake-Sleep (Bornschein and Bengio, 2015), and a similar work “Revisiting Reweighted Wake Sleep…” (Le et al., 2018)); it’s been argued in these that the resulting neural posterior estimates tend to be overdispersed or conservative as a result. 1) How does this impact the relevance of this work if the baseline is usually already conservative? 2) The NPE row of Figure 1 doesn’t seem to reflect this intuition in practice; is there any explanation for why this is so? Is NPE indeed the same as sleep-phase in reweighted wake-sleep in the experiments?\n* What are the advantages and disadvantages of using either the sorting-based computation or the direct computation method? Is the use of the latter solely due to the belief that backprop through indicator functions is somehow worse? As $\\hat{\\alpha}$ itself requires indicator functions, I suppose the backprop operation can’t be too problematic. It would be interesting to see both methods implemented in the experiments. \n* The form of eq. 9 suggests self-normalized importance sampling, but as $\\hat{p}(\\theta \\mid x^*)$ is already normalized, this doesn't seem necessary. Maybe this is a misunderstanding of notation. To my reading, though, just the numerator of eq. 9 is a valid importance sampling estimator. Is so-called “standard” or self-normalized importance sampling nonetheless used?\n", "limitations": "Generally the settings where this work may be applied are clear, but as mentioned above discussing the use of the method in amortized vs. non-amortized settings may make more clear the costs associated with with this method when $\\hat{p}(\\theta_j \\mid x^*)$ can’t be computed quickly. ", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "* The proposed penalty/regularization term is relatively lightweight computationally, and can be tacked on to many existing simulation-based inference algorithms. \n* The intuition is simple; the construction of the statistics $\\hat{\\alpha}$ is clever, and the Kolmogorov-Smirnov test statistic is well-understood.\n* The method is assessed on a variety of problems from the test suite of benchmarks provided by Lueckmann et al. (2021), and performs favorably in that the resulting posteriors are either calibrated or tend to be conservative.\n", "weaknesses": "* Although synthetic likelihood and ABC approaches are mentioned in the introduction, the proposed method seems prohibitively costly in these scenarios, and seems geared toward amortized methods only, where eq. 7 can be computed rapidly.\n* For ratio estimation in particular, the non-regularized NRE is better calibrated than the corrected version in some cases (e.g., spatial SIR), suggesting that it is at least possible the penalty term can worsen calibration.\n* In practice, it appears that attempts at calibration often result in conservative rather than calibrated posteriors; while this is likely still preferable to many users to the alternative of *not* adding this correction, it suggests a bit of a misnomer. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer,\n\nThank you very much for taking the time to review our manuscript and for your comments. Below we would like to address the questions:\n\n- *Although synthetic likelihood and ABC approaches are mentioned in the introduction, the proposed method seems prohibitively costly in these scenarios, and seems geared toward amortized methods only, where eq. 7 can be computed rapidly.*\n\n   We would like to clarify that the proposed method is intended only for amortized methods. Synthetic likelihood and ABC are mentioned only to give a more complete context of SBI, we will make that more clear in the text.\n\n\n- *For ratio estimation in particular, the non-regularized NRE is better calibrated than the corrected version in some cases (e.g., spatial SIR), suggesting that it is at least possible the penalty term can worsen calibration.*\n\n   The non-regularized NRE for Spatial SIR is not conservative for the first three simulation budgets, while CalNRE is conservative for all of the analyzed budgets which is the improvement that we intend to introduce with the proposed method. Indeed, the distance of the coverage curve from the diagonal may increase when using the proposed regularization term, but this is not a violation of the objective in eq. (6) re-formulated in L127 for conservativeness.\n\n\n- *In practice, it appears that attempts at calibration often result in conservative rather than calibrated posteriors; while this is likely still preferable to many users to the alternative of not adding this correction, it suggests a bit of a misnomer.*\n\n   We would like to clarify that the results presented in the main part of the submitted manuscript (Figures 1 & 2 + Figure 5 in the Appendix) are obtained with the objective of minimizing conservativeness error (L130), not calibration error. The result of the latter, we show in the Appendix (Figures 6-8) and find it difficult to draw any firm conclusions. Finding a calibrated approximate posterior is a very difficult task and we make no claim that the proposed method is able to systematically provide it. We will emphasize this in the revised version of the manuscript.\n\n\n- *The method called neural posterior estimation (NPE) has been shown in other work to be equivalent to minimizing a forward KL divergence between exact and approximate posteriors (Reweighted Wake-Sleep (Bornschein and Bengio, 2015), and a similar work “Revisiting Reweighted Wake Sleep…” (Le et al., 2018)); it’s been argued in these that the resulting neural posterior estimates tend to be overdispersed or conservative as a result. 1) How does this impact the relevance of this work if the baseline is usually already conservative? 2) The NPE row of Figure 1 doesn’t seem to reflect this intuition in practice; is there any explanation for why this is so? Is NPE indeed the same as sleep-phase in reweighted wake-sleep in the experiments?*\n\n   The training objective of NPE is indeed the same as the one used in sleep-phase of RWS, but with instances $(\\theta, x)$ coming from a fixed training dataset (measured in the field or sampled from a simulator) in contrast to “dream samples” from the trained generative network in RWS. So the problem solved by these two methods is fundamentally different.\n\n   The NPE row of Figure 1 is reprinted from Hermans et al., 2022 where in an intense empirical study the authors identify multiple methods (including NPE) to be over-confident (non-conservative) for a number of benchmark problems in SBI – in a limited study we confirmed their results (not reported in the submitted manuscript). In our work, we use the same benchmark problems, and with the application of the proposed regularizer are able to train models that are conservative in cases where the non-regularized models were not. We do not claim that the use of our method is necessary, but rather we propose it as a solution when over-confidence is identified. In addition, we would like to note that Hermans identifies two standard solutions in machine learning as often effective in the analyzed problems, i.e. increasing the amount of the training data and ensembling. Therefore, it is enough that in the cited works in the field of RWS, the budget for simulations during training was large enough that the problem of over-confidence did not arise. Meanwhile, in SBI, obtaining a sufficient dataset can be very expensive, or even impossible.\n\n\n- _The form of eq. 9 suggests self-normalized importance sampling, but as $\\hat{p}(\\theta | x^*)$ is already normalized, this doesn't seem necessary. Maybe this is a misunderstanding of notation. To my reading, though, just the numerator of eq. 9 is a valid importance sampling estimator. Is so-called “standard” or self-normalized importance sampling nonetheless used?_\n\n   Indeed, the self-normalized importance sampling is given in eq. (9) and also used in our experiments. For NPE the denominator is completely unnecessary because the approximate posterior is normalized by construction. However, in the case of NRE, there is no such guarantee (we only hope that this is achieved at convergence) which exposes a slight abuse of notation on our side. Therefore, self-normalized IS is advisable for NRE and to have a single implementation we use it also for NPE."}, {"review_id": "HiDzFsJ5yz", "reviewer": "Reviewer_KAKy", "summary": "The paper proposes a calibration term to be used directly in the training objective of NREs and NPEs. The paper shows that the introduction of this term achieves competitive or better results in terms of coverage and expected posterior density. ", "questions": "* Line 206: Does IS provide an improvement for NPE, even when it is not needed? \n* Since the regulariser relies on importance sampling did the authors do any experimentation on how the approach scales with dimension of $\\theta$? This was mentioned in future work, but even in the experiments in the paper the theta varies across experiments. What are the dimensions of the experiments included in this paper? And did the authors see that the higher dimensional ones performed worse with the new regulariser?\n", "limitations": "* The limitation related to computational cost is sufficiently highlighted in the work.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "* The quality of the writing and presentation is high, with the structure easy to follow. In particular the related work is clearly included in the introduction.\n* The topic is of importance as it is increasingly more common to use neural estimators in SBI applications.\n* The experimentation is sufficient to show the effect of incorporating the new regulariser. It seems to show that the new regulariser leads to slightly conservative estimators that are on average better calibrated than the baselines. \n", "weaknesses": "* The paper does not appear to have any major weakness. A few minor weaknesses that seem to exist have been appropriately described in the paper. These are the computational cost of the approach and the fact that the paper’s main metric of performance is also the same one that has been used in the regulariser. If the authors could think of a different metric to use to evaluate performance then the paper would be further strengthened. However, it is appreciated that finding a new metric would potentially be a new paper in its own right. Perhaps including C2ST as a metric in the appendix might provide additional comparison. While the focus of the paper is on calibration, it is possible to get a perfectly calibrated, but badly performing estimator.\n* There is no mention of hyperparameter optimisation (although a sensitivity analysis is given).\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer,\n\nThank you very much for taking the time to review our manuscript and for your comments. Below we would like to address the questions:\n\n- *Line 206: Does IS provide an improvement for NPE, even when it is not needed?*\n\n   Intuition suggests that IS should introduce noise and worsen the results. In a limited study (please see Figure 3 of the global response PDF) we did a sanity check on two problems (Weinberg - 1D posterior; Lotka Volterra - 2D posterior) where IS is not used for NPE. The results show that IS with 16 samples (submitted manuscript) gives almost the same outcomes as directly sampling from the approximate posterior with the reparameterization trick (Figure 3 in global response PDF). We expect the comparison to look very different for moderate and high-dimensional posteriors.\n\n\n- *Since the regulariser relies on importance sampling did the authors do any experimentation on how the approach scales with dimension of ? This was mentioned in future work, but even in the experiments in the paper the theta varies across experiments. What are the dimensions of the experiments included in this paper? And did the authors see that the higher dimensional ones performed worse with the new regulariser?*\n\n   We did experiments only with the problems mentioned in the submitted manuscript. The dimensions of posteriors in the submitted manuscript are as follows: SLCP - 2D; M/G/1 - 3D; Weinberg - 1D; Lotka Volterra - 2D; Spatial SIR - 2D; Gravitational Waves - 2D. All of the studied problems are of the same order of magnitude in terms of the dimensionality of parameters (which is not the case for the dimensionality of observations but this has limited impact on the regularizer), therefore we cannot draw conclusions about scaling. We hold that this should be verified in subsequent work."}, {"review_id": "mzMDAq7pud", "reviewer": "Reviewer_JSnb", "summary": "The paper presents a method to perform calibrated simulation-based inference. To do so, the paper employs the well-known coverage and proposes a way to differentiate through this term and to use it as a regularizer during training. The authors evaluate their method on benchmark tasks and conclude that it has good coverage and (sometimes) even outperforms existing methods in terms of log-likelihood.", "questions": "Why is the sorting based-computation preferrable over the direct computation? The authors say that `we need to backpropagate through F_N(\\alpha_k)`. Is this a problem in practice? \n\nL472: what do you mean by “proposed regularizer set to calibration objective”\n\nL505: what does “outstanding values in observations” mean?", "limitations": "The authors state limitations of their method, but some limitations have to be highlighted (see comments above, compute time for large batch sizes, poor performance for calibrated (not convservative) posteriors.", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "**Originality**: \nThe method is novel and the use of a differentiable sorting algorithm for differentiating through coverage is novel and useful.\n\n**Quality**: \nThe theoretical part of the paper is done rigorously and the paper provides additional empirical results for the impact of hyperparameters and computational cost.\n\n**Clarity**: \nThe figures are clear and support the messages of the paper.", "weaknesses": "**Quality**: \nI expect that the method is very expensive if the batch size is large because in this case, GPU will not help either. This should be clarified in the paper.\n\nI found it very interesting that CA1NPE outperforms NPE in terms of log-likelihood. What exactly are the methods that the authors use to prevent overfitting? Always training for 500 epochs is clearly not something anybody would do in practice (L497). Please use a proper implementation of NRE or NPE to draw comparisons.\n\nThe poor results of Appendix Fig 6 should be mentioned in the main paper and it should be highlighted that the method can be used to produce conservative posteriors, but that it is not suitable to produce calibrated posteriors.\n\nAll tasks in the paper are very low dimensional. I would expect the NPE version of this algorithm to scale to high-dimensional asks, but this would need emprical evidence. For NRE, I could imagine that Importance sampling requires exceedingly many samples in high-d parameter spaces. Please clarify and ideally add tasks with a more high-dimensional parameter space.\n\n\n**Clarity**: \n\nTwice (L159 and L183), the authors propose alternative formulations of their method. They never empirically investigate these formulations and also do not describe why they are less good. I would appreciate if the authors either add additional details on these methods or remove them entirely to avoid confusion.\n\nThe paper introduces **many** symbols which makes the paper very tedious to read. I would appreciate if the authors would redefine symbols in new sections to make the paper easier to follow. Also, some abbreviations do not really have to be defined (rarely used, e.g. KS, SAE, STE).\n\nMany papers are listed as `arxiv` although they got published. Please fix.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer,\n\nThank you very much for taking the time to review our manuscript and for your comments. Below we would like to address the questions:\n\n- *I expect that the method is very expensive if the batch size is large because in this case, GPU will not help either. This should be clarified in the paper.*\n\n   In L262 we say that the advantage of computation on GPU is due to the parallelization. Implicitly, this means that once <parameters dimensionality X batch size X number of samples $L$> exceeds the available memory, parallelization, and consequently GPU advantage is limited. We will explicitly mention this threat in the revised version of the manuscript.\n\n\n- *I found it very interesting that CA1NPE outperforms NPE in terms of log-likelihood. What exactly are the methods that the authors use to prevent overfitting? Always training for 500 epochs is clearly not something anybody would do in practice (L497). Please use a proper implementation of NRE or NPE to draw comparisons.*\n\n   We use the same training protocol for both NPE and CalNPE that was used in Hermans et al., 2022 - literally, the same implementation that is available online. Although we run 500 epochs of training in total, only the best model based on validation log-likelihood is kept as the final model. In addition, Gradient Norm Clipping is used. We consider the implementation from Hermans et al., 2022 to be proper (which is based on some standard packages used in the field of SBI). Moreover, both the regularized and non-regularized versions use the same implementation allowing us to draw conclusions.\n\n\n- *The poor results of Appendix Fig 6 should be mentioned in the main paper and it should be highlighted that the method can be used to produce conservative posteriors, but that it is not suitable to produce calibrated posteriors.*\n\n   We will add this information in the revised version of the manuscript.\n\n\n- *All tasks in the paper are very low dimensional. I would expect the NPE version of this algorithm to scale to high-dimensional asks, but this would need empirical evidence. For NRE, I could imagine that Importance sampling requires exceedingly many samples in high-d parameter spaces. Please clarify and ideally add tasks with a more high-dimensional parameter space.*\n\n   Following the nowadays state-of-the-art evaluation proposed in Hermans et al., 2022 already imposes a huge computational effort. We see the high-d parameter spaces as an important challenge and list it as a future work direction. However, with the available computational resources we are unable to provide empirical evidence in the short term.\n\n- *Many papers are listed as arxiv although they got published. Please fix.*\n\n   In the revised version of the manuscript, we will update all the references that got published since we submitted the manuscript.\n\n- *L472: what do you mean by “proposed regularizer set to calibration objective”*\n\n   Minimizing calibration error not the conservativeness error. We will clarify this information in the revised version of the manuscript.\n\n - *L505: what does “outstanding values in observations” mean?*\n\n   The experimental protocol of Hermans et al., 2022 (which we follow) did not include data standardization. The highest value encountered in observations for the largest simulation budget of M/G/1 was 4e+7 which is the “outstanding value” in L505. However, it was unlikely to happen, and therefore for the small simulation budgets, the extreme observations were less extreme leading to a stable learning process. We will replace “outstanding value” with “extreme value” to avoid confusion."}]

Author rebuttal:
Dear Reviewers,

Thank you very much for taking the time to review our manuscript and for your comments. Below we would like to address a question that has come up in several reviews:

*What are the advantages and disadvantages of using either the sorting-based computation or the direct computation method?*

First, we would like to underline that we made no claims of the superiority of the sorting-based computation (the one used in the experiments in the submitted manuscript) over the direct computation. Sorting-based computation happened to be the first one we evaluated empirically and since the method seemed to perform well we did not investigate the alternative approach in the light of limited computational resources.

Our recent empirical results - Figure 1 and Figure 2 in global response PDF - show that using direct computation instead of sorting-based computation while keeping all the remaining hyper-parameters untouched, leads to performance degradation both in terms of coverage and log-posterior density. $D_N$ is evaluated over 128 (same as the batch size) randomly sampled (in every iteration) levels on the (0,1) interval. We also observed optimization stability issues that reveal themselves in high variance over random initialization of the expected value of the approximate log posterior density of the nominal parameters for SpatialSIR in Figure 2 of the global response PDF.

We hypothesize that the poor performance of direct computation is due to the double use of a straight-through estimator through the indicator function - first in eg. (7) and then in eq. (8) - with the same backward relaxation. We see this issue as a direction for further investigation. Moreover, differentiable sorting can be seen as a straight-through estimator of a piecewise linear function suggesting an equivalence between sorting-based computation and direct computation but with different backward relaxations.


Meta-review: This paper considers simulation-based-inference (SBI) in a setting where inference might be expected to return an inexact posterior. A novel regularizer is developed that encourages SBI to find well-calibrated posteriors despite being inexact. Reviewers agreed this was a useful and novel contribution.
Decision: Accept (poster)