L-C2ST: Local Diagnostics for Posterior Approximations in Simulation-Based Inference
Authors: Julia Linhart
Alexandre Gramfort
Pedro L. C. Rodrigues
Venue: NeurIPS 2023 poster (2023)
Keywords: machine learning
calibration
simulation-based inference
neuroscience
normalizing flows
classifier two-sample tests
TL;DR: A validation method for conditional density estimators particularly tailored for the simulation-based inference setting, with theoretical guarantees and interpretable diagnostics.

Abstract: Many recent works in simulation-based inference  (SBI) rely on deep generative models to approximate complex, high-dimensional posterior distributions. However, evaluating whether or not these approximations can be trusted remains a challenge. Most approaches evaluate the posterior estimator only in expectation  over the observation space. This limits their interpretability and is not sufficient to identify for which observations the approximation can be trusted or should be improved. Building upon the well-known classifier two-sample test (C2ST), we introduce $\ell$-C2ST, a new method that allows for a local evaluation of the posterior estimator at any given observation. It offers theoretically grounded and easy to interpret -- e.g. graphical -- diagnostics, and unlike C2ST, does not require access to samples from the true posterior. In the case of normalizing flow-based posterior estimators, $\ell$-C2ST can be specialized to offer better statistical power, while being computationally more efficient. On standard SBI benchmarks, $\ell$-C2ST  provides comparable results to C2ST and outperforms alternative local approaches such as coverage tests based on highest predictive density (HPD). We further highlight the importance of local evaluation and the benefit of interpretability of $\ell$-C2ST on a challenging application from computational neuroscience.

L-C2ST: Local Diagnostics for Posterior Approximations in Simulation-Based Inference
Julia Linhart
Alexandre Gramfort"
Université Paris-Saclay, Inria, CEA Palaiseau 91120, France
Université Paris-Saclay, Inria, CEA
Palaiseau 91 120, France alexandre gramf
julia. www.arttitniia.....
fr
Pedro L. C. Rodrigues
Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK Grenoble 38000, France
roddoorrddiguestniii..f..
Abstract
Many recent works simulation-base inference (SBI) rely on deep generative models approximate complex, high- sigh-dimensional posterior distributions. How- ever, evaluating whether or not these approximations can be trusted remains a challenge. Most approaches evaluate the posterior estimator only in expectation over the observation space. This limits their interpretability and is not sufficient identify for which observations the approximation can be trusted or should be improved. Building upon the well-known classifier two-sample test (C2ST), we introduce E-C2ST, new method that allows for local evaluation of the posterior estimator at any given observation. It offers theoretically grounded and easy to interpret- - c.g. graphical diagnostics, and unlike C2ST, does not require access to samples from the true posterior. In the case of normalizing bow-based posterior estimators, E-C2ST can be specialized offer better statistical power, while being computationally more efficient. On standard SBI benchmarks, E-C2ST provides comparable results C2ST and outperforms alternative local approaches such as coverage tests based on highest predictive density (HPD). We further highlight the importance of local evaluation and the benefit interpretability of 2-C2ST on a challenging application from computational neuroscience.
Introduction
Expressive simulators the of modern experimental science, enabling the exploration of rare or challenging-to-meann events in complex systems across various fields such as population
genetics [43]. astrophysics [7], cosmology [32]. and neuroscience [28, 16, 1, 20]. These simulators implicitly encode the intractable likelihood function p(r 0) of the underlying mechanistic model, where 0 represents set of relevant parameters and Simulator (0) the corresponding realistic observation. The main objective is to infer the parameters associated with given observation using the simulator's posterior distribution p(0 .)) [4]. However, classical methods for sampling posterior distributions, such as MCMC [41] and variational inference [40], rely on the explicit evaluation of model-likelihood, which not possible when working with most modern simulators.
Simulation-based inference (SBI) [4] addresses this problem by estimating the posterior distribution on simulated data from the joint distribution. This be done after choosing prior distribution
A. Gramfort joined Meta and can be reached agramfort@mete
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

p(O) over the parameter space and using the identity p(0,x) p(z p)p(0). In light of recent developments in the literature on deep generative models, different families of algorithms have been proposed to approximate posterior distributions in SBI [4]. Certain works use normalizing flows [37] to directly learn the posterior density function (neural posterior estimation, NPE (18[) or aim for the likelihood (neural likelihood estimation, NLE [36]). Other approaches reframe the problem in terms of a classification task and aim for likelihood ratios (neural ratio estimation, NRE (22]).
However, appropriate validation remains a challenge for all these paradigms, and principled statistical approaches are still needed before SBI can become a trustworthy technology for experimental science.
This topic has been the goal of many recent studies. For instance, certain proposals aim at improving the posterior estimation by preventing over-confidence [23] or addressing model misspecification [14] to ensure conservative [8] and more robust posterior estimators [46, 27]. Another approach is the development of SBI benchmark [31] for comparing and validating different algorithms on many standard tasks. While various validation metrics exist, Lueckmann et al. [31] show that, overall, classifier two sample tests (C2ST) [30] are currently the most powerful and flexible approach. Based on standard methods for binary classification, they can scale to
well handle non-Euclidean data spaces [26, 34]. Typical use-cases include tests for statistical independence and the evaluation of sample quality for generative models [30]. Implicitly, C2ST is used in algorithms such as noise contrastive estimation [19] and generative adversarial networks [17], or to estimate wet-iiooo--t--vvidenn
ratios [22]. To be applied in SBI settings, however, C2ST requires access to samples from the true target posterior distribution, which renders useless in practice. Simulation- based calibration (SBC) [42] bypasses this issue by only requiring samples from the joint distribution. Implemented in standard packages of the field (sbi [44], Stan [2]). it has hecome the go-to validation method for SBI [31, 23] and has been further studied in recent works [33, 5, 15). Coverage tests based on the highest predictive density (HPD) as used in [23, 8], can be seen as a variant of SBC that are particularly well adapted to multivariate data distribution.
Nevertheless, a big limitation of current SBI diagnostics remains: they only evaluate the quality of the posterior approximation globally (in expectation) over the observation space and fail to give any insight of its local behavior. This hinders interpretability and can lead to false conclusions on the validity of the estimator [48, 29]. There have been attempts to make existing methods local, such as local-HPP) [48] or local-multi-PIT [29], but they depend on many hyper-parameters and computationally too expensive to be used in practice. In this work, we present (-C2ST, new local validation procedure based on C2ST that can be used to evaluate the quality of SBI posterior approximations for any given observation, without using any data from the target posterior distribution. (-C2ST comes with necessary, and sufficient, conditions for the local validity of multivariate posteriors and is particularly computationally efficient when applied validate NPE with normalizing flows, as often done in SBI literature [7, 16, 27, 46, 6, 45, 24]. Furthermore, E-C2ST offers graphical tools for analysing the inconsistency of posterior approximations, showing in which regions of the observation space the estimator should be improved and how to act upon, e.g. signs of positive / negative bias, signs of over / under dispersion, etc.
In what follows, we first introduce the SBI framework and review the basics of C2ST. Then, we detail the L-C2ST method and prove asymptotic theoretical guarantees. Finally, we report empirical results on two SBI benchmark examples to analyze the performance of E-C2ST and a non-trivial neuroscience use-case that showcases the need of local validation method.
2 Validating posterior approximations with classifiers
Consider a model with parameters E I and observations x E RRd obtained via a simulator. In what follows, we will always assume the typical simulation-based inference setting, meaning that the
likelihood function p( ()) of the model cannot be easily evaluated. Given a prior distribution p(0).
it is possible to generate samples from the joint pdf p(0, x)) as per:
O₂ p00) +
X, Simulator((.,, (e..) ~ p(x On)) = (On, Xn) ~ (((,x))
(1)
Let N, be a fixed simulation budget and (x,, I = Dtrain Deal with Dirain n Deal = V. The data from Dirain are used to train an amortized approximation q( x) ze p(0) x), e.g. via NPE [18], and those from Deat to diagnose its local consistency [48].
P..e. the approximation g( 2)) is close to p(0 2) on average for all values of x E Rd so we can quickly generate samples from the posterior for any choice of conditioning observation without redoing any training.
2

---

Definition (Local consistency), A conditional density estimator is said to be locally consistent at To with the true posterior density p if. and only if. the following null hypothesis holds:
Ho(xo) q(0 xo) p(0 Io), 40
R
(2)
We can reformulate Ho(x)) as a binary classification problem by partitioning the parameter space into two balanced classes: one for samples from the approximation (C = 0) and one for samples from the true posterior (C = 1), as in
θ (C 0) 7(0 xo) vs. θ (C 1) ~ p(0 Io),
for which the optimal Bayes classifier [21] fz. (0) argmax { dz. (0), (0)) with
d² PC(C 1 - 0:
PCC 0 0| O 0;: p( 2b
p(0 m0)
is a standard result [30, 26] to relate (2) with (3) as per
Hz(),) holds + d= (0) P(CC = 1 0: -===))= =
40
When the classes are non-separable, the optimal Bayes classifier will be unable to make a decision and we can assume that behaves as Bernoulli random variable [30].
2.1 Classifier Two-Sample Test (C2ST)
The original version of C2ST [30] uses (5) to define a test statistic for Ho( based on the accuracy of classifier fz. trained on dataset defined as
0 {(0 xo)
and
0% ~ p(0 mo) and D = (00 0)) U ((O%) 1)
(6)
C=0
C=1
The classifier accuracy is then empirically estimated over 2N, samples samples in each class) from held-out validation dataset D, generated in the same way as D;
2N.
(f.)
Σ (ff. =0) (0))
2N
Theorem (Local consistency and classification accuracy). If fz. is Bayes optimal' and N.- co, then (foc)f 1/2 is a necessary and sufficient condition for the local consistency of q at To-
See Appendix A. for proof. The intuition that, under the null hypothesis Ho(x)), is impossible for the optimal classifier to distinguish between the two data classes, and its accuracy will remain at chance-level [30]. In the context of SBI, C2ST has been used to benchmark a variety of different procedures on toy examples where the true posterior is known and can be sampled [31]. This is why we call this procedure an oracle C2ST, since uses information that is not available in practice.
Regression C2ST. Kim et al. [26] argues that the usual C2ST based on the classifier's accuracy may lack statistical power because of the "binarization" of the posterior class probabilities. They propose to instead use probabilistic classifiers (c.g. logistic regression) of the form
fa. (0) = (da, (0) > 2)) where (0) = P(C = 1 (;
(8)
and define the test statistics in terms of the predicted class probabilities dz, instead of the predicted class labels. The test statistic is then the mean squared distance between the estimated class posterior probability and one half:
N (MSR N Σ
1)²²
1 2
N. Σ n=1
1 1)²
d.. (0g)
2)
(O.)
dr.
-
(9)
N
Theorem 2 (Local consistency and regression). If fr. is Bayes optimal and N₂ - 0o, then (f 0 is necessary and sufficient condition for the local consistency of q at To.
See Appendix A.2 for a proof. The numerical illustrations Kim et al. [26] give empirical evidence that Regression C2ST has superior statistical power as compared to its accuracy- based counterpart, particularly for high-ding dige-dimensional data spaces. Furthermore, it offers tools for interpretation and visualization: evaluating the predicted class probabilities dz, (f) for any f E R informs the regions where the classifier is more (or less) confident about its choice [30, 26].
L.e it is the classifier with lowest possible classification error for the dataset D.
3

---

(-C2ST: Local Classifier Two-Sample Tests
The oracle C2ST framework is not applicable in practical SBI settings, since requires access to samples from the true posterior distribution to (1) train a classifier and (2) evaluate its performance discriminating data from 9 and p. This section presents a new method called local C2ST (e-C C2SC)T)) capable of evaluating the local consistency of posterior approximation requiring data only from the joint pdf p(e, 2) which can be easily sampled as per (I).
(1) Train the classifier. We define a modified version of the classification framework (3) with:
(0, (C 0) q(0 x)p(z)
vs. (0, (C 1) p00,z).
(10)
The optimal Bayes classifier is now f= x) = argmax {1 d" (0, x). d d0, x)) with
(0,,x)
p(0 x))
(0.0))) .x) p(0,x) q(O x)p(z> p x)) + ( x)) = d0)
(11)
where notice the direct relation with the Bayes classifier for (3), Therefore, using data
sampled as in (10), is possible to train a classifier f(0, x) and write Sa. (0) = [[0, xo) for each To-
See Algorithm I for details on the implementation of this procedure.
(2) Evaluate the classifier. Define a new test statistic that evaluates the MSE-statistic for classifier and its associated predicted probabilities d using data samples from only the class associated to the posterior approximation (C = 0):
1 Nu No
1) with 0% ~ q(O xo)
fASE (f. = Σ n=1
d(er, n: 2o)
(12)
Theorem 3 (Local consistency and single class evaluation). If is Bayes optimal and N₂. -> oo. then fMSE, v.) 0 necessary and sufficient condition for the local consistency of To- Proof. Let d be an estimator of IP(C =1 O) and f = I(d > 0.5). Suppose that f = f is Bayes optimal and let vo be a fixed observation. We have that
lim iMSE, _0, 1 f No0
d (
2) q( do)d0.
Because of the squared term in the integral and being a p.d.f., we have that d= (0) PCC 0; x.)
lim ÊMSE, (f*,
No0
This new statistical can thus be used to assess the local consistency of posterior approximation without using any sample from the true posterior distribution but only from the joint pdf. Furthermore, is amortized, so single classifier trained for (10) that can then be used for any choice of conditioning observation F., This is not the case in the usual oracle C2ST framework.
Optimal Bayes (statistics)
QDA-C2ST (power)
1.0
1.0
Figure 1: Results for the C2ST framework when p = (0, I2) and g == =O I212). Left panel por- trays the test statistics for the opti- mal Bayes classifier and the right panel shows the test's empirical power with QDA. Single-class accu- racy test (iAera) fails to detect when 7 q but IMSE, behaves correctly.
buss, (-05)
0.5
0.5
H6
0.0
0.5
1.5
d 1.0
1.5
6²
Figure illustrates the behavior of different test statistics to discriminate samples from two bivariate Normal distributions whose covariance mismatch is controlled by single scaling parameter a. Note that the optimal Bayes classifier for this setting can be obtained via quadratic discriminant analysis (QDA) [21]. The results clearly show that even though (MSE, exploits only half of the dataset (i.e. samples from class C' 0) it is capable of detecting when p and g are different (o * 1). The plot also

---

includes the results for a one-class test statistic based on accuracy values (fAcc,) which, as opposed to ÊMSE. has no guarantees for being a necessary and sufficient condition for local consistency. Not surprisingly, it fails reject the null hypothesis for various choices of
The assumptions of Theorem are never met in practice: datasets are finite and one rarely knows which type of classifier is optimal for a given problem. Therefore, the values of (MSF, in the null hypothesis (p = g) tend to fluctuate around one-half, and it is essential to determine threshold for deciding whether or not Ho) should be rejected. In &-C2ST, these threshold values are obtained via a permutation procedure [13] described in Algorithm 1. This yields Nx estimates of the test statistic under the null hypothesis and can be used to calculate p-values for any given a significance level as described in Algorithm 2. These estimates
also be used to form graphical summaries known as PP-plots, which display the empirical CDF of the probability predictions versus the nominal probability level. These plots show how the predicted class probability d(x.)) deviates from its theoretical value under the null hypothesis (i.e. one half) as well as (1 - a) confidence regions; see Algorithm B. available in the appendix for more details and Figure for an example.
Algorithm 1: L-C2ST - training the classifier data from the joint distribution
Input: posterior estimator q; calibration data Deal {{,,. Xx, } Ney; classifier f: number of samples Nu from the distribution under the null hypothesis Output: estimate d of the class probabilities; estimates (d1,.... dN,} under the null hypothesis ** Construct classificatio
training set
*/
for n 1, I..... Neat do 0 ~ q( X
2m (0,, X,): C2n 0 Sample from q() x)p(x)
W2n+1 = (On, Xn); C2n+1 = 1 /* Sample from v(0,x)
D {W,, CN
/+ Get estimate d of the class probabilities
Train the classifier f on D
d + Prrabability
/* Estimate under the null hypothesis via permutation procedure for h -...... Nu do
Randomly permute labels C, in D Train the classifier f on new D dn Sprobability
return d; (d...... ]N
Algorithm 2:
- - test statistics and p-values for any To
Input: Observation Fo; estimates and (di...., dN } obtained in Algorithm
Output: test statistic fMSE, (zo); p-value p(ro)
Generate No samples 6 ~ q(0) zo) with predicted probabilities /+ Compute test statistics
and (0g, 2o) */
fMSE (Ia) +- N, Σ. (d(Og, x.)) - 2)²
for h 1,...,, NH do
E (xo) + 1 N. Σn (dn (0q, xo) - =
/* Compute p-value
p(x.)) No Σx ((In INsSTo (zoo)
return ÉMSE, (Io), p(zx)
3.1 The case of normalizing flows
The t-C2ST framework can be further improved when the posterior approximation is conditional normalizing flow [37], which we denote qg. Given Gaussian base distribution u(>) (0, and a bijective transform (=(:z) with Jacobian Jr, (::xz) we have
(0 ()) u(2) det JI (2;z)1 0 To(2;)) m
(13)
5

---

In other words, normalizing flows (NF) are invertible neural networks that define map between latent space where data follows a Gaussian distribution and the parameter space containing complex posterior distributions. This allows for both efficient sampling and density evaluation:
(00 mm) = 6 = (2 ~ ( x)
(14)
1o [[01x)] x)) - u (0;x((;;)) (0;z)) det JT.- (0; (0:z)]
(15)
Our main observation is that the inverse transform T5 can also be used characterize the local consistency of the conditional normalizing flow in its latent space, yielding much simpler and computationally less expensive statistical test for posterior local consistency.
Theorem (Local consistency and normalizing flows). Given a posterior approximation qø based on normalizing flow, its local consistency at To be characterized as follows:
p(0 Øo) = (0 2o) = p(Tx (0; 20) xo) u(z). v0 E -
(16)
Proof. Let + p(0 so). Following (14), have that θ (o Ho) if, and only if, = zo) with Z ~ No Jm). Applying the inverse transformation of the flow gives us TV (0, TT (2; [o; = Z ~ NO, Im). which concludes the proof.
Based on Theorem 4 we propose a modified version of our statistical test named É-C2ST-NE. The new null hypothesis associated with the consistency the posterior approximation Os at 2o is
H (z.) p(TT x.) [ao = Im),
(17)
which leads to new binary classification framework
(Z, X) (C 0) ~ NNO0, (p)=
(C= (C=1)~<(T'(00a)000
(18)
Algorithm 3 describes how to sample data from each class and train classifier to discriminate them. The classifier then evaluated on N, samples Zn ~ N(O, Im) which are independent of Ta-
A remarkable feature of E-C2ST-NF is that calculating the test statistics under the null hypothesis considerably faster than for E-C2ST. fact, for each null trial h ., N₂ we use the dataset Deal only for recovering the samples X,, and then independently sample new data Z, NN, L..). As such, it is possible to pre-train the classifiers without relying on a permutation procedure (cf. Algorithm 4), and to re-use them to quickly compute the validation diagnostics for any choice of To or new posterior estimator ds of the given inference task. It is worth mentioning that this is not possible with the usual & E-C2ST, depends on q and it would require new simulations for each trial.
Algorithm 3: (-C2ST-NF training the classifier on the joint distribution
Input: NF posterior estimator qa: calibration data Deal = {On, X, ) classifier f Output: estimate d of the class probabilities /* Construct classification training set
Noa
for in ....., Neat do
Zn (O, mm); Z = T (0m; /* inverse NF-traneffrmatioo
Wen (Z, Xa); C2n
W2n+1 = 74 Xn); C2n+1 =1
D + {W,, Cn
2Naai 2N
/+ Get estimate d of the class probabilities
Train the classifier f on D
Sprabability return d

---

Algorithm 4: (-C2ST-NF-- -CCSS--NE-preoompucc the null distribution for any estimator
Input: calibration data Deal {O,, Xn )y classifier f; number of null samples Nu Output: estimates [di,..., dNz) of the class probabilities under the null for h in NN do
/* Construct classification training set forn = I,..... 2N
Sample Zn ~ N/ [m)
D 1-Z2n (,),,; X), 0) Nay U ((Z2n+1, 1) Nay
/* Get estimate d of the class probabilities
Train the classifier f on D
dn c Sprohailiity
return (di...., dNu
Experiments
All experiments were implemented with Python and the sbi package [44] combined with Py- Torch [38] and nflows [12] for neural posterior estimation 4, Classifiers on the C2ST framework use the MLPClassifies from seikit-lear [39] with the same parameters as those used in sbibm [31].
4.1 Two benchmark examples for SBI
We illustrate 2-C2ST on two examples: Two Moons and SLCP. These models have been widely used in previous works from SBI literature [18, 36] and are part of the SBI benchmark [31]. They both represent difficult inference tasks with locally structured multimodal true posteriors in, respectively, low (0 E R², x E I2) and high (0 E IR5, x R**) dimensions. See [31] for more details. To demonstrate the benefits of f f-C2ST-NF, all experiments
neural spline flows [11] trained under the amortized paradigm for neural posterior estimation (NPE) [37]. We use implementations from the sbibm package [31] to ensure uniform and consistent experimental setups. Samples from the true posterior distributions for both examples are obtained via MCMC and used to compare (-C2ST(NFF) to the oracle-C25) framework. We include results for local-HPP implemented using the code repository of the authors of [48] with default hyper-parameter and applied to HPD.
First, we evaluate the local consistency of the posterior estimators of each task over ten different observations No with varying Ntrain and fixed Noal = 10 The first column of Figure 2 displays the MSE statistics for the oracle and C-C2ST frameworks. As expected, we observe a decrease in test statistics as Nuaio increases: more training data usually means better approximations. For Two Moons the statistic of E-C2ST decreases at the same rate oracle-C2ST. with notably accurate results for E-C2ST-NF. However, in SLCP both E-C2ST statistics decrease much faster than the oracle. This is possibly due to the higher dimension of the observation space in the latter case, which impacts the pracedug-proeeuure of 2-C2ST on the joint distribution.
We proceed with an empirical analysis based on 50 random test runs for each validation method and computing their rejection-rates to the nominal significance level of a 0.05. Column 4 in Figure 2 confirms that the false positive rates (or type errors) for all tests are controlled at desired level. Column 2 of Figure 2 portrays the true positive rates (TPR), i.e. rejecting Ho(x,) when p + q. of each test as function of Nusin Both 2-C2ST strategies decrease with Nuran as in Column 1, with higher TPR for (-C2ST-NF both tasks. The oracle has maximum TPR and rejects the local consistency of all posterior estimates across all observations at least 90% of the time. Note that SLCP is designed to be a difficult model to estimate meaning that higher values of TPR are expected ( A in Column 1). In Two Moons, the decreasing rejection rate can be seen as normal the convergence of the posterior approximator (i - in Column 1): Nirsin increases, the task of differentiating the estimator from the true posterior becomes increasingly difficult.
reflects
We fix Nusin 10 (which yields inconsistent qe in both examples) and investigate in Column 3 of Figure how many calibration samples are needed to get a maximal TPR in each validation method.
*Code is available at attps://githul com/ htt//Jllaala. aLinhartnnnntc
'The average run- un-times for each validation method are provided Appendix C. SSLCP= = simple likelihood, complex posterior
7

---

SLCP
opee erur (FPR)
1.0
lest stsistie
(176
power (TPR)
power (TPR)
Two Moons
0.25
1.00
1.0
type empor (FPR)
1.0
1e1t statssle
power (TPP)
power (TPR)
orade CISI (heee Taod
PCCSST :NF (has.)
0.12
0.5
0.5
0.5
Aeaa/HPD
0.00
10
10'
100
10
10a
10
16
10
100
500
1000 2000 5000
10
500
1000 2000 5000 16
Muai 1N
Mon (a = 10*)
New (N = 10*)
Naa (Nus = 10*)
Figure 2: Results on two examples from the SBI benchmark: SLCP and Two Moons. We compare (-C2ST and E-C2ST-NF (dashed) to the oracle C2ST and local-HPD. Columns I and 2 display the test statistic and empirical power as function of Nurain- while Columns and show the empirical power and type error for varying Neal. The (-C2ST(-NF) statistics are comparable to the oracle, as their decreasing behaviour reflects the convergence of NPE to the true posterior for large training datascts. We also note that E-C2ST-NF is uniformly better than t-C2ST (i.e. higher power for all Ntrain and increases faster with Neat), and both reach maximum statistical power with smaller Neal than local-HPD. All Type I errors controlled different observations To (mean and std) and Columns 2-4 used additional 50 random test runs.
0.05. Experiments were performed over 10
SLCP is expected to be an easy classification task, since the posterior estimator is very far from the true posterior (large values of i in Column 1). We observe similar performance for (-C2ST-NF and E-C2ST, with slightly faster convergence for the latter. Both methods perform better than local-HPD, that never reaches maximum TPR. Two Moons represents a harder discrimination task, as is already pretty close to the reference posterior (see Column 1). Here, L-C2ST-NF attains maximum power at Neal 2000 and outperforms all other methods. Surprisingly. the regression- based oracle-C2ST performs comparably to local-HPP converging to TPR. 1 at Neat = 5000.
4.2 Jansen-Rit Neural Mass Model (JRNMM)
We increase the complexity of our examples and consider the well known Jansen & Rit neural mass model (JRNMM) [25]. This is neuroscience model which takes parameters O = (C, 0, 9) R as input and generates time series x E R 1024 with properties similar to brain signals obtained in neurophysiology. Each parameter has a physiologically meaningful interpretation, but they not relevant for the purposes of this section; the interested reader is referred to [3] for more details.
The approximation qs of the model's posterior distribution is conditioned masked autoregressive flow (MAF) [35] with 10 layers. We follow the same experimental setting from [3], with a uniform prior over aaaiooogically-releaaa parameter values and simulated dataset from the joint distribution including N Nirain = 50 000 training samples for the posterior estimator and Neal = 10 000 samples to compute the validation diagnostics. An evaluation set of size Non 10 000 is used for (-C2ST-NE.
We first investigate the global consistency of our approximation, which informs whether do is consistent (or not) on average with the model's true posterior distribution. We use standard tools for this task such as simulation-based calibration (SBC) [42] and coverage tests based on highest predictive density (HPD) [48]. The results are shown in the left panel of Figure 3. We observe that the empirical edf of the global HPD rank-statistic deviates from the identity function (black dashed line), indicating that the approximator presents signs of global inconsistency. We also note that the marginal SBC-ranks are unable to detect any inconsistencies in qo-
We use 8-C2ST-NF to study the local consistency of do on set of nine observations (4) defined as
( = (R(M with 0(*) (135, 220, 2000, gf and 969 E I-
(19)
The right panel of Figure 3 shows that the test statistics of E-C2ST-NF vary in a U-shape, attaining higher values as Io deviates from zero and at the borders of the prior. The plot is overlaid with the 95% confidence region, illustrating how much the test statistics deviate from the local null hypothesis.
Note that the uniform prior for g is defined in [-20, +20] when training de with NPE.

---

1.0
Figure 3: Results for global and lo- cal tests on JRNMM. Left: PP-plots for the marginal SBC and global HPD rank statistics. Right: Test statistics for 8-C2ST-NF observa- tions with varying Go- SBC fails to detect any inconsistency of qa. while HPD only provides a global assess- ment, unlike E-C2ST which locally explains the inconsistencies in qø-
SBC -0
0.10-
E-C2ST-NE (MSE.)
SBC 02
SBC- -
empiriall CDF
SBC 04
Best statistic
95% CR
0.5
0.05
HPD(O)
0.00
0.0
0.0
0.5
1.0
20
-5
10
15
20
rank
gain (ga)
Figure 4: Graphical diagnostics of 8-C2ST-NF for JRNMM. Top right panel displays the empirical CDF of the classifier (blue) overlaid with the theoretical CDF of the null hy- pothesis (step function at 0.5) and 95% confidence region of estimated classifiers under the null displayed in gray. The pairplot displays his- tograms of samples from qs within the prior region and dashed lines indicate values of 0g used gen- erate the conditioning observation Too The predicted probabilities are mapped on the colors of the bins in the histogram. Blue-green (resp. or- ange-red) regions indicate low (resp. high) predicted probabilities of the classifier. Yellow regions corre- spond to chance level, thus Uo Rs p.
Local PP Plot
10
empirical
400
200
0.0
0.0
0.5
1.0
0.8
3500
1000
0.5
20
0-
-2
0.2
100
250
200
400
10003500
20
20
C
We demonstrate the interpretabilityy of the results for E-C2ST-NF with focus on the behavior of qo when conditioned on an observation for which Jo 10. The local PP-plot in Figure 4 summarises the test result: the predicted class probabilities deviate from 0.5, outside of the 95%-CR, thus rejecting the null hypothesis of local consistency at 9. 10. The rest of Figure 4 displays ID and 2D histograms of samples @ ~ (0 xo) within the prior region, obtained by applying the learned T7 to samples Z (0, I.). The color of each histogram bin is mapped to the intensity of the corresponding predicted probability in É-C2ST-NF and informs the regions where the classifier is more (resp. less) confident about its choice of predicting class 0. This relates to regions in the parameter space where the posterior approximation has too much (resp. not enough mass) w.r.. to the true posterior: 9s > (resp. Q p). We observe that the ground-truth parameters indicating positive bias for y and a and negative bias for g in the 1D marginal. also shows that the posterior is over-dispersed in all 2D marginals. See Appendix D for results on all observations zo
often outside of the red regions,
5 Discussion
We have presented E-C2ST, an extension to the C2ST framework tailored for SBI settings which requires only samples from the joint distribution and is amortized along conditioning observations. Strikingly, empirical results show that, while L-C2ST does not have access to samples from the true posterior distribution, it is actually competitive with the macle-C2ST approach that does. This comes at the price training a binary classifier on potentially large dataset to ensure the correct calibration
*Specifically, we compute the average predicted probability of class 0 for data points Z ~ 14) corresponding to samples T (Z; ~ (0 2o) within each histogram bin.

---

of the predicted probabilities. Should this be not the case, some additional calibration step for the classifier can be considered [10].
Notably, E-C2ST allows for a local analysis of the consistency of posterior approximations and is more sensible, precise, and computationally efficient than its concurrent method, local-HPD. Appendix F.4 provides a detailed discussion of these statements, based on results obtained for additional benchmark examples. When exploiting properties of normalizing flows, f-C2ST can be further improved as demonstrated by encouraging results on difficult posterior estimation tasks (see Table 2 in Appendix F). We further analyze the benefits of this -NF version in Appendix F.3. e-C2ST provides necessary. and sufficient, conditions for posterior consistency. features that are not shared by other standard methods in the literature (e.g. SBC). When applied to widely used model from computational neuroscience, the local diagnostic proposed by E-C2ST offered interesting and useful insights on the failure modes of the SBI approach (e.g. poor estimates on the border of the prior intervals), hence demonstrating its potential practical relevance for works leveraging simulators for scientific discoveries.
Limitations and Perspectives
Training a classifier with finite data. The proposed validation method leverages classifiers to learn global and local data structures and shows great potential in diagnosing conditional density estimators. However, it's validity is only theoretically guaranteed by the optimality of the classifier when N,, + co. In practice, this can never perfectly be ensured. Figure in Appendix F.2 shows that depending on the dataset, E-C2ST can be more or less accurate w.r.t. the true C2ST. Therefore, one classifier and make
should always be concerned about false conclusions due to a far-from-optime sure that the classifier is "good enough" before using it as diagnostic tool, e.g. via cross- cross-validatie validation. Note, however, that the MSE test statistic for L-C2ST is defined by the predicted class probabilities and not the accuracy of the classifier, thus one should also check how well the classifier is calibrated.
Why Binary Classification? An alternative to binary classification would have been to train second posterior estimate in order to assess the consistency of the first one. Indeed, one could ask whether training a classifier is inherently easier than obtaining a good variational posterior, the response to which is non-trivial. Nevertheless, we believe that adding diversity into the validation pipeline with two different ML approaches might be preferable. Furthermore, building our method on the C2ST framework was mainly motivated by the popularity and robustness of binary classification: it is easy to understand and has a much richer and stable literature than deep generative models. As such, we believe that choosing a validation based on binary classifier has the potential of attracting the interest of scientists across various fields, rather than solely appealing to the probabilistic machine learning community.
Possible extensions and improvements. Future work could focus leveraging additional infor- mation of g while training the classifier as in [47]. For example, by using more samples from the posterior estimator q or its explicit likelihood function (which is accessible when g is normalizing flow). On a different note, split-sample conformal inference could be used to speed up the y-value calculations (avoiding the time-consuming permutation step in Algorithm 1).
In summary, our article shows that E-C2ST is thcoretically valid and works on several datasets, sometimes even outperforming local-HPD, which to our knowlegde is the only other existing local diagnostic. Despite facing some difficulties for certain examples (just like for other methods as well), an important feature of E-C2ST is that can directly leverage from improvements in binary classification to adapt and enhance it for any given dataset and task. This makes É-C2ST a competitive alternative other validation approches, with great potential of becoming the go-to validation diagnostic for SBI practitioners.
Acknowledgments
Julia Linhart is recipient of the Pierre-Aguilar Scholarship and thankful for the funding of the Capital Fund Management (CFM). Alexandre Gramfort was supported by the ANR BrAIN (ANR-20- CHIA0016) grant while in his role at Université Paris-Saclay, Inria.
10

---

References
[1] Sean R. Bittner, Agostina Palmigiano, Alex T. Piet, Chunyu A. Duan, Carlos D. Brody, Ken- neth D. Miller, and John Cunningham. Interrogating theoretical models of neural computation with emergent property inference. eLife, 10, 7 2021. ISSN 2050084X. doi: 10.7554/eLife.5622
[2] Bob Carpenter, Andrew Gelman, Matthew Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language. Journal of statistical software, 76(1), 2017.
[3] Pedro Luiz Coelho Rodrigues, Thomas Moreau, Gilles Louppe, and Alexandre Gramfort. HNPE: Leveraging Global Parameters for Neural Posterior Estimation. In NeurIPS 2021, Sydney (Online), Australia, December 2021. URL htt.h:t/p///aaa science/hall -a3te/3911100319116611
[4] Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences (PNAS), 117:30055-30063 2020. ISSN 0027-8424. doi: 10.1073/pnas.19 12779117.
[5] Niccolo Dalmasso, Rafael Izbicki, and Ann Lee. Confidence sets and hypothesis testing in likelihood-free inference setting. In Hal Daumé III and Aarti Singh, editors, Proceedings the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 2323-2334. PMLR, 13-18 18 Jul 2020.
[6] Maximilian Dax, Stephen Green, Jonathan Gair, Jakob H. Macke, Alessandra Buonanno, and Bernhard Schölkopf. Real-time gravitational wave science with neural posterior estimation. Phys. Rev. Lett., 127:241 103, Dec 2021. doi: 10.1 103/PhysRevLeet.... 127. 173.21103. URL https:
//link. aps. htt///oo1111011 1103/ 1wt.:/111///h/e.eeeeee whtp//hyeyteeetteeeeeeeee 127 241103.
[7] Maximilian Dax, Stephen R Green, Jonathan Gair, Michael Deistler, Bernhard Schölkopf, and Jakob H. Macke. Group equivariant neural posterior estimation. In International Con- ference Learning Representations, 2022. URL https: htap::::///oerreiew..eeee www.eeeee//tmmmmtmmmmm id= u6s8dSpor08.
[8] Arnaud Delaunoy, Jocri Hermans, François Rozet, Antoine Wehenkel, and Gilles Louppe. To- wards reliable simulation-base based inference with balanced neural ratio estimation. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Informa- fion Processing Systems, 2022. URL https ttttpps:/oooae /ttp////ppp/eee/////////// wwtprreeeeweeeeeefrreeeeeeeeeeee met//oormmiddo////ooooo 0722m 4XK.
[9] Biprateep Dey, David Zhao, Jeffrey A. Newman, Brett H. Andrews, Rafacl Izbicki, and Ann B. Lee. Calibrated predictive distributions via diagnostics for conditional coverage. 2022. doi: (0.48550/arxiv.2225. 14466. URL httts://aaxxxi http://aaii...rraa//////// /ttp://bbb/222255555522222
[10] Victor Dheur and Souhaib Ben Taieb. A large-scale study of probabilistic calibration in neural network regression. In Proceedings of the 40th International Conference on Machine Learning, 2023. To appear.
[11] Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios. Neural spline flows. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Bue, E. Fox, and R. Garnett, editors, Advances Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[12] Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios. nflows: normalizing flows in PyTorch. November 2020. doi: 10.528
[13] Bradley Efron and Trevor Hastie. Institute of mathematical statistics monographs: Computer age statistical inference: Algorithms, evidence, and data science series number Cambridge University Press, Cambridge, England, July 2016.
[14] David Frazier, Christian Robert, and Judith Rousseau. Model misspecification in approximate Bayesian computation: consequences and diagnostics. Journal of the Royal Statistical Society: Series B, 82(2):421-444 2019. doi: 10.1111/tssb. 12356.
[15] Andrew Gelman, Aki Vehtari, Daniel Simpson, Charles C. Margossian, Bob Carpenter, Yuling Yao, Lauren Kennedy, Jonah Gabry, Paul-Christian Bürkner, and Martin Modrák. Bayesian workflow, 2020.
11

---

[16] Pedro J. Gonçalves, Jan Matthis Lueckmann, Michael Deistler. Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William F. Podlaski, Sara A. Haddad, Tim P. Vogels, David S. Greenberg, and Jakob H. Macke. Training deep neural density estimators to identify mechanistic models of neural dynamics. elife, 9:1-46, 2020. ISSN 2050084X. doi: 10.7554/ELJPFE.56261
[17] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial networks. Communications of the ACM, 63:139-111 134,144, 2014. ISSN 15577317. doi: 10.1 10.1145/342262..
[18] David Greenberg, Marcel Nonnenmacher, and Jakob Macke. Automatic posterior transformation for likelihood-free inference. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97, pages 2404-2414, PMLR, 09-15 Jun 2019.
[19] Michael U. Gutmann and Aapo Hyvärinen. Noise-contrastive estimation of unnormalized statistical models, with applications to natural image statistics. Journal of Machine Learning Re- search, 13(11):077-361, 2012. URL http: /mlrr org/ wttpooatrs//poeoeoeoooooe.seeeet pappps//ppss gutt/nnt 1a.. html.
[20] Meysam Hashemi, Anirudh N. Vattikonda, Jayant Jha, Viktor Sip, Marmaduke M. Woodman, Fabrice Bartolomei, and Viktor K. Jirsa. Simulation- -based inference for whole-braid network modeling of epilepsy using deep neural density estimators. medRxiv, page 0t..2.06.02.2227555 6 2022. doi: 10. 10110002222202 000.022..6.02..22000 02.277666.. URL https:/ htapps://wwwww http://wwwwwwx/i//////////////pt medrxiv org/ org/cootent/... content/
10. 1...101/2022.00 06. 02. 22275860v1.
[21] Trevor Hastie, Robert Tibshirani, and J H Friedman. The elements of statistical learning. Springer series in statistics. Springer, New York, NY, 2 edition, December 2009.
[22] Joeri Hermans, Volodimir Begy, and Gilles Louppe. Likelihood-free MCMC with amortized approximate ratio estimators. Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 4239-4248. PMLR, 13-18 Jul 2020. doi: 10.48550/arxiv. 04057.
[23] Joeri Hermans, Arnaud Delaunoy, François Rozet, Antoine Wehenkel, Volodimir Begy, and Gilles Louppe. A crisis simulation-based inference? beware, your posterior approximations can be unfaithful. Transactions Machine Learning Research, 2022. ISSN 2835-8856. URL https: h/tp::://////eviie/e// net wwt/fformmeee (it.orrmiidLLHHHHtttaata
[24] Maëliss Jallais, Pedro L. C. Rodrigues, Alexandre Gramfort, and Demian Wassermann. Inverting brain grey matter models with likelihood-fres inference: a tool for trustable eytoarchitecture measurements, Machine Learning for Biomedical Imaging. 1:1-28, 2022. ISSN 2766-905X.
URL htttps://mmmllii journal. org/ 2::2:002:0000 010.
[25] Ben H. Jansen and Vincent G. Rit. Fleetroenceepallogam and visual evoked potential generation in mathematical model of coupled cortical columns. Biological Cybernetics 1995 73:4, 73: 357-366, 9 1995. ISSN 1432-0770. doi: 10. 10.000/0000119977..
[26] Ilmun Kim, Ann B. Lee, and Jing Lei. Global and local two-sample tests via regression. Electronic Journal of Statistics, 33:5253-5305 12 2018. ISSN 19357524. doi: 10.48550/arxiv, 1812.08927.
[27] Pablo Lemos, Miles Cranmer, Muntazir Abidi, ChangHoon Hahn, Michael Eickenberg, Elena Massara, David Yallup, and Shirley Ho. Robust simulation-based inference in cosmology with bayesian neural networks. Machine Learning: Science and Technology, 4(1):01LT0,, feb 2023. doi: 0008//2622-2153acbb533 URL https: (wa.xxddi..oog///// 1088/26322 2153/ acbb53.
[28] Paula Sanz Leon, Stuart A. Knock, Marmaduke Woodman, Lia Domide, Jochen Mersmann, Anthony R. Mcintosh, and Viktor Jirsa. The virtual brain: simulator of primate brain network dynamies. Frontiers in neuroinformatics, 7, 2013. ISSN 1662-5196. doi: 10.3389/FNINF. 2013.00011 URL https: /twp////beeeeddee/////// ncbi nlm. nih gtt.....////881199/////////
12

---

[29] Julia Linhart, Alexandre Gramfort, and Pedro L. C. Rodrigues. Validation diagnostics for SBI algorithms based on normalizing flows, 2022.
[30] David Lopez-Paz and Maxime Oquab. Revisiting classifier two-sample tests. 5th International Conference on Learning Representations, ICLR 2017, 10 2016. doi: (0.48550/axii.11100.6666 URL https: http://xxivivio// http//ass1111...55555 06545v4.
[31]
Jan-Matthis Lueckmann, Jan Boelts, David S Greenberg, Pedro J Gonçalves, and Jakob H Macke. Benchmarking simulation-based inference. Proceedings of The 24th International Conference on Artificial Intelligence and Statistics (PMLR), 130:343-355 4 2021. doi: 10.
48550/anxiv2201..4455
[32] Lucas Makinen, Tom Charnock, Justin Alsing, and Benjamin D. Wandelt. Lossless, scalable implicit likelihood inference for cosmological fields. Journal of Casmology and Astroparticle Physics, 2021, 2021. doi: (0.1088/1475-7516/0000
[33] Martin Modrák, Angie H. Moon, Shinyoung Kim, Paul Bürkner, Niko Huurre, Kateñina Faltejsková, Andrew Gelman, and Aki Vehtari. Simulation-based calibration checking for bayesian computation: The choice of test quantities shapes sensitivity. 11 2022. doi:
10.488550/rxxv.2211.. .02383. URL https: /ttp//xxvvvrrg....a attpaasss.....//////////aaaaaaaa 22t.///b///...2388311
[34] Teodora Pandeva, Tim Bakker, Christian A. Naesseth, and Patrick Forré. E-valuating classifier two-sample tests, 2022.
[35] George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation. Advances in Neural Information Processing Systems (NeurIPS), pages 2339-2348, 12 2017. ISSN 10495258. doi: 10.48550/arxiv1170555 170.000077.000.705
[36] George Papamakarios, David Sterratt, and lain Murray. Sequential neural likelihood: Fast likelihood-fii frree inference with autoregressive flows. 89:837-848, 16-18 Apr 2019.
[37] George Papamakarios, Eric Nalisnick. Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of 111500axxiv.112..22222
Machine Learning Research, 22:1 1-1,,,, 2021. ISSN 15337928. doi:
[38] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An Imperative Style, High-Performance Deep Learning Library. Advances in Neural Information Processing Systems (NeurIPS), page 12, Vancouver, BC, Canada, 2019.
[39] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, Prettenhofer. R. Weiss, V. Dubourg. J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-2830 2011.
[40] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In Francis Bach and David Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 1530-1538, Lille, France, 07-09 Jul 2015. PMLR. httpttps: /ttppeed////ccedingg.. mlr, press vv7 (007/rezen rreeenezze... html.
[41] Christian Robert and George Casella. Monte Carlo statistical methods. Statistics. Springer, New York, NY, edition, July 2005.
Springer Texts in
[42] Sean Talts, Michael Betancourt, Daniel Simpson, Aki Vehtari, and Andrew Gelman. Validating bayesian inference algorithms with simulation-based calibration. 2018. doi: 10.48550/arxiv. 1804.06788.
[43] Simon Tavaré, David J. Balding, R. C. Griffiths, and Peter Donnelly. Inferring coalescence times from DNA sequence data. Genetics, 145:505-518, 1997. ISSN 00166731. doi: (L.109//GENETICS//45522.00
13

---

[44] Alvaro Sejero-Cantern Jan Boelts, Michael Deistler, Jan- Man-Maathis Lueckmann, Conor Durkan, Pedro J. Gonçalves, David S. Greenberg, and Jakob Macke. sbi: A toolkit for simulation- based inference. Journal of Open Source Software, 5(52):2505, 2020. doi: 10.21 105/joss.02500
[45] Vasist, Malavika, Rozet, François, Absil, Olivier, Mollière, Paul, Nasedkin, Evert, and Louppe, Gilles. Neural posterior estimation for exoplanetary atmospheric retrieval. A&A, 672:A147, 2023. doi: 202245263.
URL httpp:://ool. org/ 1051/0004- /0061//000061////////
[46] Daniel Ward, Patrick Cannon, Mark Beaumont, Matteo Fasiolo, and Sebastian M Schmon. Robust neural posterior estimation and statistical model criticism. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022. URL Utt.attpp://operrrriee.. http:////opprrreeeeeeeefftonnnne htt.://ww..eeefoommddddtafaa Miteiii27dMBBaaa
[47] Yuling Yao and Justin Domke. Discriminative calibration. 2023.
[48] David Zhao, Niccolò Dalmasso, Rafael Izbicki, and Ann Lee. Diagnostics for conditional density models and bayesian inference algorithms. In Cassio de Campos and Marloes H. Maathuis, editors, Proceedings of the Phirty-Seventy Conference on Uncertainty Artificial Intelligence, volume 161 of Proceedings of Machine Learning Research, pages 1830-1840. PMLR, 27-30 Jul 2021. URL https: http:///e/opcseeeeeeeeee mlr. wwl.sprees///eeme. /vtsss/vv11////// v161/zhao21bbbbb zhb html.
14

---

Appendices
Proofs
16
A.1 Proof of Theorem
16
A.2 Proof of Theorem 2
17
A.3 Proof of Theorem 3
17
Algorithms
19
B.1 Generating PP-plots with E-C2ST
19
C Run-times for each validation procedure
20
D Graphical diagnostics (4-C2ST-NF) for the JRNMM posterior estimator
21
E On the cross-entropy loss for 8-C2ST
23
F Additional Experiments several benchmark examples
24
F.1 Scalability to high dimensions
24
F2 Accuracy of (-C2ST(-NF) w.r.t. the true C2ST
25
F3 Benefit of the -NF version
27
F44 Advantages of (-C2ST w.r.t. local-HPI)
27
15

---

A Proofs
In what follows, assume that it is sufficient for the null hypothesis Ho to hold any set C Th of strictly positive measure, rather than requiring it to hold for all points E R".". In practice, it generally has no practical implications if the posterior estimator inconsistent with the true posterior (q + p) on a set of measure zero, since those sets don't I have any real statistical significance.
For the proof' of Theorems I and 2, we will consider classifier fx. defined for fixed observation To E Rd on Sz, (0 E 18mh q(( xo) p(0 To) > 0}.
A.1 Proof of Theorem
Theorem (Local consistency and classification accuracy). If fz, is Bayes optimal and N.- -> oc, then Ace (fx.) 1/2 is necessary and sufficient condition for the local consistency of q at Zo.
Proof. As Sz, contains all data points 6% ~ q(0 xo) (C,, = 0) and Og, ~ p(0 xo) (C, = 1), the 2N. 1 is well defined (7) and empirical accuracy (f..) over the validation set D. = [(O,,. C,
(Au) - Aco(S2.) = (0) ( (0") No
(P(Jz.
= 0) + P(1z. (0P) = 1))
Let's show that if fu. is Bayes optimal, then Ho(xx) holds = Acc(fz.)= =
(+): Suppose Ho (z.) holds. The optimality of fz, implies that (09) and fr,, (0P) are Bernoulli random variables ((}) ) (see interpretation of equation (5)). and so Acc(Sm..) (1 2 2"
(+): Let's proceed by showing the contraposition: if Ho(Io) does not hold, then Acc(Ja.) + 1
Suppose that Ho (z)) does not hold, there exists a set C = {o E IR" p(0 xo) + q(0 xo) of strictly positive measure (w.r.t. p or g, which ever non zero on that set). We can decompose into the direct sum of A 0 C, x=) p(0 20)] and B 0 E q(0 x.) p(0 Io)
Either A or B is necessarily of strictly positive measure. Let's say A (see Lemma for the symmetric case).
A is exactly the region where Prob(C = 1 1 = 0) > 1/z and thus where (0) = 1; B is the region where fz. (0) O. We therefore get that:
1
(/
Acc(fzo)
q(0 x2)d0
p(0 2x)
2 100 (0)=0
1xa10)=1
1
(
L q(0 zz)d0 + p(0) [2)]d0)
2
1 (1 / p(0 I.)) (0 xx)dd)) (because / q=
2
But V0 A, 0 S q(0 ()) < p(0 I) (and A is of strictly positive measure), so the integral in the last equality is strictly positive and we get Aco(fz) > 2. which concludes the proof.
Let and g be two probability density functions defined on space S. If there exists set A C of strictly positive measure such that V0 E A, g(0) < p(0), then there exists a set B S of strictly positive measure such that v0' E B. x(0') > p(0')-
Lemma 1.
Proof. We know that SSP p Is q or equivalently, using AUB (g u
/ (8)dd + / q())do = -f p(0)do + / p())do =
By grouping the integrals over A on one side and the ones over B on the other, we get:
[ (())d0 / p(0)d0 p(0)do - [ (0))00
which non-negative because A assumed to be of strictly positive measure and q - ppp- > everywhere in A.
16

---

The integral of p over {p = q) is zero, which implies that
[ Jecp (p(0) ((0)) d0 = f
BB45PP ((00) 0
9(0)
meaning that the region {0 E E < {()) of strictly positive measure, which concludes the proof.
A.2 Proof of Theorem 2
Theorem 2 (Local consistency and regression). fa. is Bayes optimal and No. + oo, then (f.)(fz.)) = is a necessary and sufficient condition for the local consistency of q at 2o.
Proof. Let dx, be an estimator of IP region contains all the data points 0% q(O zo) (C, 0) and Og p(0 zo) (Cn 1), the mean squared error (f.. over the dataset D {(O,, Cn) is well defined (9) and
O æo) defined on Sa, such that fz. = -e, As this
(So)(S2.) N.
(ss)( =
d2. (0) 2 q(0) xO) + p(0 vo) m.)))a
1
This integral zero if, and only if (dz. (0) ) 0 for every 0 Sa.. (all terms are non-negative and q(0) | x.) + p(0 x.) > 0)", which is equivalent d (0) 1/2 for every E Sr., Assuming
Jz, = f2. is optimal, we have that dz, (0) d= (0) P(C = 0; v.) and we conclude the proof using the result from equation (5) (and knowing that outside of S... p = q = 0):
(MSE((Z.)
(0) P(C 11
v0 Sr. 2
Ho(xo)
holds
A.3 Proof of Theorem 3
Theorem 3 (Local consistency and single class evaluation). f is Bayes optimal and N₂ + co, then
bMSE, (f, xo) is necessary and sufficient condition for the local consistency of q at Vo-
Proof. Let d be an estimator of P(C = 6 X) and f = - & the associated classifier, both defined on S (w- (0,z) E R" Rd q0(0,z)) + p(0, x) > 0) Suppose that f f is optimal and let To Rdd be a fixed observation. As explained in section 3, we have that
d' d' 2.) P(C 1 0: x.) d= (0) and (8, ft, (0), VO E Sz,
Consider the support Sy..r. (0 E RT 0g ~ q( x.) from our single-class validation set Dro = (0)) N test statistic MMEE, is well defined (12) and
q(0 xo) > 0} C Sy, containing all data points
Noo Therefore our single-class
MSE,, (f)) Noo+t0 IMSE, (x,) = (1/ d (0)
1
<(
x.)d0
2
With the same arguments as the proof Theorem A.2, we get that
1
IMSE, (f* x.)
(0) PPC 11 2¹
A0 Sq.
Ha (x.)) holds
where the second equivalence is true because Sq.a. = Sz, for p(. xo) = q(. 2o). Therefore, tMsE, (f*, =.) 0 is necessary and sufficient condition for Ho (x.).
N.B. Single-class accuracy (Acco) does not provide sufficient condition for Ho (x.).
"ooe that this integral can also be zero if Sz. is of measure zero. But as mentioned at the beginning of this appendix, this has generally no practical implications.
17

---

Proof. Following the proof of Theorem in A.1, we have that
1 2 > (0)) = (0))) = =1
Ho(
holds - Acc(Sz.)
This means that Acco Tn)) = P(f (09) 0) can only be sufficient condition for Ho P(f2, (09) = 0) = PP([2, (OP) = 1), which is not generally true. In conclusion, evaluating
Acco(f" no) only provides a necessary condition for the local null hypothesis (see (=) in A.
18

---

B Algorithms
B.1 Generating PP-plots with (-C2ST
Algorithm 5: (-C2ST local PP-plots for any To
Input: evaluation data set Deo: an observation Zo: estimate d of the class probabilities; estimates {di.... dxx } under the null; grid & of PP-Ievels (0,1); significance level
Output: empirical CDF-values {F(I; 2o) hee of predicted class-0 probabilities;
(1 a>-coniideecc bands {Lo( Ui(xo) heg
Predict class-0 probabilities {du(u, zo) = d(e, (o),:ooeeee hEDD. / d is an estimate of
class for in 9 do
/* Compute empirical CDFs at
Fl): +- N.o Σ Idoooooo)s for h = N4 do (a + 1
Noo [l
/* Compute confidence bands at 1 Lı(z)), Ui(zo) + q
*/
No
/* quantiles
return (F(1 heg, (x.)). (x,0))ee
19

---

C Run-times for each validation procedure
To complete the benchmark experiments from Section 4.1, we analyze the run-times of each validation method to compute the test statistic 10 for an NPE of the SLCP-task, obtained for different values of Nirain. Results are displayed Table and computed on average over all given observations To and for Valuue that ensure high test power. We here focus solely on the SLCP-task, as the higher dimensional observation space allows to illustrate the differences between local methods (trained on the joint data space, (0,z) E R and the oracle (trained on the parameter space only, 0 E IR.).
As expected, the computation time increases with the sample size Neal (left vs. right part of Table 1). While the oracle C2ST has close to constant run-time for fixed Neat. local methods become faster with increasing Nirain. We observe that (-C2ST(-NF) has comparable run-times the oracle C2ST: the amortization cost is negligible, in particular for difficult tasks involving "good" posterior estimators (high Nirain). However, this is not the case for local-Hie Haa- which the most expensive method. Indeed, it involves (1) the costly computation of the HPD statistics on the joint and (2) the training of several (default is n 11) classifiers, both of which increase with the sample size and dimension of the data space.
Neal 5 000
Neal
10
000
Nurain
102
103
104
103
10a
10a
104
10
oracle C2ST
5.47
4.52
5.95
7.56
16.36
18.03
23.3
15.33
E-C2ST
5.92
5.09
1.78
1.81
27.62
34.06
17.9
3.65
1-C2ST-NF
6,98
6.84
6.38
1.72
43.99
25.01
18.56
7.62
local-HPD
282.19
282.85
279.38
282.5
956.91
938.21
682.92
530.45
Table I: Run-time (in seconds) to compute the test statistic for the SLCP task (mean over observations). C2ST has close to constant run-time for fixed Neal- Local methods become faster with increasing
Nuaain and E-C2ST(-NF) stays comparable to the oracle C2ST, even for Neal 10 000. While the amortization cost of & 8-C2ST is not an issue, Local-HPD always at least 30 times slower.
Note that in order compute to compute p p-values, we need to compute the test statistic NN times under the null hypothesis. The number of classifiers we need to train depends on how many we need to compute the test statistic (Nvx for É-C2ST vs. NN n. for local-HPD). In summary, if E-C2ST is more efficient in computing single test statistic, it will also be more efficient compute Nx test statistics.
20

---

D Graphical diagnostics (4-C2ST-NF) for the JRNMM posterior estimator
The following figures present additional results for the interpretability of f-C2ST applied to the JRNMM neural posterior estimator. They complete Figure 4. Results are shown for all given observations associated ground -otuhh gain values Jo -20, 15,, 10, -5, 0, 5, 10, 15, 20.
In each Figure, the top right panel displays the empirical CDF of the classifier (blue) overlaid with the theoretical CDF of the null hypothesis (step function 0.5) and 95% confidence region of estimated classifiers under the null displayed in gray. The pairplot displays histograms of samples from the posterior estimator 4e within the prior region and dashed lines indicate values of 0, used generate the conditioning observation To. The predicted probabilities are mapped on the colors of the bins in the histogram. Blue-green (resp. orange-red) regions indicate low (resp. high) predicted probabilities of the classifier. Yellow regions correspond to chance level, thus 23
Local PP-Plot
Local PP-Plot
1.0
1.0
2+
emyorical CDF
copimiaa
0.5
0.5
400
400
200
200
0.0
0.0
0.0
0.5
1.0
0.0
0.5
1.0
0.8
3500
3500
1000
1000
0.5
0.5
20
20
-20
20
0.2
100
250
200
400
1000
-20
100
250
200
400
1000
-20
0
20
C
07
8
Local PP-Plot
Local PP-Plot
10
empiriial CDF
copinial
0.5
0.5
400
400
200
200
0.0
0.0
1.0
0.8
0.8
3500
1000
1000
0.5
0.5
20
20
0-
-20
0.2
20
100
250
200
400
1000
20
100
250
200 400
1500
-20
20
C
21

---

Local PP-Plot
Local PP-Plot
1.0
1.0
80-33
coptrical CDT
empiriial CDD
0.5
0.5
400
400
200
200
0.0
0.0
02.5
1.0
0.0
0.5
1,0
0.8
3500
3500
1000
0.5
0.5
20
20
-0
-20
100
250
200
400
3500
20
0
20
100
250
200
400
10003300
-20
20
C
6
Local PP-Plot
Local PP-Plot
1.0
1.0
80
ooriinnaa
0.5
0.5
400
200
200
0.0
0.0
0.5
1.0
0.0
0.5
1.0
3500
3500
1000
1000
0.5
0.5
20
20
0-
-20
0.2
-20
0.2
100
250
200
400
1000
3500
-20
0
20
100
250
200 400
1000 3500
20
0
20
#
C
a
1.0
ampinial CDF
400
200
0.0
0.5
1.0
0.8
3500
1000
0.5
%6
0-
-20
0.2
100
250
200 400
15003500
20
20
#
22

---

E On the cross-entropy loss for 2-C2ST
This section aims at facilitating the understanding of E-C2ST by proving the result of equation (11).
As detailed section 3, the first step in &-C2ST consists in training a classifier to distinguish between Neal data points (+n, Xx. and (0% Xx.) from the joint distributions p(0, 2)) and q(0, x) respectively. Here, the same conditioning observations {_xn Waa are used to construct the data samples for each
class (see Algorithm 1). We show that this does not affect the objective function and convergence of the classifier.
The theoretical cross-entropy loss function to distinguish between data () from class C' = 1 and class C = 0 is defined by
(ce(d) == - (eex (log (d(o, ( Dog (1 d( X))]
(20)
Note that the we have equal marginals X =C p(z) = (C 0)." This allows us to take the expectation over X and approximate (20) via Monte Carlo for only one set of
conditioning observations {X, J and with data points 6, and O respectively associated to class C = 1 and C 0 for given X,:
(cE(d) Ex
- Log
(log (d(O, - -
[log d(o,
(21)
Newl
2
log (d(On, X,)) + log (1 - d(0,,
2N cal
n=1
Therefore, we can train classifier to minimize (20) using data from the joint distributions with same conditioning observations. The obtained estimate d arg min(Icr(d)) of the class probabilities p(9.z)
defined for every (0, x) E R Rd by d(0, x) p(0,=)+q(0,z)) As p(z) q(z), we recover the class probabilities of the optimal Bayes classifier the conditional data space for any given Rd
p(0 x)) d' (0, p(0 x) + {()) x.))
(0)
(22)
For an example, see works related to neural ratio estimation (NRE) [22, 8]: these algorithms implicitly use a classifier to distinguish between the joint and marginal distributions p(0,x) x)) and p(0)p(x). Like
in our case, both classes are modeled using the same observations Xn, obtained via the simulator.
"Teee joint distributions p(0,2) 2²) and q(0, 2) are both modeled using samples X ~ p(x () obtained from the prior è p(0), which implies that the marginals p(z) and q(x) are both defined by J p( ()p(O)d0
23

---

Additional Experiments on several benchmark examples
The results on the two benchmark examples in Figure 2 give first intuitions about the validity and behaviour of (-C2ST(-NF), our proposal, w.r.t. the oracle C2ST and the alternative local-HPD methods. In this section, Figure extends Figure with results on additional benchmark examples and more detailed results on the correlation between the statistics of é (-C2ST(-NF) and the oracle C2ST for different conditioning observations are shown in Figure 6 and Table 3. We refer the reader Table 2 for description of all benchmark examples (data dimensionality, posterior structure, challenges) and a summary of the main results comparing local methods.
While investigating the scalabilitty of the algorithms high dimensional data spaces, these addi- tional experiments help to further analyze how well (-C2ST(-NF) captures the true C2ST, when outperforms local-HPD, and better understand the benefit of the -NF version. These points further detailed in the following subsections. The goal is to create first guideline for when our proposal can and should be used, while raising awareness of its limitiations (i.e. when it should not be used or least be adpated for improved performance).
Dimension (0,z)
Posterior structure
Challenge
Better local method
SLCP
low
4 symmetrical modes
complex posterior
8-C2ST-NF
(5, 8)
Two Moons
low
bi-modal, erescent shape
global and local structure
local-HPD (-C2ST-NF
(2,
Caussian Mixture
low (2,2)
2D Gaussian
large vs. small s.t.d. in GMM
all similar
Gaussian Linear Uniform
medium (10, 10)
multivariate Gaussian
dimensionality scaling
local-HPD / C-C2ST
Bernoulli GLM
medium (10, 10)
unimodal, concave
dimensionality scaling
(-C2ST(-NF)
Bernoulli GLM Raw
high (10, 100)
unimodal, concave
raw observations (no summary stats)
& C-CST-NF
Table 2: Description of benchmark examples and summary of main results for local methods.
F.1 Scalability high dimensions
First all, note that the specific case SBI, the dimension the parameter space is typically of order 109 to 10a and m z 10a is already often considered as high dimensional. The observation space, however, can of higher dimension (e.g. d As 10 for time-series), but summary statistics are often used to reduce the dimension of the observations to the order of d R 102. In Section 4 we analyze the results obtained for rather low-dimensional benchmark examples (Two-Moons.ng SLCP: 5, d 8). As extension of Figure 2, Figure includes additional benchmark examples with low and medium dimensionality: Gaussian Mixture (m = 2, = 2) Gaussian Linear Uniform (m = 10, d 10) and Bernoulli GLM (m 10, d 10). Furthermore, the Bernoulli GLM Raw task allows us to analyze how our method scales to high dimensional observation spaces only (without parameter-span / task variability): it considers
2, =
observation data with d = 100, as opposed to sufficient summary statistics in the Bernoulli GLM task.
Column 3 in Figure shows that 2-C2ST requires more data to converge to the oracle C2ST (at maximum power TPR 1) as the data dimensionality increases: Neal 2000 for Two Moons and SLCP, but Nea1 & 5000 for the Bernouilli GLM task. Note that the Gaussian Mixture and

---

Gaussian Mixture 1.0
tppo crper (PPR)
1.0
les statssic
1196
power (TPR)
powerr (TPR)
oolondi 015
oradle CIST Oua Red
0.12
0.5
c-C157
4000
0.01
0.0
0.0
Gaussian Linear U Uniform
test statsslic
0.25
()Hb
power (TPR)
1.00
power (TPR)
1.0
type erter (FPP)
1.0
eelovd: 015
0.12
0.5
0.5
Bemoulli
1.0
type erer (FPR)
tox statnsic
(TPR)
power (TPR)
o levell 0.0
0.12
0.5
Bemoulli GLM Raw
1.0
1.0
sype errer (EPR)
1.0
1/
(TPR)
powee (TPP)
e-loves
0.12
0.5
3.5
10
10
10
10
10
10
10
10
100
500
1000
2000 5000
10
100
500 1000 2000 5000 10
(Na. 10
(N 10")
Nai (Noan =10")
Noc (N = 10")
Figure 5: Results additional sbibn bm benchmark examples using the same experimental setup for Figure in Section 4 (50 test runs, mean / std over 10 different reference observations).
Gaussian Linear Uniform tasks were not included in this analysis, as here, the difficulty of the classification task has more impact on statistical power than the data dimensionality 12,
Interestingly, we observe in the Bernoulli GLM Raw task, that E-C2ST-NF scales well to the high- dimensional observation space (faster convergence maximum TPR compared to Bernoulli GLM), while the normal E-C2ST and local-HPD significantly lose in statistical power. It should also be noted that local-HPD performs significantly worse in medium dimensions (cf. Bernoulli GLM or even SLCP) than in low dimensions (Caussian Mixture and Two Moons), though this could be because of the complex posterior structure.
F.2 Accuracy of (-C2ST(-NF) w.r.t. the true C2ST
Figures and 5 compare our method to the oracle C2ST, but only in terms of statistical power, as the local analysis is limited to the averaged results over 10 different reference observations xg- We here provide a more detailed local analysis that examines how the (-C2ST(-NF) test statistics correlate with those from the oracle C2ST, by plotting them against each other. The results obtained for the above mentioned 10 reference observations are shown in Figure 6a. To allow for more robust conclusions, we also show results obtained for total of 100 different reference observations (cf. Figure 6b), as well as quantitative results on the statistical significance of the correlation indices in Table
Overall, the scattered points are not too far from the diagonal, which indicates that there is some correlation between the test statistics for (-C2ST(-NF) and those from the oracle C2ST. This cor- relation becomes weaker when Ntrai becomes larger, since the test statistics in these cases tend to zero and start be confused with noise. This observation is consistent with the results in Table 3, showing the p-values of the Pearson test, a standard tests for the statistical significance of the correlation indices between the scores.
Another general trend is that the scattered points tend to be below the diagonal, indicating that (-C2ST(-NF) is less sensible to local inconsistencies than the oracle C2ST. This behaviour was expected, as C-C2ST is trained on the joint and thus less precise. Interestingly this doesn't 't apply to
Local methods TPR < 1 at Nuxin = 1000 (see Column 2), which means that the classification task harder and requires more data: local methods never reach maximum TPR = 1, and even the oracle C2ST (with MSE test statistic) needs Neal 2000, least four times more than for the other tasks.
25

---

the Gaussian Mixture task. This could be due to big variability in the local consistency of while trained on the joint (-C2ST(-NFF) could overfit on the "bad" observations, resulting in higher test statistics for observations where the true C2ST statistic would be small.
Finally, across all benchmark examples observe results that consistent with the ones from Figure 5: higher correlation means higher statistical power.
Two Moons
SLCP
Gaussian Mixture Gauss. Linear Unif.
Bemoulli GLM
Bermoulli GLM Raw
0.25
(-C2ST (ba)))
1000 T0000 1X0X0)
0.125
(-C2ST-NE (hos.)
1000
10000
0.125
TOUNID
0.125
0.25
0.125
0.25 0
0.125
0.25
0.125
0.250
0.125
0.25
0.125
0.25
oracle C2ST (bwr)
oracle CIST (bar)
oracle C2ST (uerr
ocacle C2ST (var)
oraclo C2ST (fum))
oracle C2ST fasr)
(a) 10 observations (same as in benchmark experiment
Figures and 5)
Two Moons
SLCP
Gaussian Mixture
Gauss. Linear Unif.
Beroulli GLM
Beroulli GLM Ruw
(-C2ST (bas)
0.125
F-C2ST-NE (ba))
0.25
0.125
0
0.125 oracle (sse)
0.25
0.125
0,250
0.125
0.25 0
0.125
0.25.
0.125
0.25
0.125
0.25
oracle (bxe)
ocacle (uaz))
C2ST (sae)
C2ST (om)
C2ST (iese)
(b) 100 ohservations (generated via sbihm package)
Figure 6: Accuracy correlation of E-C2ST(-NF) w.r.t. the oracle C2ST. We show seatter plots for all sbibm bm examples on (a) 10 and (b) 100 different reference observations. Each point corresponds to one observation and represents the MSE test statistic obtained for the oracle C2ST (x-axis) vs. our (-C2ST(-NF) method (y-axis). The diagonal represents the case where (-C2ST(-NF) C2ST. The closer points are to the diagonal, the more accurate E-C2ST w.r.t. C2ST.
Ntrain
10a
103
10
105
SLCP
10-4 /0.12
10-3/0.00
0.31 0.82
0.21 0.40
Two Moons
10-27 10
10-3
10-16 10-11
0.052/ / 10-5
Gaussian Mixture
10-8 / 10
10-0
0.006 0.35
10-14 / 0.006
Gauss. Linear Unif.
10-13/10-11
407/1109
0.42 0.002
0.68 0.87
Bernoulli GLM
10- 10-5
10-10/10-4
0.67 0.18
0.39 0.31
Bernoulli GLM Raw
0.03 / 0.37
10-8/10-1
10-4 /0.04
0.92 0.06
Table 3: P-values of the Pearson test of non-correlation between the E-C2ST( / -NF) and the oracle C2ST MSE test statistic. Obtained for 100 observations (as plotted Figure 6b) using scipy.star -ttssttttsppessonrr..e. pearsonr. Blue values indicate the cases for which the Pearson test rejects the null hypothesis of non-correlation with 95% confidence (i.e. there is significance evidence for correlation).
26

---

F.3 Benefit of the -NF version
The results of all benchmark experiments indicate that the -NF version of E-C2ST works better when the (true) posterior distribution of the model is "more complicated" than Gaussian distribution (see Table at the beginning of Appendix F). This for example the case for the Two Moons and SLCP tasks: the posterior distributions are globally multi-modal and locally structured. We observe in Column 3 of Figure 2 in the main paper that the -NF version requires less samples (i.e. lower Neal) to reach maximum power/TPR, This is also the case for the additional Bernoulli GLM task (see Column 3 of Figure
Appendix F.1). In contrast, for Gaussian Mixture and Gaussian Linear Uniform, where the posterior is Gaussian, the normal &-C2ST is powerful even better than its -NF counterpart.
Finally, we refer the reader to the analysis in Section F.1 to point out an interesting observation: for the Bernoulli GLM. the -NF version scales much better to high dimensional observation spaces than the normal E-CSST. Note that this does not allow us to make any general conclusions, but might be worth further investigating this result.
F,4 Advantages of (-C2ST w.r.t. local-HPD
First of all, it is important to mention that having uniform HPD-values is not sufficient condition for asserting the null hypothesis of consistency (see end of Section 3.3 in [48]). This is clear disadvantage compared to our proposal, which provides necessary and sufficient proxy for inspecting local posterior consistency.
Furthermore, the HPD methodology summarizes the whole information concerning 0 into a single scalar, while in E-C2ST handle the f-vector in its multivariate form. In medium-high O-dimensions (m > 2 as in Bernoulli CLM) or for complex posterior distributions (SLCP), such summarized information might discard too much information and not be enough to satisfactorily assess the consistency of the posterior estimator. Indeed, the tasks where loca/-HPD has similar statistical power E-C2ST either in low dimensions (Two Moons) and or have Gaussian posterior (Gaussias Mixture/ Gaussian Linear Uniform). See Table 2 for a summary of those results w.r.t. data dimensionality and posterior structure. For a detailed analysis of the scalability to high dimensional data spaces see Appendix F.I.
Finally, loC in its naive implementation is much less efficient than L-C2ST (see Appendix C). Note however that a new "amortized" version of local-HPD has recently been proposed [9] and could be interesting to look at.
27


Reviews:
[{"review_id": "xeks3BNEm1", "reviewer": "Reviewer_Esye", "summary": "An algorithm for evaluating amortized posterior estimators $q(\\theta|x)$ in simulation-based inference (SBI) is proposed. In the setting considered, we have the ability to sample $p(\\theta)$ and to sample $p(x|\\theta)$ but not to evaluate its density. The algorithm modifies the common classifier two-sample test by learning to discriminate between samples $(\\theta,x)$ sampled ancestrally (though simulation) from $p(\\theta,x)$ and those where $\\theta$ is then resampled from $q(\\theta|x)$. An extension to normalizing flows is also considered, where the discriminator works in the latent space rather than in the theta space. There are strong results on several SBI benchmarks and real-world problems in terms of statistical power and runtime.", "questions": "Please see above.", "limitations": "yes", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Please note that this is a relative outsider's point of view, as I do not directly work on SBI.\n\nI find the paper well-written and was able to understand the main points and the math. In particular, the first few pages are a good introduction to the problem area that will be accessible to any reader who is familiar with Bayesian inference and hypothesis testing.\n\nThe results and their significance are explained well and give intuitions about why the proposed algorithm works well on the chosen problems.\n\nThe provided code is a helpful addition.", "weaknesses": "I do not see any major weaknesses, but have a few questions.\n- Could you please comment on the scalability of the algorithm? What difficulties do you foresee when scaling to high-dimensional theta or problems when the prior is wide? (Same question about high-dimensional observations. Is more risk of overfitting to the class split of $x$?)\n- It would be interesting to understand better the benefit of the -NF version of the algorithm beyond the reasons described in L180, which seem to explain only computation cost improvements. Is it expected that the classification boundary is smoother in the latent space than in the theta space? Maybe this could be illustrated, e.g., on the two moons example.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for finding our paper \"clear and well-written\", as well for having checked our code. We now answer the two main questions raised by the reviewer. Please note that we have added results on more SBI benchmark examples that you will find in Figures 1 and 2 of the PDF attached to the rebuttal. These examples are described and the results summarized in the global response to all reviewers.\n\n\n**Question 1: Scalability of the algorithm to high-dimensional theta or x**\n\nAn obvious answer would be to say that in higher dimensions more samples are needed: larger $N_{\\mathrm{cal}}$ for $\\ell$-C2ST(-NF) to converge to the oracle C2ST. In our setting, this would mean better convergence to the *optimal Bayes classifier*, which is known to be the same for $\\ell$-C2ST and oracle C2ST. Note, however, that the MSE test statistic for $\\ell$-C2ST is defined by the predicted class probabilities and not the accuracy of the classifier. Therefore, adding a regularizer or maybe a calibration step (e.g. using methods from the conformal inference literature) might be necessary when working with data defined in high dimensions, so as to prevent overfitting and overconfident predictions.\n\nIn the specific case of SBI, the dimension of the parameter space ($m$) is typically of order $10^0$ to $10^1$ and $m \\approx 10^2$ is already often considered as high dimensional. The observation space, however, can be high-dimensional (e.g. time-series), but summary statistics are often used to reduce the dimension of the observations to the order of $d \\approx 10^1$.\n\nIn our initial submission, we show results for rather low-dimensional tasks ($m=2, d=2$) for Two-Moons and $m=5, d=8$ for `SLCP`). They illustrate how our method behaves in terms of the difficulty of the inference task (e.g. complex posteriors as in `SLCP`). To demonstrate how the results change with respect to the dimensionality of the data, we've added more tasks of low and medium dimensionality: `Gaussian Mixture` ($m=2, d=2$) `Gaussian Linear Uniform` ($m=10, d=10$) and `Bernoulli GLM` ($m=10, d=10$). To analyze how our method scales to **high dimensional observation spaces** only (without parameter-space / task variability), we’ve also added the `Bernoulli GLM Raw` task.  It considers *raw* observation data ($d=100$), as opposed to *sufficient summary statistics* in `Bernoulli GLM` ($d=10$). We refer the reader to the global response for all reviewers for a description of all benchmark examples (in terms of dimensionality, posterior structure, challenges) and a summary of all experimental results. \n \nIn the attached PDF , Figure 1 shows results obtained for the new benchmarks, therefore extending Figure 2 of the initial submission. We can see in Column 3 of both Figures that  for $\\ell$-C2ST to converge to the oracle C2ST (at maximum power $\\mathrm{TPR} = 1$), less samples are needed in the low-dimensional `Two Moons` and `SLCP` tasks than for the medium-dimensional `Bernoulli GLM` task: $N_{\\mathrm{cal}} \\approx 2000$ vs. $N_{\\mathrm{cal}} \\approx 5000$. This confirms our intuition stated above. \n\nNote that we did not include the `Gaussian Mixture` and `Gaussian Linear Uniform` tasks in this analysis, as they are not comparable to the other tasks: the TPR of $\\ell$-C2ST at $N_{\\mathrm{train}}=1000$ is smaller than 1 (see Column 2). The classification task is thus harder and more samples are required to converge to reach maximum TPR: the oracle C2ST now requires $N_{\\mathrm{cal}}=2000$ and *local* methods never reach maximum TPR (see Column 3). Here, the difficulty of the classification task has more impact on statistical power than the dimensionality of the task (e.g. the convergence to maximum TPR is slower in `Gaussian Mixture` than in `Bernoulli GLM` of higher dimension).\n\nInterestingly, we observe in the  `Bernoulli GLM Raw` task, that $\\ell$-C2ST-NF scales well to the high-dimensional observation space (faster convergence to maximum TPR compared to the `Bernoulli GLM` task), while the normal $\\ell$-C2ST and *local*-HPD significantly lose in statistical power. \n\n\n**Question 2: Benefit of the -NF version** \n\nThe numerical illustrations in our manuscript indicate that the -NF version of our statistical test works better when the (true) posterior distribution of the model is \"more complicated\" than a gaussian distribution. This is the case for the `Two Moons` and `SLCP` tasks: the posterior distributions are globally multi-modal and locally structured (cf. task descriptions in the global response to all reviewers). We observe in Column 3 of Figure 2 in the main paper that the -NF version requires **less samples** (i.e. lower $N_{\\mathrm{cal}}$) to reach maximum power/TPR. This is also the case for the additional `Bernoulli GLM` task (see Column 3 of Figure 1 in the attached PDF). In contrast, the additional `Gaussian Mixture`and `Gaussian Linear Uniform` tasks, where the posterior is a gaussian distribution, the normal $\\ell$-C2ST is as powerful or even better than its -NF counterpart (see Column 3 of Figure 1 in the attached PDF).  \n\nWe have also observed that $\\ell$-C2ST-NF yields test statistics that correlate more closely with the oracle C2ST in situations where the true posterior distribution can be sampled with MCMC. The experiments and illustrations related to this finding can be found in Figure 2 of the attached PDF and are described in our response to Reviewer `ukES02`.  \n\nFinally, we refer the reader to the previous question to point out an interesting observation: for the `Bernoulli GLM`, the -NF version scales much better to high dimensional observation spaces than the normal $\\ell$-C2ST. Note that this does not allow us to make any general conclusions, but it might be worth further investigating this result.\n"}, {"review_id": "mM8maV4Hof", "reviewer": "Reviewer_ukES", "summary": "The authors present a new diagnostic tool for simulation-based inference. The method learns the C2ST without knowledge of the true posterior distribution. If the density estimator is a normalizing flow, the authors propose to perform the classification in latent space.", "questions": "- Do I understand correctly that, for large N_H, one has to train the classifier multiple times? If yes, couldn’t this make the method less efficient than local-HPD?", "limitations": "The authors state limitations of their method.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "**Originality**: The method tackles and important issue for simulation-based inference. The method is new and potentially very useful. The extension for normalizing flows is interesting and creative. If L-C2ST works well, I could see this be adopted by the community.\n\n**Quality**: The method is rigorously derived.\n\n**Clarity**: Everything is explained well and with sufficient detail. Figures are clear.", "weaknesses": "**Quality**: \nMy main concern with this paper is that it does not sufficiently demonstrate how well the method actually works. In particular, I believe that the following crucial questions remain unanswered:\n\n- Does the method indeed capture the true C2ST for different `x`?\n\nFor example, the (decent) TPRs shown in figure 2 could also be obtained if L-C2ST converges onto the average C2ST across x_o. The central claim of the authors is that L2CT is local, but this claim needs substantially more empirical evidence.\n\nTo address this concerns, I would (for example) suggest to add a scatter plot of L-C2ST vs true C2ST on a task where the ground truth is available by MCMC or analytically\n\n- How well does the method scale with parameter and data dimensionality\n\nSince C2ST is trained on theta and x as input, I could imagine it to not scale well to high-D x or theta (or require **even more** simulations). All tasks in the paper are low-D though. I suggest that the authors add an additional benchmark task with higher data and parameter dimensionality.\n\n- Is the method actually better than local-HPD\n\nOn one task, L-C2ST performs better than local-HPD and worse on another one. Which method should be preferred? The paper also claims that local-HPD is much less efficient than C2ST but this is never empirically shown. In particular, while a naive implementation of local-HPD might be slow, amortizing over the confidence alpha should be trivial and make local-HPD as fast as C2ST. Please correct me if I am wrong or, otherwise, clarify this in the paper.\n\nTo address this, I recommend adding an analysis of the runtimes of the algorithms and ideally have more than only two benchmark tasks.\n\nI believe that the above concerns are critical to the usefulness of the method, but I would be willing to significantly increase my score if these things are addressed.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for recognizing the potential and usefulness of the method that we propose, as well as appreciating its extension to the case of normalizing flows. We also thank the reviewer for the very deep and interesting questions related to our method. We have followed the reviewer's suggestions and have added results on more SBI benchmark examples, as well as new experiments for analyzing our results (cf. Figures 1 and 2 in the attached PDF). These examples are described and the results summarized in the global response to all reviewers.\n\n**Question 1: Does the method indeed capture the true C2ST for different x?**\n\nFigure 2 in the submitted paper compares our method to the oracle C2ST, but only in terms of statistical power, as we limited the *local* analysis to the averaged results over 10 different reference observations. Following the reviewer's suggestion, we've added scatter plots to examine how the values of the $\\ell$-C2ST(-NF) test statistic correlate with those from the oracle C2ST. The results are presented in Figure 2.a of the attached PDF for the same 10 reference observations $x_o$ initially used. To allow for more robust conclusions, we've also included scatter plots using 100 reference observations (cf. Figure 2.b). Please note that we have also considered additional benchmark examples.\n\nOverall, the scattered points are not too far from the diagonal, which indicates that the test statistics for $\\ell$-C2ST(-NF) correlate quite well with those from the oracle C2ST. Also, as $N_{train}$ increases, the scattered points become closer to zero and are more concentrated. This result comes with no surprise: when the posterior approximate is consistent, the test statistics should be close to zero. If the approximate is poorly estimated, the statistics should deviate from zero and the points can take on different values. We observe that the -NF version seems to be slightly closer to the oracle C2ST than the normal $\\ell$-C2ST.\n\nIn examples `Two Moons`, `Gaussian Linear Uniform`, and `Bernoulli GLM`, the scatterplots display points that are close to the diagonal. In tasks `SLCP` and `Bernoulli GLM Raw`, most scattered points are below the diagonal, especially for higher $N_{train}$ values. This means that $\\ell$-C2ST has lower test static values, meaning that the null hypothesis will likely not be rejected: it is, therefore, less sensitive to local differences between $p$ and $q$ than the oracle C2ST, which is consistent with Figure 1 of the attached PDF. Intuitively, this can be explained by the fact that $\\ell$-C2ST is trained on the joint pdf and is thus less precise. \n\nThe results for the `Gaussian Mixture` task however deviate from the general trend, as in Figure 1 of the attached PDF. Maybe due to big variability in the *local* consistency of $q$. Unlike the true C2ST, $\\ell$-C2ST is trained on the joint data space and could therefore overfit on the \"bad\" observations, resulting in higher test statistics for observations where the true C2ST statistic would be small. \n\n**Question 2: How well does the method scale with parameter and data dimensionality**\n\nPlease see our response to Question 1 from Reviewer `Esye`.\n\nIn our initial submission, the performance of C2ST on the SBI tasks that we considered was mainly impacted by the structure of the parameter space and the corresponding shape of the posterior distribution, since the examples were of rather low dimensionality. We have added more SBI tasks with varying sizes of observation/parameter space to give a clearer picture of the performance of our method.\n\nIt should be noted that local-HPD performs significantly worse in medium dimensions (cf. `Bernoulli GLM` or even `SLCP`) than in low dimensions (`Gaussian Mixture` and `Two Moons`), though this could be because of the complex posterior structure. Also, $\\ell$-C2ST-NF scales well to the high-dimensional observation space of `Bernoulli GLM`, while *local*-HPD significantly loses statistical power. \n\n\n**Question 3: Is the method actually better than local-HPD?**\n\nFirst of all, it is important to mention that having uniform HPD-values is not a sufficient condition for asserting the null hypothesis of consistency. This is mentioned by Zhao et al. (2021) (see end of section 3.3) and is a clear disadvantage compared to our proposal, which provides a necessary and sufficient proxy for inspecting local posterior consistency.\n\nFurthermore, the HPD methodology summarizes the whole information concerning $\\theta$ into a single scalar, while in $\\ell$-C2ST we handle the $\\theta$-vector in its multivariate form. In high $\\theta$-dimensions (`Bernoulli GLM`) or for complex posterior distributions (`SLCP`), such summarized information might discard too much information and not be enough to satisfactorily assess the consistency of the posterior estimator. Indeed, the only task where local-HPD outperforms $\\ell$-C2ST is `Gaussian Mixture` (low dimension and gaussian posterior).\n\nFinally, as mentioned by the reviewer, local-HPD in its naive implementation is much less efficient than $\\ell$-C2ST (see Appendix A.5 in submitted paper). A new version of *local*-HPD with amortized $\\alpha$ has recently been proposed by DEY et al. (2023). However, we were not able to verify how well it works in time for the submission deadline.\n\n**Question 4: For large N_H, one has to train the classifier multiple times? (...) less efficient than local-HPD?**\n\nN_H is the number of times we compute the test statistic under the null hypothesis in order to compute p-values. The number of classifiers we need to train depends on how many we need to compute the test statistic ($N_H$ for $\\ell$-C2ST vs. $N_H \\times n_{\\alpha}$   for local-HPD). In summary, if $\\ell$-C2ST is more efficient in computing a single test statistic, it will also be more efficient to compute $N_H$  test statistics. "}, {"review_id": "4Sp2e5rjXh", "reviewer": "Reviewer_Ppab", "summary": "This paper defines a new diagnostic to check posteriors learned from simulation-based inference.  The main issue in this setting is that one does not have access to likelihoods, and so one needs diagnostics that do not rely on being able to calculate likelihoods.  Here the authors present L-C2ST, which essentially tries to train a classifier to distinguish the _joint_ space over both parameters _and_ observations into whether a given point was more likely under the true joint distribution, $p(x, \\theta)$ or under a variational approximation $q(\\theta | x) p(x)$.  It turns out that a Bayes-optimal (probabilistic) classifier will assign equal weight to the two distributions for all $\\theta$ for a given $x_0$ if and only if $q(\\theta | x_0) = p(\\theta | x_0)$ (up to some technical considerations).  As a result, the authors propose training such a classifier and then seeing how much its weights deviate from $0.5$ and using that as a measure of how well $q$ matches $p$.  Finally, the authors develop a hypothesis testing framework (to account for randomness in the trained classifier) to determine how significant a given deviation is.", "questions": "* Using the same $X_n$'s in Algorithm 1 for draws from $p(\\cdot | X_n)$ and $q(\\cdot | X_N)$ introduces dependencies in the training data for the classifier.  Can that somehow be leveraged when training the classifier (similar to how in linear regression, we would want to perform GLS instead of OLS to increase the power when data points have shared structure like this).  Along similar lines, if it is much cheaper to sample from $q(\\theta | x)$ than it is to sample from $p(x, \\theta)$, then could it help to somehow incorporate more samples from $q(\\theta  | X_n)$ in the algorithm (or perhaps the whole distribution somehow)?  Obviously this would require accounting for the dependencies (perhaps via some kind of downweighting) in the training of the classifier.  This is obviously not necessary to implement for this paper, but I am curious to hear the authors thoughts.\n* Algorithm 1 also seems computationally intensive due to the permutation step.  Would it be possible to do something like split-sample conformal inference to speed up the p-value calculations?  Again, obviously not necessary for the paper, I'm just curious.\n", "limitations": "I don't see any potential for direct societal harm from this work.  The authors may wish to include a limitations section that highlights failure modes of their method.", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 4, "contribution": 3, "strengths": "* Trustworthy checks of how good a given variational posterior is for real settings where we do not have access to the true posterior are sorely needed, and this work is an important step in that direction.\n* The paper is well-written, well-motivated, and clearly explained.\n* The theoretical results complement and motivate the proposed algorithm well (but see below for some minor technical concerns).\n* The proposed algorithm is extremely simple to describe, making it practical and elegant.", "weaknesses": "* My major concern is that I am not totally convinced that learning the proposed classifier is any easier than learning the variational posterior.  As such, I would be concerned that diagnostics showing a good posterior fit could be driven by a far-from-optimal classifier.  My intuition for why this should be a hard problem is as follows: for _any_ amortized posterior $q(\\theta | x)$,  if we can learn a Bayes-optimal classifier $d^*$, then one could recover the true posterior for any $x$ via $p(\\theta | x) = q(\\theta | x) / (1-d^*(\\theta, x))$.  As a result, learning $d^*$ for _any_ $q$ is equivalent to learning the true posterior.  This is somewhat borne out in the third column of figure 2, where one needs roughly the same number of samples to learn a good posterior as one needs in order to learn a good classifier.  I'm wondering if somehow the proposed approach is essentially equivalent to training two different amortized variational posteriors on independent data and then seeing how well they match.  Some further simulations showing that one problem is easier than the other would alleviate some of my concerns.  I should also say that I think that even if my concerns are valid, the diagnostic is still interesting and useful.\n* I felt that the benchmarking could have been improved slightly.  In particular, some empirical relationship between the proposed measure of fit and existing measures would be nice.  For example, the authors could take the true posterior in a toy model where it is known, generate $q$ by perturbing that posterior in a systematic way and compare how increasingly large perturbations affect both the L-C2ST measure as well as things like KL, reverse KL, TV, and so on.  The results in Figures 1 and 2 hint at this, but there are a lot of moving parts (especially in Figure 2) making it difficult to determine whether to contribute power (or lack thereof) to failure of the SBI (i.e., $q$ is bad) or the classifier (i.e. $r$ is good).\n* A minor technical point throughout the proofs of the theorems is that they rely on the null hypothesis $\\mathcal{H}_0$ holding on a set of strictly positive measure.  That is, if $\\mathcal{H}_0$ fails on a set of measure zero then the proofs would break down.  In particular, in that case $\\mathcal{A}$ would have measure zero and the final integral after line 454 would be exactly zero.  I tend to feel that failures on sets of measure zero have no practical implications, and so this is a minor point. But, if one is going to state theorems and proofs, it would be good to be rigorous.\n* Similarly, Lemma 1 should not be phrased in terms of \"points\" as those would have zero measure, and instead should say \"If there exists a set $\\mathcal{A} \\subseteq \\mathcal{S}$ with positive measure such that $p(\\theta) > q(\\theta)$ for all $\\theta \\in \\mathcal{A}$, then there exists a set $\\mathcal{B}$ with positive measure such that $p(\\theta') < q(\\theta')$ for all $\\theta'\\in \\mathcal{B}$.\n* Finally, there are similar  (and equally minor) issues appearing in the proofs of Theorems 2 and 3.\n\nTypos:\n* Line 103: \"may lack of\" --> \"may lack\"\n* Line 240: \"in average\" --> \"on average\"\n* Line 251: Missing period before \"The local PP-plot\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for saying that our paper is well written and clearly explained. We also appreciate that he found our method \"simple and elegant\" and recognized that we aim at an important practical question for simulation-based inference. Here below we address some of the remarks and questions raised by the reviewer, which were very interesting and pointed towards exciting perspectives that we shall investigate in the near future.\n\n**Remark 1: My major concern is that I am not totally convinced that learning the proposed classifier is any easier than learning the variational posterior.**  \n\nThis is an interesting and true remark. I am not sure whether one tool could be considered as strictly easier than another. The results in Figure 2 of the initial submission indeed show that the sample-sizes are roughly the same for the `Two Moons` task. For the `SlCP` task however, the oracle C2ST statistic never reaches close-to-zero values, implying that much more samples (N_train) are needed for the posterior estimator to converge, than for $\\ell$-C2ST to reach maximum TPR (N_cal). Results for additional benchmarks can be found in Figure 1 of the attached PDF: for the `Gaussian Linear Uniform` example classification appears harder than posterior estimation, for the `Bernoulli GLM` task however the posterior estimator appears to converge at $N_{\\mathrm{train}}=10^5$ (see oracle C2ST curve in Column 1), while the maximum TPR is reached for $N_{\\mathrm{cal}}=10^4$ (see $\\ell$-C2ST-NF curve in Column 3). \n\nBeing concerned about false conclusions due to a far-from-optimal classifier is totally justified. One should always make sure the classifier is “good enough” before using it as a diagnostic tool. If it’s not, using another diagnostic tool would be better. One way to do this would be to use cross-validation as a tool to select the classifier with highest accuracy for the observed data (the accuracy should never be exactly at 0.5) . Note, however, that the MSE test statistic for $\\ell$-C2ST is defined by the predicted class probabilities and not the accuracy of the classifier. Therefore one should also check how well the classifier is calibrated. \n\nWe have thought about this aspect of training a second posterior estimate instead of using a classifier to assess the consistency of the first posterior estimator. However, we decided that this would be redundant and preferred using two different tools to play against each other.  Furthermore, the choice of a binary classifier was not really motivated by the “easy to learn” aspect, but mainly the popularity of binary classification, the fact that it’s easy to understand and has a much richer and stable literature than deep generative models. We think that this way our method will be more convincing for people from many different fields, instead of just the sbi / bayesian inference community. \n\n\n**Remark 2: I felt that the benchmarking could have been improved slightly. In particular, some empirical relationship between the proposed measure of fit and existing measures would be nice (...).**\n\nWe duly note this remark and thank the reviewer for pointing it out. Indeed, it would be nice to have some easier to understand and interpretable experiments, where less factors have an influence on the results. Unfortunately we did not have the time to design a new experiment.\n \n\n**Remark 3: A minor technical point throughout the proofs of the theorems is that they rely on the null hypothesis holding on a set of strictly positive measure (...)**\n\nWe will correct the demonstrations in our camera ready version. We thank the reviewer for pointing this out!\n\n**Question 1: Using the same $X_n$'s in Algorithm 1 for draws from $p(.|X_n)$ and $q(.|X_n)$ introduces dependencies in the training data for the classifier. Can that somehow be leveraged when training the classifier (...). Along similar lines, if it is much cheaper to sample from $q(\\theta|x)$ than it is to sample from $p(x, \\theta)$, then could it help to somehow incorporate more samples from $q(\\theta|X_n)$ in the algorithm (or perhaps the whole distribution somehow)? (....)**\n\nUsing the same $X_n$’s for draws from $p$ and $q$ is a sensible and interesting topic. In Appendix A.7 in the initial submission we address this and show that in this case the theoretical cross-entropy loss is the same as if the $X_n$’s were drawn independently. \nHowever, we have thought about adding more samples from $q$ to improve our method as a follow-up to this work. When considering Normalizing Flows as posterior estimators, we would even have access to the likelihood. That could definitely help. My concern is that the classification problem becomes unbalanced and we might overfit on data from $q$ w.r.t to data from $p$. Interestingly the recent work Discriminative Calibration  from Yoa and Domke (2023) proposes an extension of C2ST for SBI that tries to leverage additional information of $q$, but without being *local* (main contribution of $\\ell$-C2ST). \n\n**Question 2: Algorithm 1 also seems computationally intensive due to the permutation step. Would it be possible to do something like split-sample conformal inference to speed up the p-value calculations?**\n\nThis is a very interesting question. We have taken a look at conformal inference to include a calibration step in the predict-prob from the classifiers. Indeed, using a  technique with finite-sample guarantees based on ordered statistics could probably allow us to avoid the rather time consuming permutations in the null hypothesis. However, it should be noted that conformal prediction is known to have difficulties in providing intervals with fixed conditional coverage, so for our local procedure this could be a problem. We will take a further look on this in the near future.\n\n**Limitations:** The authors may wish to include a limitations section that highlights failure modes of their method.\n\nWe will try to include one in the final camera ready version of this article.\n"}]

Author rebuttal:
We thank the reviewers for their thorough reading of our manuscript and their very interesting questions. We have had very positive remarks concerning the clarity of our text and the potential of our method to the SBI community. We are very grateful for all this. Some interesting questions that were raised concerned how our method behaves when the dimensionality of the data increases, whether the test statistics of $\ell$-C2ST really capture the true test statistics (i.e. those from the oracle C2ST) for different values of conditioning observation $x$, and whether our method is indeed superior to local HPD. We address all these questions (and some more) below.

We have followed suggestions and have added results on further examples from the SBI benchmark. In the attached PDF, Figure 1 extends Figure 2 of the submitted version, with the intent of better validating our method on varying dimensions for observations and parameters. We have also included results for a new experiment in Figure 2 of the PDF, showing the correlation between the test statistics of $\ell$-C2ST(-NF) and those from the oracle C2ST for different conditioning observations. 


**Description of the used SBI benchmark examples from `sbibm`  and summary of main results:**

`Two Moons`:

- **Dimensions ($\theta, x$):** (2,2)
- **Posterior structure:** bi-modal, crescent shape 
- **Challenge:** globally and locally structured

- **$\ell$-C2ST(-NF) vs. local-HPD**  (power, Figure 1): faster convergence to max power, but worse for high $N_{train}$
- **$\ell$-C2ST-NF vs. L-C2ST** (power, Figure 1): better (less samples for higher power)
- **$\ell$-C2ST(-NF) vs. C2ST** (relative position to diagonal, Figure 2): good

`SLCP`:

- **Dimensions ($\theta, x$):** (5,8)
- **Posterior structure:** 4 symmetrical modes
- **Challenge:** posterior designed to be complex

- **$\ell$-C2ST(-NF) vs. local-HPD**  (power, Figure 1): better
- **$\ell$-C2ST-NF vs. L-C2ST** (power, Figure 1): better
- **$\ell$-C2ST(-NF) vs. C2ST** (relative position to diagonal, Figure 2): lower

`Gaussian Mixture`:

- **Dimensions ($\theta, x$):** (2,2)
- **Posterior structure:** 2D gaussian
- **Challenge:** one of the gaussians in the mixture has much broader covariance than the other.

- **$\ell$-C2ST(-NF) vs. local-HPD**  (power, Figure 1): slightly worse
- **$\ell$-C2ST-NF vs. L-C2ST** (power, Figure 1): same
- **$\ell$-C2ST(-NF) vs. C2ST** (relative position to diagonal, Figure 2): lower, then higher

`Gaussian Linear Uniform`:

- **Dimensions ($\theta, x$):** (10,10)
- **Posterior structure:** multivariate Gaussian
- **Challenge:** scaling to higher dimensions

- **$\ell$-C2ST(-NF) vs. local-HPD**  (power, Figure 1):  similar
- **$\ell$-C2ST-NF vs. L-C2ST** (power, Figure 1): worse
- **$\ell$-C2ST(-NF) vs. C2ST** (relative position to diagonal, Figure 2): good

`Bernoulli GLM`:

- **Dimensions ($\theta, x$):** (10,10)
- **Posterior structure:** unimodal, concave
- **Challenge:** scaling to higher dimensions

- **$\ell$-C2ST(-NF) vs. local-HPD**  (power, Figure 1): better
- **$\ell$-C2ST-NF vs. L-C2ST** (power, Figure 1): better
- **$\ell$-C2ST(-NF) vs. C2ST** (relative position to diagonal, Figure 2): good / slightly lower

`Bernoulli GLM Raw`:

- **Dimensions ($\theta, x$):** (10,100)
- **Posterior structure:** unimodal, concave
- **Challenge:** Scaling to high dimensional observation spaces

- **$\ell$-C2ST(-NF) vs. local-HPD**  (power, Figure 1):  better
- **$\ell$-C2ST-NF vs. L-C2ST** (power, Figure 1): better
- **$\ell$-C2ST(-NF) vs. C2ST** (relative position to diagonal, Figure 2): slightly lower

Meta-review: The paper introduces a diagnostic tool, L-C2ST, for evaluating the quality of posterior at any given observation in simulation-based inference (SBI). This is achieved by training a classifier to discern between samples from the true joint distribution and those from a variational approximation. The work delves into the Bayes-optimal classifier's significance in the context, and an hypothesis testing framework is developed to understand the deviation's significance.

**Strengths:**
- Novel diagnostic tool for simulation-based inference, tackling a significant gap in current SBI methods.
- Strong theoretical grounding that is explained well and complements the proposed method.
- The paper is well-written, making it accessible to readers not deeply familiar with SBI.
- Provides open-source code.

**Weaknesses:**
- Several reviewers raised concerns about the method's scalability with high-dimensional data or parameters, an aspect not sufficiently addressed in the paper.
- Concerns about the actual utility and practicality of training the classifier: Is it inherently easier than obtaining a good variational posterior?
- The paper doesn't provide enough evidence about the true performance of the method, especially in capturing the true C2ST for different x values and the comparison to existing methods like local-HPD.
- Benchmarking could be more extensive, possibly comparing the proposed diagnostic with other measures using a known true posterior.

The paper has significant potential and originality, though there remain concerns particularly in its empirical validation and scalability. Given its potential and the positive reviews, I recommend accept with a strong suggestion for revisions to address the highlighted weaknesses.
Decision: Accept (poster)