Discriminative Calibration: Check Bayesian Computation from Simulations and Flexible Classifier
Authors: Yuling Yao
Justin Domke
Venue: NeurIPS 2023 poster (2023)
Keywords: simulation based calibration
simulation based inference
Bayesian computation
diagnostics
classifier two-sample test
likelihood-free
TL;DR: We use a classifier+statistical feature approach to diagnose Bayesian computations, applicable to MCMC, variational and simulation based inference. We provide a high-power test and a divergence estimate.

Abstract: To check the accuracy of Bayesian computations, it is common to use rank-based simulation-based calibration (SBC). However, SBC has drawbacks: The test statistic is somewhat ad-hoc, interactions are difficult to examine, multiple testing is a challenge, and the resulting p-value is not a divergence metric. We propose to replace the marginal rank test with a flexible classification approach that learns test statistics from data. This measure typically has a higher statistical power than the SBC test and returns an interpretable divergence measure of miscalibration, computed from classification accuracy. This approach can be used with different data generating processes to address simulation-based inference or traditional inference methods like Markov chain Monte Carlo or variational inference. We illustrate an automated implementation using neural networks and statistically-inspired features, and validate the method with numerical and real data experiments.

Diseriminative Calibration: Check Bayesian Computation from Simulations and Flexible Classifier
Yuling Yao Flatiron Institute, New York, NY 10010. yyao@yyao. dev
Justin Domke University of Massachusetts, Amherst, MA 01002. domke@cs umass edu
Abstract
To check the accuracy of Bayesian computations, it is common simulation-bart -basad calibration (SBC). However, SBC has drawbacks: The test statis- tic is somewhat ad-hoc, interactions a challenge, and the resulting p-value is not a divergence metric. We propose to replace the marginal rank test with flexible classification approach that learns test statistics from data. This measure typically has a higher statistical power than the SBC test and returns an interpretable divergence measure of miscalibration, com- puted from classification accuracy. This approach can be used with different data generating processes address simulation-bai baasd inference or traditional inference methods like Markov chain Monte Carlo or variational inference. We illustrate an automated implementation using neural networks and statistically-insee inspiied features, and validate the method with numerical and real data experiments.
tank-based
difficult to examine, multiple testing is
Introduction
Simulation based calibration (SBC) is default approach to diagnose Bayesian computation. SBC was originally designed to validate if computer software accurately draws samples from the exact posterior inference, such Markov chain Monte Carlo (MCMC, [4, 38, 26]) and variational inference [41]. With recent advances in amortized and simulation-basee inferences [6] and growing doubt on the sampling quality [23, 18], there has also been an increasing trend apply SBC likelihood-fie free inference such as approximate Bayesian computation [42] and normalizing-flor flow-baawwdd- [31, 1] neural
posterior sampling [22, 20), with wide range of applications in science [15, 7, 34].
Bayesian computation tries to sample from the posterior distribution p(fly) given data y. We work with the general setting where it may may not be possible to evaluate the likelihood p(yl0). Suppose we have an inference algorithm or software (0ll)) that attempts to approximate p(fly), and we would like to assess if this q is calibrated, meaning if (0lg)) Simulation based calibration involves three steps: First, we draw a from the prior distribution p(0). Second, we simulate synthetic observation y from the data model p(y]0). Third, given y we draw size-M posterior sample 0 ,, from the inference engine ((0l)) that we need to diagnose. SBC traditionally computes the rank statistic of the prior draw O among the samples, i.c. .M -- =1 1(0 < ). If the inference is calibrated, then given both 0 and 0m are from the same distribution p p00ly). hence with repeated simulations of (0, 3), we should expect such rank statistics to appear uniform, which can be checked by histogram visualization a formal uniformity test.
p(0 (y) for all possible y-
>M
Despite its popularity, rank-based SBC has limitations: We only compute the rank of univariate parameters. In practice, 0 and are high dimensional. We typically run SBC on each component of separately, this creates many marginal histograms and does not diagnose the joint distribution or interactions. We may compare ranks of some one-dimensional test statistics, but there is no method to
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Unssiication cramples
simulation table
0,
(1, Φ1
permutation test (p value)
label mapping d
0,2
l2, th2
classifier
02
y/1
prediction loss
0sM
Lx, tx
divergence D and its C.I.
feature engineering log log q. rank, etc.
Our discriminate calibration frumework has three modules (a) generate simulation table (0, y. 8). and map it into classification examples with some label and feature 6, (b) train a classifier to predict labels, (c)
Figure 1:
from the learned classifier, perform hypothesis testing and estimate divergence.
find the best summary test statistic. (ii) As long as we test multiple components of 0 or test statisties, directly computing the uniformity p-values is invalid (subject to false discovery) unless we make a multiple-testing adjustment, which drops the test power (subject to false negative) in high dimensions. (ii) Often we know inference is approximate, so the final goal of diagnostic is not to reject the null hypothesis of perfect calibration but to measure the degree of miscalibration. The p-value is not such measure: neither can we conclude an inference with p .02 is better than p = .01, nor connect the p-value with the posterior inference error. The evidence lower bound, a quantity common in variational inference, also does not directly measure the divergence due to the unknown entropy.
Heuristic: calibration via classification. To address all these drawbacks, while maintaining versatility in various computing tasks, we propose discriminative calibration, pragmatic and unified framework for Bayesian computing calibration and divergence estimation. The intuitive heuristic behind discriminative calibration is to use classifier to perform similarity tests-when two several distributions are similar, we cannot distinguish samples from them so the classification error is large. In Bayesian computation, we compare conditional distributions, the true p(0ly) and inferred ((ll)), to which we have access via simulations, and in some occasions explicit densities. It is natural to try to classify the samples drawn from p v.s. from q. and the ability to distinguish them would suggest miscalibration between p and q.
To formulate such heuristics into an algorithm, we design family of "label mapping" that prepares the simulations into classification examples that contain the label and feature. In Sec. 2, we first give four concrete label mapping examples, where the rank-based SBC becomes essentially a special case. Sec. 3 states the general theory: if we train a classifier to predict the labels, then the prediction ability yields a computable divergence from p(Oly) to (0ly)) In Sec. 4, we illustrate the practical implementation to get the divergence estimate, its confidence interval, and a valid hypothesis testing p-value. We explain why the learned classifier helps statistical power. In Sec. 4.1, we discuss the classifier design and incorporate extra known information such as ranks, likelihood, and log densities, whenever available, features. We also show method is applicable to MCMC without waste from thinning. We illustrate numerical and cosmology data examples in Sec. 5. We review other related posterior validation approaches and discuss the limitation and future direction in Sec. 6.
2 Generate labels, run classifier, and obtain a divergence
As with traditional SBC, we generate a simulation table by repeatedly sampling parameters 0 and synthetic data y from the target model p. i.c. draws (0, y) ~ p(0, v). Then, for each (0,v), y). we run the inference routine q to obtain a set of M IID approximate posterior samples 0,.-.. AN ~ q(fly). We wish to assess how close, on average (over different y), the inference procedure g(0ly) is to the true posterior ((0l)). Here, we observe that classification example-sets can be created in several ways and these produce different divergences. We generalize the four examples and the claims on divergence in Sec. 3.
Input:
one simulation run: (v, 0, OI...... Oax);
Output:
M + examples:
label
features ø
0
(0, v)
(01,0)
Example 1: Binary classification with full features. An intuitive way to
(0Me,) v)
estimate this closeness between ((lly) and p(Oly) is to train a binary classifier. Imagine creating a binary classification dataset of (t. Φ) pairs, where t a binary label, and & are features. For each (8, y) simulated from M + 1 pairs of examples are created. In the first, t = 0
and o- (0, v). In the others, t = 1, and (0m. v). < M. Collecting all data across

---

1 I S, we obtain S(M + 1) pairs of (t, o) classification examples. A binary classifier is then trained to maximize the conditional log probability of t given ø. If inference q were exact, no useful classification would be possible. In that case, the expected test log predictive density could be no higher than the negative binary entropy h((=)) == w log w + (1 - w) log(1 w) of Bernoulli distribution with parameter == 1/(M + 1).
Now imagine drawing a validation set in the same way, and evaluating the log predictive density of the learned classifier. We will show below (Thm. 1) that the expected log predicted density [10] ELPD =
log Pr(t)o) of the classifier on validation data is lower bound to divergence D1 between p p(fly) and q(flv) up to the known constant h(w),
ELPD1 - h(w) = (i q) = wKL (p(0ly) r(0lv)) + (1 K) (q(0l)) |l r(0l)))
(I)
where = 1/(MM + 1) and r(0ly)- = up(00y) + (1- (0lg) is a mixture of posterior density. If the classifier is optimal (c(tlo) Pr(t]o) in the distribution) then the bound in the above equation is tight maXelassifiers ELPD; - h(w) = (p(p.)). Here (L(p(Olv)lla(0)) (g) denotes a standard conditional Kullback-Leibler divergence'. By optimizing the classifier, max ELPD; h(w) becomes a commutable divergence. and its approaching zero is a necessary and sufficient condition for perfect calibration since Pp.p.) = if and only p(0ly) = {Blly) almost everywhere.
Example 2: Binary classification without y. Similar to Example 1, from each 1
simulation draw we generate M pairs of (1) examples, except that the feature á only contains the parameters 0, not y. A binary classifier is then trained to predict & given øb The ELPD of this classifier on validation data is lower bound to a generalized divergence D2 between the prior p(0) and q(0), up to a known constant lww)
0
1
01
1
On
ELPD2 - hww) < D2(p, 9) wKL (p(0) II r())) (1 - KL (q(0) r()))
(2)
where w = (M+1)) r(0) = wp(0) - (1 w (0) is the prior mixture, and the bound is tight when the classifier is optimal. A large EL.D reveals the difference between the inference q(0ly and p(Oly).
But D2 is only a generalized divergence: D2- = 0 is necessary not sufficient for ((fly) = p(0
Example 3: Binary classification with ranks (where the classical
Output: M + examples:
SBC is a special case). Instead of classification using full (0,y), y). we construct feature; the rank statistics. From each simulation draw we generate M +1 1 pairs of (t, 6) examples. The first pair is t = 0 and é = >m 1 < 0,)), the rank statistics of the prior draw. In the others. = 1 and Φ = E 100m.
Φ > 1 S 6m)
M 1(0, S 0m))
DM Dm 110s mm))
dml) + 10mm < 8). 1 15 < .M are the rank statistics of the inferred samples. A binary classifier is then trained predict given ø. The ELPD of this classifier is a lower bound a generalized divergence D3 between p(Oly) and o(Oly) up to a known constant h(w)
ELPD3 h(w)= S D2(p.q)= = (z q) || Uniform(0, 1)), w = 1/(M + 1),
(3)
and again the bound is tight if the classifier is optimal. Here Z(p, q) is a random variable defined by Z ())))), (0,y) ~ p(0,2)) v). where Q is the cumulative distribution function of g(0ly).
Training a "generative" classifier on this rank-based babel-gcnerating map is similar to testing for uniformity in rank statistics, done in traditional SBC which estimates the distribution of rlt 0 by histograms (See Appendix A. for precise correspondence between SBC and the naive Bayes classifier). The success of SBC suggests the usefulness of ranks, which motivates us include ranks or more generally feature engineering in the classifier. However, D3 is only generalized divergence: D3 = 0 is necessary but not sufficient for p(0ly) ((0l))). If inference always returns the prior, ((0))) 9) = p(0), then Dz(p, q) = O, a known counterexample of when rank-based SBC fails [32].
Example 4: Multi-class classification. We go beyond binary
Output:
M + examples:
labeling. Given simulation run (y,0, 0,,00,,.... "a,,,,.......... .02)). we now create
an (M+1)-class classification dataset with M + 1 pairs of (t, Φ) examples. In each one, the features are ò (y. 8*) where 8* is a
(0,01, 02 Bu,v)
(0,, 0, O2 .z... 00..))
permutation of (0, O,,... Ōst) that moves 0 into a given location and /. E 0, ..., M indicates the location of 0 in the permutation
M
(0,, 02 ... Bas 0.0)
Standard notation for conditional divergences [5] is that (L(p(()))ll((ll)) = Bptr,e) log p(0 x) Condi- tional divergence is not the divergence of conditional distributions. We interpret p(y) = q(v).
((()))
3

---

(See the table on the right). We train a (M+1)-class classifier to predict from 6. The ELPD on validation set is lower bound to the following divergence D4 between p(0ly) and q(Oly) up to known constant:
M HogDD++log(MJ)) < D4(p.2) = KL ((0) II ((0)),
M Σ mm0
p(0m)
II ((xx)) Ref
(4)
M +1
Again the bound is tight if the classifier is optimal. The divergence S D1 S log(M + 1), and D4 0 if and only if p(0lv) y) - q(0l)) l), necessary and sufficient for calibration. In Theorem 3 shows that as M + 0c, Da(p. q) converges to (p((()))). q(0ly)) at an (11/M) convergence rate.
3 Theory on calibration divergence
To generalize the previous examples, we define a "label mapping" + (y, 1y,0, 0,0 ...., ,,r))
((t,,o))..... (tz, ΦL)}. that maps one simulation run (y, 0, .............................. 0as) into a K-class classification example-set containing L. pairs of labels and features Φ. The label t {0, 1........ K 1) is deterministic and only depends on I. The features di can depend on I and (v. 0, 0 0m)). We only consider & satisfying that, when p(Oly) = ((lly), & given y is conditionally independent of t. (equivalently, (y. t) has the same distribution for any t). Let F be the set of these mappings +
We train a classifier on the classification examples created by & collected across repeated sampling. The classifier performance is measured by its expected log predictive density (ELPD), or equivalently the negative cross-entropy loss. Given p. y. and the mapping O let c(() be any probabilistic classifier that uses & to predict the label I, and Pr(kipo, c) is the predicted k-th class probability. Taking expectations over features and labels with (0, y) ~ p and Om ~ ((0lg) reaches the ELPD of the classifier,
,,,)) ELPD(0,c) c) == E,, log Pr(t = k)) (t, ogy 0 ........... We then define the prediction ability D(p, $ c) to be the ELPD plus the entropy of a categorical distribution, i.c.
(5)
K-
D(p q. 2,)) ELPD(0,C)
Σ
wk log wks where we Σ 1(t k).
(6)
The optimal classifier is the that achieves the highest prediction ability in the population: Dopt (p. q, T) = max (pp q. & c), where C is the set of all probabilistic classifiers.
(7)
eEC
The next theorem is the basic theory of our method: as long as we pass the simulation draws to a label mapping & IF, and train classifier the classification examples, then Dopt 4, +) is generalized divergence between p and 4.
Theorem 1 (Prediction ability yields divergence). Given any p. 9, and feature mapping & F, the optimal prediction ability Dopt (p. q, ②) is a generalized divergence from to q in the sense that Dept (p, q, 2) 2 O, and p(0ly) q(0ly) almost everywhere implies D (p, q, 4) 0. This generalized divergence is reparametrization invariant and uniformly bounded. For any classifier e,
K-1
$ D(p.q. S Dopt(p.
(p. 1) <-
Σ
Wk log wis
(8)
K=0
v(0 ac (00 Dopt (p. q, b) 0.
That is, any label mapping & produces generalized divergence Dopt q, 4), the prediction ability of the optimal classifier. The prediction ability (pp, q, 2, e) of any classifier c estimated on validation data is a lower bound to Dept (p, D). Further, this generalized divergence D (p, 4) is always proper Jensen-Shannon
Shannon divergence in the projected feature-label space (Theorem in the Appendix).
When p(0 y) * q(fly), to increase the statistical power, we wish that the generalized divergence can be as "strong" as possible such that we can detect the miscalibration. In the four examples in Section 2, we have used D,, D2, D3, D4 to denote the (generalized) divergence they yield. The next theorem shows there a deterministic domination order among these four metrics. Moreover, D4 the largest possible classification divergence from any given simulation table.

---

Theorem 2 (Strongest divergence). For any given and any 4 E. (1) (. D2. (2) Da and D are proper divergences. They attain if and only if = ((1l)) almost everywhere. They attain the corresponding upper bound in (8) if and only p(Oly) are q(Oly) are disjoint, i.e. SA (())))(()))) ( for any measurable set and almost surely y- (3) For any p.q q and a (I can have an arbitrary example size L and class size K), Da(p. 2 (p. q, Φ).
The following result shows that the divergence D4 in Eq. (4) approaches the "mode- spanning" KL divergence in the limit that the number of posterior draws M + co. This is appealing because for many inference routes, increasing the number of posterior draws is cheap. Thus, Da provides an accessible approximation to the KL divergence that is otherwise difficult to compute from samples.
Theorem 3 (Big M limit and rate). For any p.q. generate the simulation { 0, (a....,,11 and train the sitticclas-classier then Di(p, q) - KL(p(0p) II q(0lv)) + 0, as M - oc.
If further p(0 and q00) y) have the same support, and if Bp(olv)
p(Plu) <(0
< 00 for a.s. y. then
1
Di(p.q)= kL(p(Olv) (01
x (2(0))) p(0l)) +(
2M
where (-1l)) is the conditional chi-squared divergence.
4 Practical implementation
Algorithm 1: Proposed method: Discriminative calibration
input The ability to sample from p(0, y) and q(0 y). and a label mapping +.
output:() (i) estimate of a divergence between p(0 y) and q(0 y): (ii) p-value for testing p(0ly) = q(0 v)-
for ( = 11:55 S)
Sample (0., y.) ~ p(0,2) y). and sample 0., 02.00.............. 1.22.......,,, deM ~ (())):):
b simulation table
Generate a batch
examples of (t) = ((y,, o., la..... Penr), ). ssts 1; & label mapping
Randomly split the LS classification examples (t, Φ) into training and validation sets (all L examples for a
given go to either training or validation);
Train a K-class classifier to predict label f on the training examples, incorporating useful features;
Compute the validation log predictive density .. in (9), obtain an estimate of the divergence (7) and its hootstrap confidence intervals;
divergence
for (6- = 1: B)
Randomly permute the label in the validation set within each batch;
Compute LPD on the permutated validation set;
Compute the calibration p-value = 1/B Em (LPD) 2 LPD
b. frequentist test
Workflow for divergence estimate. We repeat the simulation of (y, 0, 0,0,,... Ō2r) ) for S times. Each time we sample (y. 0) ~ p(0, y) and 0 ~ g(lly) and generate a batch of L examples through label mapping E IF. In total, we obtain ST pairs of (t, Φ) classification examples. We recommend using the binary and multiclass labeling schemes (Examples 1 and 4, where L -= M +1 and K = 2 or M + 1) such that we can obtain a proper divergence estimate. We split the classification example-set {(t, })} into the training and validation set (do not split batches) and train a K K-class classifier on the training set to minimize cross-entropy. Denote Ivai to be the validation index, Pr(l = L 1g e) to be the learned class probability for any validation example (4g, 6z), we compute the ELPD (5) by the validation set log predictive density:
ELPD( c) ~ L.D val (t.c) = Ivail" Σ JiETval
log Pr(t = t))
(9)
( -0.c)-ckke -K- wy log Wi becomes lower bound estimate of the (generalized) divergence Dept (p. q, 1) in Thm. 1, and an estimate of D (p. 9. b) itself when the classier c is good enough. In addition to the point estimate of the divergence, we can compute the confidence interval. It is straightforward to obtain the standard error of the sample mean (9). To take into account the potentially heavy tail of the log predictive densities, we can also use Bayesian bootstrap [37] that reweights the sum in (9) by uniform Dirichlet weight.
For any c, LPD val
5

---

Hypothesis testing, Our approach facilitates rigorous frequentist hypothesis testing. The null hypothesis is that the approximate inference matches the exact posterior, i.e. p(0 = q y) almost everywhere. We adopt the permutation test: We train the classifier once on the training set and keep it fixed, and evaluate the validation set log predictive density LPD ( c) in (9). Next, permutate the validation set B times: at time b, keep the features unchanged and randomly permutate the validation labels within each batch of examples (&) generates batch of L examples each time), and reevaluate the validation set log predictive density (9) on permutated labels, call it LPDg.. Then we compute the one-sided permutation p-value as p _b-1 1(LPD 2 LPD //B. Given significance level, say 0.05, we will reject the null if p < 0.05 and conclude a miscalibration.
(f.c)
Theorem 4 (Finite sample frequentist test). For any finite simulation size S and posterior draw size M, and any classifier under the null hypothesis p(0ly) = g(Oly) almost everywhere, the permutation test is valid as the p-value computed above is uniformly distributed on [0, 1].
Our test is exact when the simulation size S and M is finite, while the original SBC [4] relied on asymptotic approximation. Further, we learn the test statistic via the classifier, and our test is valid regardless of the dimension of ( and there is need to worry about post-learning testing [3], while the usual SBC rank test will suffer from low power due to multiple testing.
Our test is always valid even if the classifier is not optimal. Why does our learning step help? For notional brevity, here we only reason for the binary classification. For any p, q, we apply the binary classification as described in Example 1, = 0 or 1 and Ø = (0, y)-
Theorem 5 (Sufficiency). Let e(0, y) = Pr(t = 1 y) be the probability of label I in the optimal classifier as per (7), and let no. and =² be the mimenimenio distributions of this o(0. y) when (0, y) is sampled from p(0, y) or from p(y)
respectively, then (i) Conditional on the summary statistic ē, the label is independent of features c = (0, (ii) Under regularity conditions, there is no loss of information in divergence as the joint divergence is the same as the projected divergence in the one-dimensional ê-space D1(p,q) (T[
That is, the best prediction è entails the best one-dimensional summary statistics of the high dimen- sional 8 x y space. The enhancement of the test power from using the sufficient statistics is then assured by the Neyman Peanaonn lemma [27].
4.1 Feature engineering: use rank, log density and symmetry
Reassuringly, whatever the classifier c, the prediction ability D(p, q, +,c) is always lower bound to the corresponding divergence D (p. q, b). and the permutation test is always exact. But we still wish to train "good" classifier in terms of its out-of-sample performance, for tighter bound in divergence, and a higher power in testing. We will typically use a flexible parametric family such as a multilayer perceptron (MLP) network to train the classifier.
The oracle optimal probabilistic classifier is the true label probability, and in the binary and multiclass classification (Example and 4), the oracle has closed-fore form expressions, although we cannot evaluate:
p(0ly)
p(0,, u) (0.lv)
Pr (L = 0) = p(0lv) + () binary
Pr (410,. Oas,u) = ΣM Ex p(0,. v) (0/ w) multi
(10)
Slatisticallymmannnnee feature. Depending on the inference task, we have more information than just the sample points and should use them in the classifier. In light of the shape and component of the optimal classifiers (10), the following maatistialll-mmmninnttii features are useful whenever available: (i) The log target density log (og(()l). As proxies, the log joint density logp((, y) and the log likelihood log p(y]() are often known to traditional Bayesian models. (ii) The log approximate density log ((lly), known to volitional and normalizing flows. (iii) When the log approximate density is unknown, transformation is to integrate the density and obtain the posterior CDF, Q(0ly) = Exng(zly) 1 (x > 0), and this CDF can be approximated by the rank statistics among the approximate draws up to a rescaling r0. =
Em=1,i=1 1 (w; 1) < 0..). See Appendix Table 1 for the usage and extension of these features in binary and multiclass classifiers.
MS
Linear features. We call the log likelihood, the log approximate density, or log prior (whenever available) linear features, and denote them to be ((0,, y). For example if both likelihood and q density

---

is known, then 1(0, v) = (log ( log (8ll))). Because they appear in the oracle 10, we will keep linear features in the last layer of the network (followed by a softmax). We recommend parameterizing the binary classification in the following form:
Pr(t (0. = inv Jov_Logt(MMLP(002 + w (0 v)].
(11)
Symmetry in multiclass classifier. The multi-class classifier is harder to train. Luckily, we can use the symmetry of the oracle classifier (10): the probability of class & is proportional to a function of (0k, y) only, hence we recommend parameterizing the multiclass probability
exp(g(Ok. y))
Pr(t- k0( 91. ........ AM..y) v))
g(0. y) MLP((. v) + w=((0.g) y). (12)
k== exp(g(Ow, v))
M
where 1(0, y) is available linear features. We only need to learn a reduced function from 9 x y to R. instead of from g x to E, reducing the complexity while still keeping the oracle (10) attainable.
Zero waste calibration for MCMC. If f d ~ ((lly) are produced by MCMC sampler, typically 0 has autocorrelations. Although classical SBC originated from MCMC application, the rank test requires independent samples, so SBC can only handle autocorrelation by thinning: to subsample 6m and hope the thinned 0 draws are independent. Thinning a waste of draws, inefficient when the simulations are expensive, and the thinned samples are never truly independent. With MCMC draws 6, our method could adopt thinning as our Thm. I and 4 are valid even when M = 1.
Yet do better by using all draws. The "separable" network architecture (12) is ready to use for MCMC samples. For example, we sample (0. v) p(0. y), and sample (f1...., 0s)) from MCMC sampler whose stationary distribution we believe is g()lv). and generate examples from the multiclass permutation (Example 4). Then we run separable classifier(12) to predict l. Intuitively, the separable network design (12) avoids the interaction between 0m with 6m,, and disallows the network predict based on the autocorrelation or clustering of Ó. The next theorem states the validity of our method in MCMC settings without thinning.
Theorem 6 (MCMC). Suppose we sample (0,y) v) ~ p(0, y). and sample (ō1...., .as) from
MCMC sampler whose stationary distribution we believe is q(0ly) (i.e. marginally 02 is from g(Oly)), and generate examples ((t,, ((1,,1)).- di)..... (1M+1: oM++)) from the multiclass permutation, such that & (0o, P1.....001 Then we run an exchangeable classifier (12) in which g is any θ x y - TR
mapping. Denote D DMCMMCsss (p. q) to be the predictive ability of the optimal classifier among all separable classifiers (12), then D (p..))opppp))))pp (p.4) Da(p.q)-
5CMC,spp
Dimension reduction and nuisance parameter. Sometimes we only care about the sampling quality of one or few key dimensions of the parameter space, then we only need to restrict the classifier to use these targeted dimensions, as result of Theorem 1. For example, in binary classification, we reduce the feature ó = (0. y) to Φ (h(0). y) in the classifier, where h(0) can be a subset of 0 dimensions, then the resulting classification divergence becomes projected divergence between h(())y. A ~ p(fly) and h(0)ly, 0 ~ g(0) v)) and other nuisance parameters do not impact the diagnostics.
Weighing for imbalanced binary classification. When M is big, the binary classification (Example 1) can suffer from imbalanced labels, and the divergence in (1) degenerates: D2 (p, q) - 0 as M + co. One solution is use the multiclass classification which has a meaningful limit (Thm. 3). Another solution is to reweigh the loss function or log predictive density by the label t. If the weights of class and class 0 examples are 2M M+1 and M+1, used in both training and validation log prediction density, then regardless of M, the weighted classification is equivalent to balanced classification, and -hasennn (JS) divergence KL[p(() y) ()(()())))) + ([)(())))))r((())) where r(0 g) = (p(Oly) + q({ly)]. See Appendix A for proofs.
the resulting divergence is the symmetric
Experiments
Closed-form example. Consider a multivariate normal parameter prior 0 E Rd MVN(O, Tda) and normal data model y ~ MVN(0, E), the exact posterior p(0ly) is explicit. In the experiments,
7

---

Our purp-values: vald under null
SBC wio multiple test ad).
Our has higher test power leasier detect miscalibration
1.0
10
ouTs
test restpower
6-3
05-
SBC
0.5
0.5
0.0
0.8
Infiation of post. cov
1.0
32
0.0
0.21
bias of post. mean (log scalel
0.02
0.25
61
Figure 2: Our test is valid as it yields uniformly distributed p-values under the mull p(oly) g(0lv). We check all dimensions at once.
Figure 3: Our test has a uniformly higher power than SBC rank test. We simulate wrong inference by multiplying the posterior covariance adding biases to the posterior mean.
Binery divergence =|5D
Multclasss div. converges to KL
SLCP
Goussian linear
Gaussian mixture
esiimeted divergence
es-
Trre L50 Diss
True Bbs
estimated divergence
ISe Ciess 1
2
0.0
S sample
Tree dv. bls =1
-0.5
10
100
De Duncozed to D1
regatii est All
1000
10
100
1000
iteratioon
tteration
Ite iterations
(# of simulatons)
M of post. draw per sim.)
Figure 5: We apply binary classifier calibration to poste- rior inferences three models from the SBI benchmark, and check the sampling quality after some iterations. The x-axis is the number of MCMC iterations. The y-axis the estimated JS divergence (+: standard error) between the tree target and sampled distribution at given num- ber of iterations, indicating gradual convergence.
Figure 4: In the Goussian example with bias = I or 2, with moderately big simulation sample size S. the divergence estimated from a binary (left panel) or multiclass (right panel) classifier matches its theory quantity dashed)): the Jensen-Shannon or Kullback Leibler divergence in the big-M M limit.
neural nets train binary (11) or symmetric multiclass classifiers (12), which generally recommend. We assume the inference posterior log ((fly) is known and added to the classifier.
Legitimacy of testing. First, to validate our hypothesis test under the null, we use true inference ((0/)) - p(0ly). With 16, we simulate posterior 6m ~ p(0ly). and run our -lassy-cl and obtain a permutation p-value. We repeat this testing procedure 1000 times to obtain the distribution of the p-values under the null, which is uniform as shown in Fig. 2, in agreement with Thm. 4. Because 0 has d dimensions, a SBC rank test on each margin separately is invalid without adjustment, and would require Bonferroni corrections.
500 draws (0, y) ~ p(0, y), then generate true
Power. We consider two sampling corruptions: (i) bias, where we add scalar noise (0.01 to 0.2) to each dimension of the true posterior mean, (ii) variance, where we inflate the true posterior covariance matrix by multiplying scalar factor (0.8 to 1.2). We compare our discriminative tests to cenferroni-correcte SBC chi-squared rank-test. test. In both settings, we fix a 5% type-I error and compute the power from 1000 repeated tests. Fig. 3 shows that our method has uniformly higher power than SBC, sometimes as good as SBC with a 10 times bigger simulation sample size.
Divergence estimate. The left panel of Fig. 4 validates Thm. 1: We run a weighted binary classifier on the Gaussian simulations with a scalar bias or added to the posterior mean in z. As the number of simulations S grows, the learned classifier quickly approaches the optimal, and the prediction ability matches the theory truth (dashed line): the Jensen-Shannon divergence between p(Oly) and q(0 y). The right panel validates Thm. 3: With fixed simulation size S we increase M, the number of posterior draws from g, and run a multiclass classification. When M is big, the estimated divergence converges from below to the dashed line, which is the theory limit (L((OO))) (00l))) in Thm. 3.
Benchmark examples. Next, we apply our calibration to three models from the SBI benchmark [23]: the simple likelihood complex posterior (SLCP), the Gaussian linear, and the Gaussian mixture model. In each dataset, we run adaptive No-U-Turn sampler (NUTS) and check the quality of the sampled distribution after fixed number of iterations, varying from 2 to 2000 (we use equal number of iterations for warm-up and for sampling, and the warm-up samples were thrown away). At each given MCMC iterations, we run our classifier calibration, and estimate the JS divergence, reported
8

---

under-confieet
caliorated over-confident
conc: under-
log density:
Jundy unde wearffentt
clossifer praifferrrdicto
of
Oay (baryon denss
concentration para.
log ql
Figure 6: Scatter plats of the classifier prediction v.s. Figure 7: In the cosmology example visually check one-dimensional parameter need to check can the classifeer predicted log odds v.s. cosmology param- reveal marginal over- or under-confidence inference eter Qh and cone, and log ()) Underconfidence is q(0 from a U or inverted-U shape.
detected in cone marginally, and more so in the joint.
in Fig. 5. In all three panels, we are able detect the inference flaws early iterations and observe a gradual convergence to zero.
Visual check. Our diagnostic outputs rigorous numerical estimates, but it also facilities visual checks. We make a scatter plot of the binary classifier prediction Pr(t = 1 a proxy of p(0 y) (0ly), against any dieennenniooal parameter we need to check: If that parameter is under- or over-confident, this scatter plot will display U- or inverted-U-shaped trend. Compared with SBC rank histograms, our visualization can further check the magnitude of mismatch (how far away Pr(( 1/Ø) is from 0.5), tail behavior (small ((]ly),, and several dimensions jointly. Fig is a visual check in the Gaussian example when we multiply the posterior covariance inference q by 0.8, 1, or 1.2.
Galaxy clustering. We would like to model galaxy spatial clus- I-flow
model 1, dy=38
model 2, d,=151 model 3, dy=1031
0.015 (0.004)
0.076 (0.007)
0.068 (0.007) 0.035 (0.006)
tering observations a Lambda 5-ens.
0.013 (0.004)
0.044 (0.006)
cold dark matter framework [34], where observations y correspond to statistical descriptors of galaxy clustering, and a 14-dimensional parameter 0 encodes key cosmological information about the Uni- verse. The forward simulation y/0 involves expensive N-body simulations (each single y simulation takes around 5000 cpu hours). We consider three different physical models pe they contain the same cosmology parameter 0 but three types of observations whose dimension dy varies from 38 to 1031, reflecting three power-spectrum designs. We apply simulation based inference to sample from each of the models, where we try various normalizing flow architectures and return either the best architecture or the ensemble of five architectures. We now apply our discriminative calibration to these six approximate inferences using a weighted binary classifier trained from S 2000, A 500 simulations for each inference (that is million examples). We have added the log densities log g(ll)) as extra features in the classifier, as they are known in the normalizing-floww flows. The table above -haennnn distances and standard deviation. Compared this table with the
shows the estimated estimates from blackbox MIP classifier (Appendix Fig. 8), using statistical features greatly improves the label prediction and tightens the divergence estimate. From the table above, all inferences y have a significantly non-zero but small divergence, among which the 5-ensemble that uses a mixture of flows always has lower divergence log odds (Fig. 7) reveals that q is under-confident in the parameter that encodes the galaxy concentra- tion rate. To interpret the tail behavior, the predicted log odds are negatively correlated with log joint normalizing
is more trustworthy. Further visualization of the classification
density (0ly), suggesting g00ly) > p(0ly) in the tail of (Oly)) another evidence of q underconfidence.
We share Jax implementation of our binary and multiclass classifier calibration in Github?.
Discussions
This paper develops classifier- diagnostic tool for Bayesian computation, applicable to MCMC, variational and simulation based inference, or even their ensembles. We learn test statistics from data using classification. Through flexible label mapping, our method yields a (family of) computable divergence metric and an always valid testing p-value, and the statistical power is typically higher features are available depending
than the traditional rank based SBC. Various on the task and we include them in the classifier. Visual exploration is also supported.
hhttps://githin htttp:://ggtabb..c///aaa/a/// http:////oooo/o/s//clltttrreeeee

---

Related diagnosties. The early idea of simulation based calibration dates back to [13] who compared the prior and the data-averaged posterior, i.c., to check p(0) = (g(()))p(())d using simulations. [42] further developed techniques for the first and second moment estimate of the data-averaged posterior using the law of total variance. Our method includes such moment comparison as special cases: using the no-y binary labeling (Example 2), then comparing the empirical moments of the sample the sample can be achieved through a naive Bayes classifier or a linear discriminant analysis.
The rank-based SBC [4, 38] can be recovered by ours using binary labels and taking ranks to be features (Example 3). The rank statistic is central to SBC, which follows the classical frequentist the tail probability under the posited model quantifies the extremeness of the realized value[]1] is only convenient way to locate the observations in the reference distribution, especially in the past when high dimensional data analysis was not developed it is the use of modern learning tools that sharpens our diagnostie. Recently, SBC has developed various heuristies for designing cone-dimensional) test statistics 6(0, y). Such rank tests are recovered by our method by including the rank of & in our features (Sec. 4.1). For example, [26] proposed to test the rank of the likelihood p(0l)) in MCMC, [20] looked at the rank of proposal density q()) ly), [22] used the g-probabilitty integral transformation in normalizing flows. In light of our Theorem 5, the optimal test statistic is related to the density ratio and hence problem-specific, which is why our method includes all known useful features in the classifier and learns the test statistic from data.
Classifier two-sample test. Using classification to compare two-sample closeness is not new idea. The classifier two-sample test (C2ST) has inspired the generative adversarial network (GAN) [16] and conditional GAN [25]. [33] has developed GAN-typed inference tools for SBI. In the same vein, [19] used classifiers to diagnose multiple-run MCMC, and [24] used classifiers to learn the likelihood ratio for frequentist calibration, which is further related to using regression to learn the propensity score [36] in observational studies. The theory correspondence between binary classification loss and distribution divergence has been studied in [21, 28, 35]. This present paper not only applies this classifier-for-two- fo------ssapplee-eeet idea to amortized Bayesian inference to obtain a rigorous test, but also advances the classifier framework by developing the theory of the "label mapping". while the traditional GAN-type approach falls in the one-class cwass-per-greee per-group category and deploy binary classifiers as in Example 1. autocorrelated, and imbalanced.
Our extension is particularly useful when samples are overlapped (multiple @ per y),
KL divergence estimate from two samples. As a byproduct, our proposed multiclass-classife provides a consistent KL divergence estimate (Thm. 3) from two samples. Compared with existing two-sample KL estimate tools from the f-divergence representation [29] or the Donsker- Varadhan representation [2], our multiclass-classinee
estimate appears versatile for applies to samples with between-sample dependence or within-sample coro-correlation. It is plausible to apply our multiclass- classifier approach to other two-sample divergence estimation tasks such as the f-GAN [30], which we leave for future investigation.
Limitations and future directions. Like traditional SBC, our method assesses the difference between p(0ly) and q(O))) y). averaged over y. This is a "global" measure relevant to developing algorithms and statistical software. But sometimes the concern is how close p and q are for some particular observation
yobs. "Local" diagnostics have been developed for MCMC [14, 12, 40, 17]. variational inference [41, 8] and simulation- based inference [42, 22] that try to assess oBly)bb)) only. It would be valuable to extend our approach to address these cases. Another future direction would be to extend our approach to posterior predictive checks [11, 9, 39] that diagnose how well the statistical model p fits the observed data.
Acknowledgement
The authors thank Bruno Regaldo-Saint Blancard for his help on galaxy clustering data and simula- tions. The authors thank Andreas Buja, Andrew Gelman, and Bertrand Clarke for discussions.
References
[1] Agrawal, A., Sheldon, D. R., and Domke, J. (2020). Advances in black-box VI: Normalizing flows, importance weighting, and optimization. In Advances in Neural Information Processing Systems.
10

---

[2] Belghazi. M. I., Baratin, A.. Rajeswar, S., Ozair, S., Bengio. Y., Courville, A.. and Hjelm. R. D. (2018). MINE: mutual information neural estimation. In International Conference Machine Learning.
[3] Berk, R., Brown, L., Buja, A., Zhang, K., and Zhao, L. (2013). Valid post-selection inference. Annals of Statistics, 41(2))802-833 837.
[4] Cook, S. R., Gelman, A., and Rubin, D. B. (2006). Validation of software for Bayesian models using posterior quantiles. Journal of Computational and Graphical Statistics, 15(3):675 692.
[5] Cover, T. M. and Thomas, J. A. (1991). Elements of Information Theory. John Wiley & Sons, 2rd edition.
[6] Cranmer, K., Brehmer, J. and Louppe, G. (2020). The frontier simulation based inference. Proceedings of the National Academy of Sciences, 117(48):30055-3006
[7] Dax, M., Green, S. R., Gair, J., Macke, H., Buonanno, A., and Schölkopf, B. (2021), Real-time gravitational wave science with neural posterior estimation. Physical Review Letters, 127((4):241105
[8] Domke, (2021). easy to interpret diagnostic for approximate inference: Symmetric divergence over simulations.
[9] Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., and Gelman, A. (2019). Visualization in bayesian workflow. Journal of Royal Statistical Society: Series A, 182(2)389--022.
[10] Gelman, A.. Hwang, J., Vehtari, A. (2014). Understanding predictive information criteria Bayesian models. Statistics and Computing, 24(6):997-1016
[11] Gelman, A., Meng. X.-L., and Stem,, H. (1996). Posterior predictive assessment of model fitness realized discrepancies. Statistical Sinica, 6:733-760.
[12] Gelman, and Rubin, D. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4):457-472
[13] Geweke, (2004). Getting right: Joint distribution of posterior simulators. Journal of American Statistical Association, 99(467):799-80)
[14] Geyer, (1992). Practical Markov chain Monte Carlo. Statistical Science, 7(4):473-483.
[15] Gonçalves, P. J., Lucckmann, J.M.., Deistler, M., Nonnenmacher, M., Öcal, K., Bassetto, G., Chintaluri, C., Podlaski, W. F,, Haddad, A., and Vogels, T. P. (2020). Training deep neural density estimators identify mechanistic models of neural dynamics. Elife, 9:e56261.
[16] Goodfellow, L. Pouget-Abadie. Mirza, M.. Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio. Y. (2014). Generative adversarial networks. In Advances in Neural Information Prucessing Systems.
[17] Gorham, J. and Mackey, L. (2017). Measuring sample quality with kernels. In International Conference Machine Learning.
[18] Hermans, J., Delaunoy, A., Rozet, F., Wehenkel, A., Begy, V., and Louppe, G. (2022). A trust crisis in iimulation-basee inference? Your posterior approximations can be unfaithful. Transactions on Machine Learning Research.
[19] Lambert, B. and Vehtari, (2022). R": robust MCMC convergence diagnostic with uncertainty using decision tree classifiers. Bayesian Analysis, 17(2):553-379
[20] Lemos, Coogan, A., Hezaveh, Y., and Perreault-Levasseee L. (2023). Sampling-based accuracy testing of posterior estimators general inference. arXiv:2302.03026,
[21] Liese, and Vajda, I. (2006). On divergences and informations in statistics and information theory. IEEE Transactions on Information Theary, 52(10))43944411
[22] Linhart. J.. Gramfort, A., and Rodrigues, (2022). Validation diagnostics for SBI algorithms based on normalizing flows. arXiv:2211.09600
[23] Lueckmann, J.-M.. Boelts, J.. Greenberg, D.. Goncalves, P., and Macke, J. (2021). Benchmarking simulation based inference. In International Conference on Artificial Intelligence and Statistics.
[24] Masserano, L., Dorigo, T., Izbicki, R., Kuusela, M., and Lec, A. B. (2022). Simulation-base inference with WALDO: Perfectly calibrated confidence regions using any prediction or posterior estimation algorithm.
[25] Mirza, M. and Osindero. S. (2014). Conditional generative adversarial nets. 4rXiv:1411.184.

---

[26] Modrák, M., Moon, A. H., Kim, S., Bürkner. P., Huurre, N., Faltejsková, K., Gelman, A., and Vehtari.
A. (2022). Simulation-based calibration checking for Bayesian computation: The choice of test quantities shapes sensitivity. (exXiv2221100260
[27] Neyman, J. and Pearson, E. (1936). Sufficient statistics and uniformly most powerful tests of statistical hypotheses. In Statistical Research Memoirs, pages 1-37, Cambridge University Press.
[28] Nguyen, X., Wainwright, M. J., and Jordan, M. I. (2009). On surrogate loss functions and f-divergences. Annals of Statistics, 37(2);876 904.
[29] Nguyen, X., Wainwright, M. J., and Jordan, M. I. (2010). Estimating divergence functionals and the likelihood ratio by convex risk minimization. IEEE Transactions on Information Theory, 56(1 16584778877
[30] Nowozin, S., Cseke, B., and Tomioka, R. (2016). f-GAN: Training generative neural samplers using variational divergence minimization. Advances Neural Information Processing Systems.
[31] Papamakarios, G., Nalisnick, E., Rezende, D. Mohamed, S., and Lakshminarayanan, B. (2021). Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22( 22(1):2617-268.
[32] Prangle, D., Blum. M.. Popovic, G., and Sisson, (2014). Diagnostic tools of approximate Bayesian computation using the coverage property. Austrulian & New Zealand Journal of Statistics, 56(4):309-329.
[33] Ramesh, P., Lucckmann, J.-M., Boelts, J., Tejero-Cantero, Á., Greenberg, D. S., Gonçalves, P. J., and Macke, J. H. (2022). GATSBI: Generative adversarial training for simulation-based inference. International Conference on Learning Representations.
[34] Régaldo-Saint Blancard, B., Hahn, C., Ho, S., Hou, J., Lemos, P., Massara, E., Modi, C., Moradinezhad Dizgah, A., Parker, L.,, Yao, Y., and Eickenberg, M. (2023). SimBIG: Galaxy clustering analysis with the wavelet scattering transform. arXiv:2310.1522
[35] Reid, M. and Williamson, R. (2011). Information, divergence and risk for binary experiments. Journal of Machine Learning Research.
[36] Rosenbaum, for causal effects. Biometrika, 70(1):11555
and Rubin, D. (1983), The central role the propensity score in observational studies
[37] Rubin, D. B. (1981). The Bayesian bootstrap. Annals of Statistics, 91)::300-344
[38] Talts, S., Betancourt, M., Simpson, D., Vehtari, A., and Gelman, A. (2018). Validating Bayesian inference algorithms with simulation-basee calibration. arXiv:1804.0m 06788..
[39] Vandegar, M., Kagan, M., Wehenkel, A., and Louppe, G. (2021). Neural empirical Bayes: Source distribution estimation and its applications to iimulation-based inference. In International Conference on Artificial Intelligence and Statistics.
[40] Vehtari, A., Gelman, A., Simpson, D.,, Carpenter, B., and Bürkner, P.-.. (2021). Cank-normalization. folding, and localization: An improved ft for assessing convergence of MCMC. Bayesian Analysis, 16(2):667- 718.
[41] Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018). Yes, but did it work?: Evaluating variational inference. In International Conference on Machine Learning.
[42] Yu, X., Nott, D.J., J. Tran, M.-N., and Klein, N. (2021). Assessment and adjustment of approximate inference algorithms using the law of total variance. Journal of Computational and Graphical Statistics, 30(4))977-990.
12

---

Appendices to "Discriminative Calibration"
The Appendices are organized as follows. In Section A, we provide several extra theory results for the main paper. In Section B, we discuss our experiment and implementations. In Section C, we provide proofs for the theory claims.
A Additional theory results
Theorem 7 (Closed form expression of the divergence). Consider the label generating process in Section 4. Let (t, Φ) be the joint density of the label and features after mapping i.e. the distribution of the (t, 6) generated by the simulation process:
1. sample (6, y) p(0, y):
2. sample (01, (01....,axl 0Ms) from q(llv);
3. generate ((11,01).... (tz, (1))) = (( 0,01,
4. sample an index from Uniform(1.......
5. return o1).
We define
T& (T) T k)= (0) / (
k)dy be the y y-averaged law of ølt k. Note that y has been averaged out.
Then classification divergence defended through (7) in Theorem has a closed expression, a divergence in this projected o space.
K Σ k=0
K ll
Dopt (p. 9, D)
WL
Tx(+)
Σ wynz()) j=0
(13)
This Theorem 7 gives a connection between our divergence and the familiar Maasen-Shannon diver- gence, which is well known be linked IID sample classification error (note that the classification examples we generate are not IID for they share y).
As special case, when K Jensen-Shannon divergence, or JSD. In general, between two densities Tı (Ir and Ir (2)... the JSD is
Wn wi 1/2, recover the symmetric two-group
Jensen Shannon divergence( m2) KL
1 2
m
T1
n2
+
KL
n22
41 + 22
(14)
2
As before, the standard notation [5] of conditional divergence is defined as taking expectations over conditional variables:
Ti (x|22) T2(2122)
(xl) (x(l))
KL
(x(=2) ||
2
1 + KL
[x|2)) + (2(22)
T2(2|2)
||
(15)
2
2
Theorem 8 (weighting). The binary label mapping generates M paris label-1 examples and one
pair of label-0 examples per simulation. We can reweight the binary classifier by letting the ELPD
be E[ 1 1) logp(t 1[c(o)) normalizing constant.
CMM 1(t -M
0) log 0fc(Φ))], where C²= (M+1)² is
-
2M
That is if we modify the training utility function the validation data LPD to be:
C
1(4 = log Pr(( 1[c((.)) +
CM
1(4 = 0) log Pr(t 0lc(o)))
M 11
M +
13

---

then the resulting binary classification divergence in Thm. 1, i.e.
weighted ELPD + log2,
is the conditional Jensen Shannon divergence (15).
A.1 SBC test and generative classifier
In Example 3, stated the intuitive connection between the SBC rank test and using rank as the feature a "generative" classifier (in contrast to our "discriminate" approach). Now we make such comparison more precise: Assuming the parameter is one-dimensional. This is what SBC does: From p we obtain P.y. and for each 0 we compute the rank it the paired @ samples who shared the same y: r 2m=11 1 < ). Repeating this simulation many times, we obtain sample of and SBC will test if such is discrete- unifoomlyy distributed. This test could be done by simulating reference distribution that matches the null. We do so by generating the ranks in the samples: in -nnn (0,, I fm.) (0,, I f). To be clear, there are other ways to generate reference samples, but all it matters is that the reference distribution matches the distribution of r under the null. The uniformity test in SBC now becomes to test and have the same distribution. For example, we can do two-sample Kolmogorov-Smiroo test or permutation test.
In comparison, this is what naive Bayes would do: the naive Bayes classifier first estimates the two conditional distributions: Pr(r|(= 0) from the empirical distribution of the prior rank r, and Pr(r]t = 1) from the empirical distribution of Then test if Pr(rlt 0) Pr(rit 1) becomes the same to testing if the empirical distribution of r is the same as the empirical distribution of r-the same task as SBC, The point is that given any two-sample based test,
SBC + sample-based test
is operationally the same as only use rank + naive Bayes classifier + density estimated by histogram + sample-based test.
B Implementations
B.1 A cheat sheet of useful pre-learned features
Table summarizes how to use the statistical feature (whenever available) in both the binary and multiclass classifiers.
binary (0,z)
multiclass
full feature MCMC VI likelihood-free
(0....., 0K. p(0,. y) y), g( K p(Or. y) g(0, K p(O:), q
p00 y), r(0, y) p(0, y), ((0l)) p(0), (0
Table 1: A cheat sheet of useful pre-learned features
B.2
Code and data
We share the python and Jax implementation of our binary and multiclass calibration code in https:// wtt.ab/.coomyoooyy//a/aaaaaa.aaa
In the MLP training we include a standard L weight decay (i.c., training loss function cross entropy loss + tuning weight 2 penalization). We tune the weight of the decay term by a fold cross-validation in the training set on a fixed grid (0.1, 0.01, 0.001, 0.0001}.
In the cosmology experiment in Section 5, we use mee-hidee--laye layer MLP with 64 nodes to parameter- ize the classifier with the form (11), with additional pre-learned features such as (0gg(0)) y)) added as linear features. The approximate log density log ((0ly) is known (i.e., we can evaluate the for any 0 and y, or at least least up to constant) in either the normalizing flow or the ensemble of the normalizing flows. One classification run with roughly one million examples took roughly two hour
14

---

model (d,=38)
model 2 (d,=151)
model 3 (d,==101))
1-flow
0.008 (0.003)
0.05(0.006)
0.0001 (0.0002)
5-ensemble
0.001 (0.002)
0.008(0.003)
0.001( 0.001)
Figure 8: The estimated divergence of two simulation based inferences (either using one flow or ensemble of five flows) three cosmology models. These three models have the same parameter space and involve expensive simulations. Here apply our classifier approach and estimate the Jensen-Shannon divergence (and standard deviation) using weighted binary classifiers from S 2000, M 500 simulations: that is 2000 501=1 million classification examples. In this table, we run this classifier by the black-box multilayer perceptron (MLP) on the full (0,y) vector, while the table have in the main paper further adds log g(0 y) as a pre-learned feature since it is known in the normalizing flows. In this black box estimate table, we would pass a l-test and not even detect misealculation. Adding a statistically menainggff feature greatly improves the classifier prediction, and tightens the bound of the divergence estimate.
cpu time on a local laptop. It would be more efficient to run the classification on GPU, thought the classification cost is overall negligible compared with the simulation step y|8 which is pre-computed and stored.
In the closed-form Gaussian experiment in Section 5, we consider the easiest setting: 0 R MVN(O, Ida) and a normal data model y ~ MVN(0, Ida). Indeed, we have kept both the true distribution and the sampling corruption mean-field to make the calibration task easier for the traditional SBC rank-test. The power of the traditional SBC rank-test would be even worse if we add interaction terms in the corrupted samples (0ly)) while our method leans the joint space sampling quality by default. The exact posterior p(0ly) is known but we pretend we cannot evaluate it. We set d 16 to be comparable with the real data example. The right panel validates Thm. fixed the simulation size S 5000 and vary M. We have tried other fixed M 1000 but seems the resulting graph is very similiar. The confidence band in Figure and Figure is computed using 1.96 times standard error from repeated experiments.
Table 8 shows the naive classifier: using the bandbox MLP without the linear features. The estimate rather loose compared with the table we obtained in Section which used statistically-meaninghhi features.
C Proofs of the theorems
We prove the theorems in the following order: Thm. 7, Thm. 1, Thm. 8, Thm. 6, Thm. 3, Thm. 5, Thm. 4, and finally Thm.
Contents
C.1 The classification divergence
15
C.2 Reweighting (Theorem 8) leads symmetric divergence
18
C33 MCMC without thinning (Theorem 6)
19
C.4 Large sample limit and rate as M - oc (Theorem
................................
20
C.5 Valid hypothesis testing (Theorem 4)
22
C.6 Why does classifier help test? (Theorem 5)
23
C.7 The maximum discriminative generator (Theorem 2)
C.1 The classification divergence
To prove our Theorem 1, we first state an equivalent but a more general theorem. Note that we are allowing non-IID classification examples for the shared variables in this theory.
15

---

Theorem 9
(The general relation between divergence and classification error), Define a simulation process (2, k)) through the following process, among the three random variables, E RR,, IR"" can have any dimension, and (1,2...., K] is an integer.
Sample z ~ T2))
2. Sample k Categorical(w), where simplex.
3. Sample r ~ TT<z)z))
4. Return (z, k, x)
Define U(T) to be the expected log probability or cross-entropy of an optimal classifier trained to predict k from (x, 2), ie..
U(x) max Ex U c(x,,)z
where Ukk, c) = log c(k). Then
K Wh log Wk
K Σ wy (x)()))) |l (<]2))
U(x)
k=1
where xx]22) Σk= k= (x|2) is the mixture over all groups.
We define D(x) U(T) E wi log wk. Then this has the following properties:
1. (lower bound). D(T) 2 O. It achieves zero if and only if all pı are the same, iee. Pr(X A) Pr(P E A) for any measurable set A.
2. (upper bound). D(T) S En=n we log wx. This maximum is achieved is are disjoint, iee. Pr(X* E A) Pr(X' A) = for any measurable set A.
3. (alternative upper bound when the difference small). maxz, (x
D(T) S
4. (reparametrization invariance). For any fixed bijective transformation + g(æ). let p be the law of g(x²), ), then Diy.....,, = D(pk........)
5. (transformation decreases divergence). For any fixed transformation x + g(x), let ph be the law of g(x²), then D(p < D(pk.....Ppy
Proof. The (unknown) optimal conditional classifier is the Bayes classifier:
Pr(t k| 2) TT (x)2)
(x(2)
Let r Σ5n (xTr( be mixture density marginalized over the group index k.
We emphasize that are using the notation of conditional KL divergence:
KL (x/2) m()]==)) = / p(z) /
[2(2)z) r(x|2) dr dz.
(=(=)) log
16

---

The expected log predictive density of the optimal classification is then ELPD(Pr(tt k//z, 2)) = Ex r(x22)
K
(
(xTT((1)2))) r(x|z)
dæ )
log
r(z|2)
E
Wkk
[<(()|2) log
r(x=22
log wk do
K
wk Es /
(x((rl2)
E,
wy log(
+
The xx(2²2 log
dx r(<z|2)
n=1
Wh log ww) +
wp KL (Tx())2), r(zl2)).
1m
K When ELPD Σx= wy log(ur) 0. Σx= T& (x]22) are the same almost everywhere.
WR (x)2) = 0, hence these K conditional densities
The reparametrization invariance is directly inherited from KL diverge.
The upper bound is the consequence of ELPD < 0 for category classification and the lower bound is that KL divergence 2 0.
Let's 's emphasize again that we are following the notation of conditional divergence in [5], such that for any two joint density p(thcta, y) and 98. y), the conditional KL divergence (pp())))))l((())) = Ep(v) Ep 0) log q(yloj. an expectation is taken over y.
p(a10)
The procedure Theorem can reflect many different generating processes. Let's recap the examples we show in Section 2.
In Example of the main paper, the binary case, K and we generate M samples from g. T1 = p(0 y). ((0l)), w1 = 1/M, w2 = (M 1)/M. Our simulation process in example is equivalent
1. Sample marginal p(y),
2. Sample Categorical(w), where w1 1/M, w2 = (M 1)/M.
3. Sample Ti. (x)2). where गi p(0ly), .T2 <(lly).
4. Return k.
As direct consequence of Theorem 9, the resulting divergence in the binary case Di(p.) 9) w1 KT. (p(Ol)) (00l)) we KT. (g(Olv) ((0l)))
where wip + W2g.
In Example 2: Binary label with y. generate M samples from K = 2, p()), π2 9(0).
Wi 1/(M + 1). w/2 = M/(M + 1). Our simulation process in example is equivalent to
1. Sample Categorical w), where w1 1/M, w2 (M 1)/M.
2. Sample ~ (x:), where म = p(f), 72 = 9()).
3. Return (k:, .))
From Theorem 9, the resulting divergence reads
Dz(p, 4) = w KL (p(0) r(0)) + w2KL (4(0) r()))
In the multivariate Example 4, K M 1, and the individual density is
F (01.k) = p(Oxlv) II
( 1v)

---

From Theorem 9, the resulting divergence reads
K Σ
Da(p, = KL
(p(Orla) II m=
g000 (v) ||
K p(Olv))
II (((mll))) mifk
give more results on this D4 in Section C.4.
The discriminative calibration is designed to detect and gauge the difference between p(0lg) and g(Oly). From now on we use long vector 01:M+1 (the index starts by to represent the simulation draw. Each simulation draw first generates one 0₃ from the prior, one data point from the sampling distribution ply 11). and finally M draws from the approximate inference M ((fly).
The "label generating process" d takes a simulated draw as its input and returns a vector of features and labels t.
: (v. M+1)) > {(01, ti),.., (Øx tx))-
The classifier will then see this feature and labels (the index & is disregarded) across all simulations.
We only ask that the map & be within in set IF: satisfying that:
F Gly is independent of if p(Oly) = g(0ly), a.e.
Proof of Theorem 1.
Proof. The classifier sees features and labels t. Given at, (é, y) IID draws from p(y)m(olv, In case y is not observable by the classifier, olt are IID draws from J, ply)p(oly,t)
To use Theorem 9, let x* {(1b) kb and wi = Pr(t k). The divergence from Theorem 9: reads
K
K
Dopt (0, ())
Σ
Wx KL
wx
(T)))
Σ
WE (0)
==1
1=1
in which [p( p(olt k) -(, essentially Theorem 7.
k)dy is the y-averaged 6 margin. This is
If p = q, from IF, the law of o, y is independent t (y is always independent t marginally), so that (Ø) = f,, ((ly), k)dy is also independent of l, as it only depends on y and ø, hence the divergence remains O.
C.2 Reweighting (Theorem 8) leads to symmetric divergence
We now prove Theorem The binary label mapping generates pairs label-1 examples and one pair of label-0 examples per simulation. We can reweight the binary classifier by letting the ELPD
CM 1( = 0) log p(t 0}c(ø))], where C (M+1)² be E M+1 1( 1) log p(t = 1[c(o)) + M+1 normalizing constant. We want to prove that after weighting, the classification divergence symmetric Jensen shannon distance, as if we have balanced classification data.
C
2M
Proof. After the label-reweighing, the optimal classifier is the balanced-class classifier
Pr(p]/
ctt 1 Pr = + Pr(p]t 0)
The weighttd ELPD in the population is
CM
ELPD [r(P 0/0)
log c(t = 0ld) + Pr(t= 0 010)
log c(t 1[0)]
M +1
M +
p(0l))) + Mg
p(0lz)
CM
log p(0 y)
Mq(0lg)
C
(00g))
Ey
log
M +1
p(0 v) + Mq(0ly) M +
2r(0
p(0w) + Mq(0lp) M +
2r (g)
2M
p(0 y)
log p(0
((0 (g)
((1y)
B.y C
log 2r 2r(0
(M + 1)2
2
2r (0 y)
2
18

---

With an appropriate normalizing constant C
(M+1)²
2M
p(fly)
p(0ly) r(0ly)
ELPD By 2 p(((y)) log
2 ( log
log 2.
r(0l))
That is, the reweighted ELPD log 2 is the conditional Jensen Shannon divergence (15) between p(Oly) and g(fly).
C.3 MCMC without thinning (Theorem 6)
To diagnose MCMC, we sample v) ~ p(0, v). and sample (potentially autocorrelated) draws (fI...., 0m) from a MCMC sampler whose stationary distribution we believe is ((0l)) (i.e., marginally 0g is from g(fly)), and generate examples ((41, (1), (IM+1, (bM+)))+ from the multi- class permutation (see definition in Example 4), such that
(0o, .0..0.....)).. 0M).
Then we run an exchangeable classifier parameterized by
exp(g (fx. v)) k'= exp(g(Ow: v))
Pr(t = (00, 01, ..... Pas, y))
(16)
where g any θ + E mapping to be learned.
Note that here Pr(t = k: (0p. ...... 0M v)) is the classifier model we restrict it to be, not the true population In general (t = (6o, ...... AM. y)) = Pr(t (6,, . (M.y)) Roughly speaking, the optimal classifier should be the Bayes classifier projected to the restricted space (16), and we need to proof that this projected restricted solution turns out to be the same the IID case 10, which is not trivial in general.
Proof. Intuitively, this separable network design (16) avoids the interaction between Om with 0. and disallows the network to predict based on the autocorrelation or clustering of 0.
Because of the permutation design and we define q to be the marginal distribution 0 first, we know that
T(0xly, t = k:) = p(Oxl)), (0x y. * Is) = g(0ulv),
in the MCMC population π (there is no need to address the joint for now).
From this we obtain the conditionals the population
p(Oxlv)
π kl0z. v)
p(0xly) + Mq(0ly)
and further the ratios for any m, l withtin the index set {0, .0, M) and m + k:
p yy (xl
T(t = k: Ok. Om. y) p(0xlv) (v)) + (xxx yy ( 0z) yy) + (M 1) 0
Here (12(0k, On) ) is the joint distribution of out of M draws from the MCMC sampler (0ly)) which would often not be the same as ((0xlv)q(( chains. We do not need specify the form of this joint. The key blessing is that when 02 is from p(Oly) and 6. is from ((lly). then they are independent (the numerator).
because of the autocorrelation of Markov
Next, from the line above we obtain the ratio estimate in the true population
( k: O,, Om: y) π m 02, Pm. v) q(Oxlv)p(PAlb)
p(0n p(Oulw)g(Oalv)) p(0x ()/(()))))))
(17)
p(0m. ()/q(0mlv)))
Now we project the restricted classifier (16). Intuitively, the separable classifier "almost" only depends on Or, except for the normalizing constant. We can remove the dependence on the normalizing constant by specifying the ratio
Pr(1 = (00, OI,..., 0M, 9))
exp(g(0E, n))
Pr(t = m (Oo. 01...., 0M, v))
exp(g(0m. m, 8)
19

---

Marginalizing out all other components, we obtain
Pr(t (0pp 6m 8)) Pr(t = mf mm u))
exp(g(Ok.1) y)) cxp(g( (m v))
(18)
Matching the restriction (18) with the true population (17), the restricted projection is attainable and only
exp(o(Ot, 3) p(0,, ((t,,)//(0elvv)
so that the optimal classifier needs to be
p(O,. (0g.v))/((.1p))
Pr(t|0 0K. E v)/4(0xl))
happens that this MCMC restricted optimal classitier matches the IID optimal classifier (10). follows from the proof in Theorem 9 that the classification divergence using the restricted classifier is still (p,
(p.gg 4), as if (fm I are IID samples from g(0ly).
C.4 Large sample limit and rate as M + oc (Theorem 3)
In the multivariate Example 4, from each draw 8, d1.M- we generate M + examples from K = .M + classes; one label from each. We will use the index starting from 1 in this subsection. For the k-th example, 4 = k, and the permutation for classifier reads Φs a vector including y and K copies of 0, where the k-th copy is the prior draw, and the remaining 0 are from ((ly). Slightly abused the notation, in this subsection, we call this long feature vector as 01:k ; it is vector in the M++ space. We prove Theorem 3.
Proof. First, we write the true conditional label probability in this process:
p(0ely) 11
p(0.lu)/((0Al)))
(t|01.k. E. p(0r (v) + p(0, y) Σ, ((0,ll1)/4(00
Plug it as the classifier, we obtain the optimal ELPD or negative cross entropy: ELPD log (1)0
ELPD E p(Orly)/y(C.
Σ 1= p02 (O/)()
p(0,
E log
Elog Σ p(Oolu)/q(Ok).
q(0.lv)
The first term above is simply
E Elog p(0rly) q(,, y)
= KL (p(Ol)))l<(BW))
According to definition (4), the divergence is ELPD offset entropy term,
Da = ELPD + log K.
We now derive the limit of D4 - KI, (p (1)))) when M - xo (or equivalently K = M+1 00)
A = D4 KL (p(0)v))la(0)))
E.PP log K KT. (p(00)))))
p(Ox
log K E Elg
((<x v)
K
p(Oxlu) (v)
= - Eog
K
((Oxlv)
A
20

---

Given any label value 1 < K, 0, ~ p(-l)), and all the remaining fg ~ a(-lv) for j & t.
Let
K
p(Okly)
p(Orl))
p(0rly)
X K
((Oxlu)
K ((0ll))
Kep ((01l))
The first term
p(Oelu) ((eely)
- O.
K
The second term, call it A2 is the mean of IID means as 02 g(-l)) for # the law of large number yields
K
p(O.lu)
p (<<ly)
A =
A + (erg()))))
1.
K hett
(00ll))
This proves that
A 0, as K 00.
Hence,
D4 KL (((0)))))((0))) - 0,
which finished the first part of the proof.
Now let's derive the rate. Under regularity conditions, for example, the variance of density ratio is bounded, i.e., there exists constant C' co, such that for all y,
p((ele) (0, y)
Varo
<C,
then CLT holds, such that the second term above has normal limit,
VK(A2 - 1) normal(O, o)), in distribution,
where
0² Var. p(Ocl1)
((0rly)
(p(Oly)
2
-1) (p(8l))
((flg)
Ey (o)g((lp)
q(0 v)
=x (p(()) v) ((ly)))
which is the definition of the conditional chi-squared divergence.
Consider a Taylor series expansion, log(1 + z) express the log function and the Slutsky theorem to ignore the zero term, we get
f.² + o(2²). Using the Delta method to
K Eog 0² 2
Plug this in the definition of A we obtain the desired convergence rate,
D4 KL (p(0W))l((0w)) 2M x2 (p(0lv) ((0))) + o(1M))
This proves the claims in Theorem 3.
21

---

p-values under null
our permutation
naive permutation
0.5
0.5
Figure 9: Our designed permutation vs the naive permutation.
Table 2: Batch permutation of labels
Original batch from
a batch of permuted labels label
another batch of permutation
label t features b
label t
(0,
0
0
(01,9)
0
(02,9)
0
0
(0M.8)
0
C.5 Valid hypothesis testing (Theorem 4)
not trivial to have valid permutation test since we do not have IID examples. One naive permutation is permutate all the labels t (across batches i ). The resulting permutation testing is not uniform under the null (Figure 9).
In contrast, when we permute the label, we only ask for within-batch permeation (permute the index (4..... (r. each batch). See Table for an illustration. Let's prove this permutation is valid.
Proof. Under the null, p(0ly)= (0ly)) almost everywhere. According to the label mapping & F, label independent of given y. We first show that label independent á.
In general conditional independence does not lead to unconditional independence. But because here design the generating process by F(y, I, o) = (y)m,()) (Ø) y, 0). Under the null we have π (y, L, Φ) = ne ( (#lu)) Hence needs to be independent of Φ.
For any classifier because c(Ø) is a function of o, c(Ø) and t also independent. Let (o) be the marginal distribution of in π.
Now we are computing the cross entropy (with respect to the population T). LPD Σn (c((),), tn) where U is the log score U(P,z) log P(x). It is a random variable be- cause we have finite validation-set. Now we are conducting permutation of label I, the permuted label & is an independent draw from the same marginal distribution t, ग. (t). Because of the indepen-
d dence, Es.; (c(()).t) t) E, Eg (c((),,)) t). hence the permuted LPD, LPD, where LPD, the
computed LPD in the 6-th permutation. Therefore Pr(LPD < x)) = Pr(LPD, S æ). In expectation (with an infinitely amount of repeated random permutation), this p-vuuue will be uniform on [0,1].
In practice, the computed p-value from finite permutation can only take values on a finite set {0, 1/B,..., 1)/ B, 1}. More precisely, under the null hypothesis, this permutation test p-value is discreetty-umiform on this set. This is because for any 0 S m S B
['
( mm) pm (m) B-m dp 1/(B + 1), Vm.
Pr(p Ptpp=//)
Hence, the permutation p-value based on B random permutations is uniformly distributed on the set {0, 1/B...., (B 1))BB11).
22

---

C.6 Why does classifier help test? (Theorem 5)
Let's recap the binary label-generating process:
1. Sample marginal p(y),
2. Sample Bernoulliw)) (w), where w1 = 1/M, w2 (M 1)/M.
3. Sample (T)(0))), where To p)0 गι q(0ly).
4. Return (z, 2))
Theorem state that the optimal classifier has a sufficiency propriety that (a) let è be the probability of label in the optimal classifier as per (7), and let лı and x be the one-dimensional distribution of this ē()) when (0, y) sampled from p(0, y) or from (())q(0l)) respectively, then (i) Conditional
on the summary statistic è, the label is independent of all features = (0, y). (ii) There is no loss of information in divergence as the joint divergence is the same as the projected divergence, D1 (p(p,)) = D1 (= Tz)-
Proof. There are three random variables, 0, y. and t. The optimal classifier è is the probability in this true joint:
(0, = Pr(t = ())
To show sufficiency or conditional independence, all we need to show is that, conditional any given value of č, Pr(t = (0, y). c) does not depend on (0, y) (in other words, Pr(t 1|(0, v). is function of only). This becomes obvious
Pr(t = (6 v), c(0,0) v)) Pr(t = 1|(0, y)) = 2(0, v).
Now we prove that there is "no loss" in divergence.
(1(p.g) w (L(p(0, ())l + (1- w) KL(g(),)))
We express the first term in ē
KL(p(0, v) || (0,(0,10 v) + (1 (1-w),(0,1))
KL( glt 0) ((0,8))
KT](T([( = 0) (()) KT 1)
This steps uses the chain rule of KL. divergence: KLIpp(,, v) q(.,)))))))))))) kL(p(z) g(x)] + KLIP(1)
x) | 44y [)]]
Using conditional independence:
x (0, 1 = 0) =0, =)
Hence KT.(T(0,, v/ē, 0) (1. 0. Therefore,
(TLpp(, = KI(=(E]t = 0) x())
where (c)) (ā/t= 0) + (1 (1=
Similarly,
KT(g(0,)))1r(00.))) (I((<([(== (())
x(0)
This proves Qp) Di(=%, T2).
C.7 The maximum discriminative generator (Theorem 2)
We save the proof of Theorem 2 in the end for its length.
The generator & contains a few degrees of freedom: the number of classes K, the number examples
L, and how to design and the label-feature pairs. In the binary labeling: 4 = 1, Φ (0,, y) and the
remaining 4 = 0, Øk = (0x, y) for 2 S A S M + 1. The multi-class 4h assigns labels 1:K as
46 tR k. cb (Perm ......
K +1
(19)
23

---

Before the main proof that the multiclass permutation creates the largest divergence, let's first convince that the multi-class classification produced higher divergence than the binary one.
In the binary classification with M = (one draw from and one draw from g)
1 (1 7) KL p(O)y), p(0ly) <(01u)
((00l))) p(0 v) + g
KL + 2
2
In the multi-class classification with M
D4(2.0) 4) KL (p(0'v)a(001)) p(0'l))9(011)))) g(0 (v))))))))))))
Use the joint KI. > marginal KI., we have
D1(p.g) 2 KL (p(o'l)). p(0'))) g(0 lv)
Likewise,
Di(p.g) KL ((00)u)) p(0 ()) + g((Pl)))
Add these two lines we obtain
Da(p, q) 2 (p.g.))
To prove that this multi-class permutation produced the uniformly largest divergence (across M, K, p, q) optimal, we organize the proof into lemmas to 9, For notation brevity, denote M M + 1 in these lemmas to avoid using index M
Lemma For an arbitrary integer any given output spuce V, and any input space X that has at least L elements, if there are functions mapping X². to y
f1, (II,...,, .......,.TL............... -yE-
satisfying the following propriety:
for any probability distribution πt X. when TI..... .... In are L iid random variables with law n. f1(Z1........ has the same distribution as fI(
then there must exist permutation of denoted by o(1 : L), such that ...,,
fa(x1
Proof. For any L distinct values 01, aL, ag E X, let be mixture distribution of L delta functions:
8(mm Pmm
where the m-th mixture probability is
Pm C
c-1=
Σ
1
2
m=1
is chosen such that Σ Pm 1.
Now that T1, ....... delta functions, For any sequence of input indices (uI...., uz) E (1:1)²
......,2 LL are IID random variables from this n. f1 .(. . is also a mixture of
II (1C/24m Ej ((1s=m))
Pr(f((1,,...,, (i(au)).... auz)
c 1 Dinni (E/= 1(0, (my=m)).m-1)
(20)
2

---

in which the power index can be written as
I Σ Iuu .).....
Σ
Σ ==1
Iuu = L)
=
Σ /=1
1(u m)
m=1
as an L-decimal-integer.
Next, study the law of
Pr(fa(#1....., x1) = fa(a)....,..)))) =
Because (2( xr.) and (ii F1) have the same distribution, (2, to match the value at which fi fi(11,..., LL) has probability (3..)).......)) to attain. Comparing with (20), this probability is only attained when Σ/= 1(uy m) 1, Vm. That is, there exists Ø, a permutation L. such that Wiy...; ur
an.) needs
Matching the value of f2(x),...,, TT.) we obtain
fe(a = fi( ao())).
Because the choice of vector ( arbitrary, we have
f2(x11...... xL) fi(Fo()) ... (o()))
Because augmentation increases divergence, for the purpose of finding the largest divergence, to produce the largest divergence, we only need to consider an augmented generator that includes all
dang (y, ()...)..)) + t(.),,),t1).... ................................ ((ŌK, 9), tx))
It is enough to consider the generator of which oR ob (b,,....... 01) are K functions of (01,..., Lemma 11. For any augmented generator Haro satisfies IF, the null, there must exists K permutations (1 L),..., (L)), with the convention @ L) = 1: such that
(x.......,02) d1 (O.n(1)))))))
Proof. Use Lemma (10) for (K times.
Lemma 12. For any augmented generator garg satisfies IF, all feature-label generator can be replaced by a permutation, &z oz(O.....0.) = (0or(1).... (oa(L)), while the divergence does not decrease.
Proof. From the previous lemma.
2x ( 02) = &1 (1(o
The augmented feature is now (Φr (1.)) y), transformation of (0., (1:))) v). Using the raw feature (0o.(1:))) y) keeps the divergence non-decreasing.
Now that we only want to consider permutation- generators: there exists a K, and K permuta- tions (1 L).
s bang
(v, .......) + (((1), 4),, ((ox, ())) &n (0,,...... 0r,) = (o.(1:))
Given a o q, any permutations 0g (1:L) contains one copy from p and (L - 1) copies from g. It suffices to only consider those permutation (o1(1:)) whose distributions are unique.
Lemma 13. Civen any p + and L fixed, assuming all (1.1) has different distributions, then D(g, 4erg) is increasing function on K.
25

---

Proof. Using that joint KL is always not smaller than the KL divergence of sub- sub-coordinates.
Lemma 14, Among all permutations o01,....0p the maximum number of distinct distributions are K L.
Proof. The total number of permutations is L!. Because 0..1. are IID given y. the permutation of index 2: I. does not change the distribution. When p & q, the total number of distributionally distinct permutations is LL - 1)! = L M + 1.
clear that the proposed multi-clas permutation 4' attains this maximum number of distinet distributions, which proves its optimality among all generators from all previous lemmas, thereby the proof of our Theorem 2.
26


Reviews:
[{"review_id": "KjJWV6CZIg", "reviewer": "Reviewer_eYcd", "summary": "This work generalizes the well known Simulation-Based Calibration method to a setting where the tests on the posterior approximation consistency are based on classifiers. The paper has a nice balance between theoretical discussion (Section 3) and more practical issues (Section 4) making it easy for readers to appreciate different aspects of their proposal.\n\nThe paper is solid and tells a nice story from start to finish, but I would be tempted to point out the limitation of assessing only the quality of the posterior approximation **in average**. The authors mention this downside but don't give any insight regarding in which situations this might be a problem and when it could be OK to abandon local information.\n\nFurthermore, it comes somewhat as a surprise the fact that the authors use a **discriminative** approach for assessing the quality of the posterior approximation, but make absolutely no mention to the well known classifier two-samples test (C2ST).\n\nIn all, I'm rather satisfied by the paper as a contribution to the SBI literature, but I'm not certain of how wide is its scope for the NeurIPS conference. I'm giving it, therefore, a score of 6 (borderline accept), and would not be surprised if in the end it was not accepted to the conference main track.", "questions": "- Why is your Figure 2 **after** the set of Figures 3-6? This is very disturbing...\n\n- How useful it is to have a diagnostics to $q(\\theta \\mid y)$ which is only valid in average over $y$? Can you tell a bit more about this?\n\n- The class imbalance in your Example 1-4 will probably make your classifier have a lot of trouble to discriminate one class from the other. Indeed, it will become more of an \"anomaly detection\" kind of thing than an actual classification. You don't make much comments on this but I think it would be very important to explain whether this is a difficulty or not.\n\n- Can you add a reference to your claim in Line 82?\n\n- Line 188: What do you mean with \"classifier $c$ is good enough\" ?", "limitations": "I've mentioned the limitations in the previous fields.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The authors managed to recast the simulation based calibration (SBC) method in a discriminative framework. This is very nice as it helps better understanding such method and see how it compares to other approaches.\n- The three theorems presented in Section 3 are of great practical utility and ensure that even though we only estimate a classifier between samples from a fixed-size dataset, we can be sure that its performance will serve as a lower bound to other theoretical quantities. In other words, the fact of actually having to abide with the \"real life\" setting of training a classifier from limited data is not a problem.", "weaknesses": "- The mathematical notation can be sometimes a bit confusing. For instance, it is not always easy to discern what is a scalar input variable to a function with what is a random variable samples from a pdf.\n\n- There's no comment on assessing the quality of the classifier trained on calibration data and then used to build the statistical tests. \n\n- It could have been of interest to consider more examples on simulated data to illustrate the procedure. For instance, using examples from the `sbi-benchmark` which are becoming more and more used in the simulation-based inference literature.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are grateful for your careful reading and constructive comments. We have now addressed your concerns to improve and clarify the manuscript. Please find below a detailed point-by-point response to all of your comments and questions. \n\n**Figure 2 index**\n\nSorry for this confusion. We will reindex.\n\n**How useful it is to have a diagnostics to  which is only valid in average over?**\n\nThank you for suggesting this important aspect! Simulation-based diagnostics are averages. Our method is consistent with the rest of the simulation based calibration (SBC) literature (e.g., the seminary papers of Cook et al. 2006 and Talt et al. 2018) in computing a measure that is an average over y. (This is often called a “global” diagnostic in the SBC literature.) SBC was originally designed to validate if *computer software* accurately draws samples. For this application, the “in-average” assessment is a feature, not a bug, especially for amortized inference. This is our main motivation for using such a measure.\n\nThat said, the “local” diagnostic can be a useful goal as well: when the inference is only intended to run once, and the user is only interested in the inference given one  $y$. Our method is still relevant because KL$(p(\\theta|y), q(\\theta|y))=0$  if and only if $p(\\theta|y) =  q(\\theta|y)$ almost everywhere (Theorem 1). With modern inference routines such as MCMC or neural posteriors, we do expect that the “global divergence” could often achieve zero, which guarantees a zero local divergence with probability 1. \n\nLastly, our classifier could be adapted to “local” calibration (last paragraph on page 9): it is enough to look at the classification performance around a small neighborhood of the observed data $y_obs$ in the simulation table. It underscores the flexibility of our framework, and we leave it for future work.\n\n**The class imbalance in your Example 1-4 will probably make your classifier have a lot of trouble to discriminate one class from the other.**\n\nWe apologize for this confusion, but our novelly developed multiclass label-mapping and multiclass classifiers are always *balanced* (example 4, where we create M+1 labels from M+1 classes from each simulation draw). Hence increasing the number of inference draws (M) is not a difficulty for us. We have developed asymptotic theory and empirical validation to this big M limit. \n\nWe refer to our **Share Response #3** on the detailed discussion on the balanced multiclass classifier.\n\n\n**Can you add a reference to your claim in Line 82?**\n\nLine 82 states that, when the inference is exact ($p(\\theta|y) = q(\\theta|y)$ almost everywhere), and we train a binary classier to distinguish $(\\theta, y)$ and $(\\tilde \\theta, y)$, then the expected test log predictive density of a binary classifier in a 1:M imbalanced labels could be no higher than the negative binary entropy h(w) := w log w + (1 - w) log(1- w) of a Bernoulli distribution with w := 1/(M + 1). Intuitively, under the null, the best classifier (in terms of the expected log predictive density) is a Bernoulli(w) distribution. This statement is a direct consequence of our Theorem 1, Equation 8.  More detailed derivation can be found in the supplement, starting with line 589. We will clarify this point in the text.\n\n\n**Line 188: What do you mean with \"classifier is good enough\"?**  \n\nThis statement was unclear and we apologize. We will revise Line 188. What we meant was that (a) for any classifier, the divergence estimate from test data is always a lower bound of the divergence and (b) when the classifier is “good”, such classification performance is further a “good” estimate of the divergence. Here a classifier being “good” means it achieves a low test data error or equivalently a high expected log predictive density.\n \n\n**It comes as a surprise that the authors make absolutely no mention to the well known classifier two-samples test (C2ST).**\n\nIt was an oversight that we did not cite the phrase “classifier two-samples test” in the related literature, for which we apologize. We will add citations to the general C2ST framework in the revision. We stress that our paper is not simply applying C2ST to SBI. We generalize C2ST by incorporating a label mapping, which allows autocorrelation and multiclass classification for the two sample tests. We would like to refer to our **Shared Response #2** for a detailed discussion of our method in relation to the original C2ST.\n\n \n **There's no comment on assessing the quality of the classifier trained on calibration data and then used to build the statistical tests.** \n\nAt a high level, the techniques used to try to make the classifier generalize well are the same as any classification problem where the goal is to generalize to test data. Our Section 4.1 gives practical recommendations on network training and feature engineering. \n\nBesides, even when the learned classifier is not optimal (as would be common in practice) the divergence we estimate is always a valid lower bound and hypothesis testing is valid. Indeed, for any classifier, the estimated divergence is always a valid lower bound of the actual divergence (Theorem 1), and the proposed hypothesis test is always valid (Theorem 4). A good classifier only helps increase the power and reduce ​​false negatives (false positives are always controlled). We do observe in finite sample examples that the power from our test is higher than others. \n \n**It could have been of interest to consider more examples of simulated data to illustrate the procedure.**\n\nThank you for your suggestions. To showcase our method, we now add three simulation data examples from this `sbi-benchmark` repo, the Gaussian linear, the Gaussian mixture, and the simple likelihood complex posterior (SLCP) example. We run our calibration method on these three datasets with varying inferences settings and it picks up positive divergence estimates when the inference is not exact. Please see our **Shared response #1** and attached pdf file."}, {"review_id": "ghvslpppKN", "reviewer": "Reviewer_TMBR", "summary": "In this paper, the authors consider conducting frequentists tests for simulation-based inference. Specifically, the authors want to test if an inference engine q(theta | y) produces samples from the true posterior p(theta | y). The key idea is to use classifiers to determine similarities between samples generated from the true and approximated posterior (similar to GAN). \n\nThe foundation of the algorithm is laid in section 3, where the authors show the negative cross-entropy is associated with a generalised divergence between p and q. In particular, D4 becomes KL divergence between p and q as the number of posterior draws goes to infinity. Empirical study validate the theoretical results and shows the effectiveness of the algorithm. \n\n", "questions": "If one wants a divergence KL(p(theta|y, q(theta|y)), one can directly approximate it without using a classifier, right? \n\nKL(p(theta|y), q(theta|y) = KL(p(theta, y), q(theta|y)p(y)). \n\nSuppose you can sample from p(theta, y), q(theta|y) and p(y). Then we can directly approximate above KL using f-GAN,  and perform a permutation test or compute CI using bootstrap.  \n\nThis approach sounds like a simpler way to obtain a divergence measure between p and q, so I wonder why the authors did not explore this idea nor use it as a benchmark method. \n\n----------\n\nRevision after author's response. \n\nThis approximation cannot be easily applied due to the IID-ness and the sequential nature of SBI as authors mentioned. \n", "limitations": "Yes", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This paper is generally well-written and the idea is intuitive and sound.\n2. This paper addresses an understudied problem (SBC) in the literature, and the expansion from a single-dimensional rank statistic to a multidimensional test statistic seems to be a significant improvement\n3. The algorithm provides theoretical support (Section 3) to the proposed algorithm. ", "weaknesses": "1. From the methodological point of view, using a classifier to compare distributions is not a new idea (as the authors have discussed in Section 6. \n2. When training a classifier, it seems authors are facing very imbalanced classes and training an imbalanced classifier itself can be difficult. I am wondering if this would cause issues in the tests performed later.  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\nWe are grateful for your careful reading and insightful comments. We have now addressed your concerns to improve and clarify the manuscript. Please find below a detailed point-by-point response to all of your comments and questions. We wish our explanations to help bring clarification. \n\n**From the methodological point of view, using a classifier to compare distributions is not a new idea (as the authors have discussed in Section 6.** \n\nUsing the classifier two sample test (C2ST) to perform two sample test is not new. However, the traditional classifier two sample test does not directly apply to the SBI problem because of the sequential sampling and the autocorrelation. Our results not only formulate the SBC into a classifier-two-sample-test task and interprets the divergence in posterior distribution, but also develop the result to incorporate autocorrelation and extra log density information. We also develop a more general label mapping framework. This framework and the novel multiple-class classifiers approach may be applicable to other two sample test problems. \n\nPlease see our **Shared Response #2**  for details on why our paper develops the classifier two sample test approach.  \n\n\n**When training a classifier, it seems authors are facing very imbalanced classes and training an imbalanced classifier itself can be difficult. I am wondering if this would cause issues in the tests performed later.**\n\nWe apologize for this confusion, but our novelly developed multiclass label-mapping and corresponding multiclass classifiers are always *balanced* (Example 4 in Section 2, where we create M+1 labels from M+1 classes from each simulation draw; one label from each class). Hence increasing the number of inference draws (M) is not a difficulty for us. We have developed asymptotic theory, convergence rate ($M\\to \\infty$) and empirical validation to this big M limit. \n\nIn comparison, the naive one-to-M binary classification could run into the imbalanced classification challenges you have raised, underscoring the benefit of our generalization of label mapping. That said, our method still supports binary classifiers with reweighting (Theorem 8 in Appendix). In our experiment, we find that even under highly imbalanced binary labeling (M=1000), with appropriate reweighting, the binary classifier is still able to detect the inference flaws and output accurate divergence estimates.  \n\nWe would like to refer to our **Share Response #3** for a more detailed discussion on why our multiclass method leads to balanced classification. \n \n\n**[f-GAN]  If one wants a divergence KL(p(theta|y, q(theta|y)), one can directly approximate it without using a classifier, right? \nSuppose you can sample from p(theta, y), q(theta|y) and p(y). Then we can directly approximate above KL using f-GAN, and perform a permutation test or compute CI using bootstrap.**\n\nThank you for your insightful remark. We note that the f-GAN cannot be directly applied to obtain a useful divergence for similar reasons as the C2ST discussed in the shared response.\n\nBackground: Nguyen et al. (2009) developed an M-estimation to compute the f-divergence from two samples. This estimate was used in f-GAN to replace the binary classifier. Belghazi et al (2018) developed another tighter sample-based M-estimation for KL-divergence.\n\nThe sample-based M-estimation in Nguyen et al. (2009) or Belghazi et al (2018) is not directly applicable to the simulation based calibration (SBC) setting because:\n(1) the method in Nguyen et al. (2009) and f-GAN applies to IID samples, but the SBI joint draws $(\\theta, y)$ and $(\\tilde \\theta_m, y)$ are not IID draws from two simulators since all the $\\theta$ and $\\tilde \\theta_m$ are paired with an identical $y$, and\n(2) when the inference $q(\\theta|y)$ is MCMC, there is additional auto-correlation in $\\tilde \\theta_m$ draws.\n\nSimilar to how we can see our method as a generalization of the C2ST, it is likely that we could adapt Nguyen et al. (2009) or Belghazi et al (2018)'s KL estimate to the SBI setting. Indeed, it is the contribution of our paper to formulate the traditionally rank- and histogram-based SBC problem into a sample-based joint-space discriminative task. Our framework is flexible to incorporate the inclusion of other sample-based divergence estimates and loss functions, such as the aforementioned sample KL estimates, Wasserstein distances, and the integral probability metrics. However, we leave this as future work."}, {"review_id": "h1GTM9aicR", "reviewer": "Reviewer_gi2H", "summary": "In this paper the authors focus on the challenge of comparing two conditional distributions p, q from their samples. In particular, this is useful as a check for bayesian computations.\nTo achieve this, the authors propose the use of a probabilistic classification approach where they create a new dataset combining samples from p and q and labels related to the particular distribution. Here, the authors propose different approaches and make connections to previous work.\nGiven the combined dataset, a classifier tries to predict the label from the features. Failing to do so, indicates that the classifier could not discriminate between the two distributions.\nThe performance of the probabilistic classifier is then used to estimate the divergence between the distributions and for hypothesis testing. \nIn the empirical section, the authors compare against SBC and show improvements in power given a reduced number of samples.", "questions": "I'm thinking of your paper as an approach to test whether the labels and the generated features are connected. This reminds me of [1] where they evaluate the mutual information via neural networks. I'm wondering both if it makes sense in your setting to estimate the divergence via the Donsker-Varadhan representation, and in their setting if your feature generation approach is beneficial. Do you have any opinions on that?\n\nThe same has applications for other approaches on independency testing via classifiers where your feature generation approach also has potential.\n\n[1] Belghazi, Mohamed Ishmael, et al. \"Mutual information neural estimation.\" International conference on machine learning. PMLR, 2018.", "limitations": "There are no potential negative societal impacts from this work. The authors mention limitations in a future work section. ", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper analyses different approaches for the generation of samples to be used by discriminators. They then present theoretical developments and demonstrate their results empirically. The sampling approaches and theoretical treatment are novel, of high quality and well written.\nThe main advantage here is that even for a sub-optimal classifier, one gets improvements in data efficiency. ", "weaknesses": "The main weakness of the paper is that in the empirical evaluation, the authors do not explore the case where the samples present more challenging behaviors such as auto-correlations, different types of imbalances, etc.\n\nOther minor suggestions:\n- The colors of the thetas in the introduction are hard to see. \n- Describing why the bounds are tight could further help the reader in section 2.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\nWe are grateful for your careful reading of our manuscript and your constructive review comments. We have now addressed your concerns to improve and clarify the manuscript. Please find below a detailed point-by-point response to all of your comments and questions. \n \n**The main weakness of the paper is that in the empirical evaluation.**\n\nThank you for your suggestions on the experiments. To showcase our method, we now add three simulation data examples from this sbi-benchmark repo, the Gaussian linear, the Gaussian mixture, and the simple likelihood complex posterior (SLCP) example. We run our calibration method on these three datasets with varying inferences settings and it picks up positive divergence estimates when the inference is not exact. Please see our Shared response #1 and attached pdf file.\n\n**Describing why the bounds are tight could further help the reader in section 2.**\n\nThank you for your suggestion. We use Section 2 to present four examples to help the reader quickly see the heuristics and intuition: why we can formulate the simulation based calibration into a joint space classification, and various label mapping could lead to different divergence estimates. We state the rigorous theory in Section 3, but we felt this was too abstract without having some examples first. We provide detailed derivation of all theorems including the tightness of the bound in Supplement C.    \n \n**I'm thinking of your paper as an approach to test whether the labels and the generated features are connected. This reminds me of [1] where they evaluate the mutual information via neural networks. I'm wondering both if it makes sense in your setting to estimate the divergence via the Donsker-Varadhan representation, and in their setting if your feature generation approach is beneficial.**\n\n\nThank you for your insightful remark. It is an interesting connection, which we did not anticipate! Nguyen et al. (2009) developed an M-estimation to compute the f-divergence from two samples using conjugate dual function theory. Belghazi et al. (2018) developed another tighter sample-based estimate of the two-sample KL-divergence using the Donsker-Varadhan representation. Both estimates have been applied in GAN to replace the traditional binary classifier. \n\nFirst, let us clarify what makes the simulation based calibration (SBC) setting different from the generic two-sample divergence computation—The SBI joint draws $(\\theta, y)$ and $(\\tilde \\theta_m, y)$ are not IID draws from two simulators since all the $\\theta$ and $\\tilde \\theta_m$ are paired with an identical $y$. Furthermore, when the inference $q(\\theta|y)$ is MCMC, there is additional auto-correlation in $\\tilde \\theta_m$. Both the across-class and within-class dependence violates the IID sampling assumption in Nguyen et al. (2009) and Belghazi et al. (2018). In comparison, our proposed multiclass classification (example 4) solves both issues. Even better, our multiclass classification is always balanced. Indeed, the main motivation of our multiclass classification development is to address imbalanced and auto-correlated samples, not to purely develop a sample based KL divergence estimate, and any finite M produces a divergence metric as well. \n\nSecond, similar to how we can see our method as a generalization of the C2ST, it is likely that we could adapt Nguyen et al. (2009) or Belghazi et al (2018)'s KL estimate to the SBI setting. Indeed, it is the contribution of our paper to formulate the traditionally rank- and histogram-based SBC problem into a sample-based joint-space discriminative task. Our framework is flexible to incorporate the inclusion of other sample-based divergence estimates and loss functions, such as the aforementioned sample KL estimates, Wasserstein distances, and the integral probability metrics. When some part of density information p or q is known such as in SBI and GAN, it is straightforward to include them as features in the learning of the Donsker-Varadhan representation as we have done in Section 4. We leave this extension for  future work. \n \nLastly, the “byproduct” of our calibration method is a novel sample based KL divergence estimate using multi-class classification (example 4), which can be used in other applications. A direct consequence of our Theorem 3 is that, in a generic two-simulator setting, the divergence estimate obtained from any multi-class classifier and any simulation size M is always a lower bound to KL divergence from p to q—and this estimate becomes tight when the classifier is optimal and M goes to infinity. In the paper we have proven its convergence, yielded the convergence rate (theorem 3, where the constant is further the chi-squared divergence), and shown the empirical evidence. Perhaps the most interesting direction for future investigation (although unrelated to this paper) is to compare this new multi-class-classifier estimate to its orthogonal counterpart, the Donsker-Varadhan representation and the f-dual, in the IID sampling case (which is NOT SBC, but includes the GAN, two-sample test, and independence test) such that all three methods are applicable."}, {"review_id": "r03GM4wvl5", "reviewer": "Reviewer_aBBC", "summary": "This paper presents a classifier based approach to provide a measure of miscalibration in Bayesian computation, including for methods such as Approximate Bayesian Computation (ABC) and Simulation-Based Inference (SBI) methods like neural posterior estimation. The method enables the test statistic to be learned from data and  provides an interpretable divergence measure. The method is a form of two sample testing applied in the amortized and simulation based inference settings. Beyond standard approaches to using classifiers for two sample testing, the authors develop several classifier based approaches which use a form a \"label mapping\", and provide theoretical work on the validity and efficacy of these approaches. Experiments are provided to test the method for posterior inference in cases where the posterior is known, and provide empirical verification of the theoretical results. The paper also compares to the commonly used simulation based calibration method, and show clear improvements with the methods proposed in the paper. Experiments on cosmological data are also explored.", "questions": "Could you provide more details on the posterior model in the multi-variate Gaussian example?\n\nCould you provide experiments on more complex but controlled examples, for instance the SLCP problem from Papamakarios et. al? This could provide more evidence on the utility of the method in SBI settings.\n\nHow do these methods scale with feature and parameter dimension? Does this affect the ability to attain a useful diagnostic?\n\nYou mention that the presence of nuisance parameters do not affect the quality of the diagnostics. However, the Neyman-Pearson lemma does not guarantee that the likelihood ratio is the uniform most powerful test for a given size in the presence of nuisance parameters. Perhaps I am not fully understanding the text, could you explain if this has any ramifications for the proposed tests?\n\nCould you add references to other simulation based inference work using classifiers for diagnostics? In addition, recent work on using neural networks for two sample testing (e.g. Grosso et. al, https://arxiv.org/abs/2305.14137) may be relevant to discuss in the related work. ", "limitations": "The authors provide clear discussions of limitations. It would be interesting to also know if there are limitations dependent on the dimension of the features or parameters, and how this may impact the number of samples needed to attain a useful test.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "This is a good paper, it provides:\n-  a clear description of the challenge they hope to address, specifically providing a better and more statistically interpretable measure of miscalibration over simulation-based calibration\n- a clear description of the proposed methods, including how to use label mapping to develop calibration based diagnostics, and how this differs from other approaches\n- theoretical grounding for the calibration based measures, including the expected behaviors in large sample limits\n- useful discussions on implementation details\n- useful discussion on the legitimacy and power of the tests.", "weaknesses": "While the idea of using classifiers for two sample testing is not new, the authors do develop new methods based on label mapping. While the new classification approaches are discussed and will likely be quite useful for simulation based inference, some attempts at using classifiers for testing posterior inference quality have been performed before (albeit not directly with some of the new classification methods discussed in the work). For instance, Vandegar et. al \"Neural Empirical Bayes: Source Distribution Estimation and its Applications to Simulation-Based Inference\" AISTATS 2021, use the AUC of a classifier as a diagnostic in a simulation based setting.  \n\nFurther development of the experiments would be quite useful. A very simple experiment is shown, allowing an exploration of the methods, but some experiments more clearly showing the quality of the method in a controlled but more complex simulation based inference setting could provide useful insights to the readers. For instance, a common example is the Simple Likelihood, Complex Posterior (SLCP) problem from Papamakarios et. al, \"Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows\", AISTATS 2019. As it stands now, it is difficult to asses the quality of the method in the more complex cosmological data set experiment.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are grateful for your careful reading and constructive reviews. Please find below a detailed point-by-point response to all of your comments.\n\n**Further development of the experiments would be quite useful/more complex but controlled examples, for instance the SLCP problem?**\n\nThank you for your suggestions on the experiments! We find the SLCP example interesting. To showcase our method on a wider range of problems, we now add three additional simulation data examples from the `sbi-benchmark` repo: the Gaussian linear, Gaussian mixture, and the SLCP example.  We run our calibration on these three datasets with varying inferences settings and we obtain positive divergence estimates when the inference is not exact. Please see our **Shared Response #1** and attached pdf file.\n\n**Could you provide more details on the posterior model in the multivariate Gaussian example?**\n\nWe would like to refer to our Supplement B.2 for experiment details---It was an oversight that we did not state more clearly that this information was there and we apologize. We will state this clearly in the revision.\n\nIn the closed-form Gaussian experiment in Section 5, we consider an easy example: a normal likelihood  $y\\sim$MVN($\\theta$, Id) and a normal prior $\\theta \\sim$MVN(0, Id). The true posterior is a closed-form Gaussian.  We consider a sequence of corrupted inference by adding a bias shift to the posterior mean or multiplying the posterior covariance matrix by a bias factor. Notably, we have kept this true posterior distribution and the sampling corruption mean-field to make the calibration task easier for the traditional SBC rank-test. \n\n**How do these methods scale with feature and parameter dimension? Does this affect the ability to attain a useful diagnostic?**\n\nFirst, traditional rank-based SBC faces challenges in high dimensions due to (1) the difficulty of checking interactions between dimensions and (2) the expense of multiple testing. Running SBC in high-dimensional problems is one of our initial motivations. In the real data cosmology example, the input dimension ($\\theta \\times y$ joint space) of the classifier input we trained is ~1000.\n\nStrictly speaking, high dimensions in parameters do not pose any conceptual difficulty for our framework. However, to get useful calibration it is necessary that the classifier is able to distinguish the two or multiple classes—if the classification problem is too challenging, the learned classifier could be far from the optimal one, which would mean our divergence estimate is a weaker bound and more false negatives would occur in testing. Our Section 4 offered practical recommendations to help dimension scaling, including using (1) pre-learned features for dimension reduction, (2) statistical features such as log p and log q whenever available, and (3) network symmetry.\n\n**the presence of nuisance parameters do not affect the quality of the diagnostics. However, the NP lemma does not guarantee that the likelihood ratio is the uniform most powerful test for a given size in the presence of nuisance parameters.**\n\nThank you for raising this interesting question! In the paper, we discussed the nuisance parameter when the full parameter space can be partitioned into two parts, say $\\theta$, the parameters of interest, and $\\phi$, the nuisance parameters.  We are only interested in the marginal sampling quality of the $\\theta$, that is if $p(\\theta|y)= q(\\theta|y)$.  We argue that we only need to restrict the classifier to use the target dimensions $\\theta$ and $y$, while other nuisance parameters can be discarded from the classifier. As a consequence of our Theorem 1, the estimated divergence from this classifier is the conditional divergence between $p(\\theta|y)$ and $q(\\theta|y)$. The nuisance parameter has no ramifications.\n\nThe Neyman-Pearson lemma does not apply to nuisance parameter problems as the null and alternative are composite. Here we do not have a composite hypothesis since we just look at marginal posterior distribution $\\theta|y$. To be clear, when looking at the $\\theta$ margin, the density ratio to be learned from the classifier is $p(\\theta|y)/q(\\theta|y)$, not the original full density ratio $p(\\theta, \\phi|y)/q(\\theta, \\phi|y)$, which is what the likelihood ratio means in the NP lemma. \n\n**Could you add references to other simulation based inference work using classifiers for diagnostics?  Recent work on using neural networks for two sample testing may be relevant to discuss in the related work.**\n\nThank you for your suggestions on references, which we will add.\n\nFirst, using the classifier two-sample test (C2ST) is not a new idea. However, the traditional classifier two sample test does not directly apply to the SBI problem because of the sequential sampling and the autocorrelation. Our present not only formulates the SBC into a classifier-two-sample-test task and interprets the divergence in posterior distribution, but also develops the traditional C2ST to incorporate autocorrelation and extra log density information, as well as a general label mapping framework. Please see our **Shared Response #2**  for details of our in relation to C2ST. \n\nIt was an oversight that we did not cite Vandegar et. al (2021) who used ROC AUC to compare (1) samples from $p(\\theta|y)$ and $q(\\theta|y)$, and (2) the data y and the simulated posterior predictive $\\int p(y|\\theta) q(\\theta|y) d\\theta$. We apologize for the oversight we will add in the revision. However,  Vandegar et. al (2021) addressed a different task than ours. Assessment (1) requires samples from the exact posterior $p(\\theta|y)$, which is different from our intended application in diagnostics. Assessment (2) is related to the posterior predictive check (e.g., Gelman 1996),  a way to examine the model specification $p(y|\\theta)$ rather than (purely) computation accuracy. Our diagnostics are only intended to diagnose the computation, not to check the model correctness (line 377)."}]

Author rebuttal:
We would like to thank all the reviewers for their careful review and insightful comments. In addition to the point-by-point response we make to each reviewer individually, below we will give three shared responses, including additional experiments (in the submitted .pdf file). 

### 1. Additional experiments

As suggested by reviewers, to showcase our method on a wider range of examples, we now add three simulation data examples from the `sbi-benchmark` repo:
1. the simple likelihood complex posterior (SLCP) example,
2.  the Gaussian linear model,
3. the Gaussian mixture model.

The attached pdf file summarize the experiments. Here we run our calibration method on these three datasets. In each dataset, we sample prior draws from the default prior, and run adaptive  No-U-Turn Sampler (NUTS). We want to check the quality of the sampler after a fixed number of iterations. To this end, we run the sampler with various numbers of iterations from 2 to 2000 (for each point, we use the same number of iterations for warm-up and for sampling, and the warm-up samples were thrown away. So 1000 in the $x$-axis means 2000 total MCMC iterations were run and the last 1000 were kept as the inference output), and for each number of iterations, we run our classifier calibration. The y-axis is the estimated divergence at the given MCMC iterations, and we visualize the plus and minus one standard error from our method. A positive divergence indicates a mismatch between the true posterior $p(\theta|y)$ and the inference $q(\theta|y)$, while the divergence being zero means that the posterior inference is exact, i.e., $q(\theta|y)= p(\theta|y)$ almost everywhere. 

From the figures attached in the pdf file, in all three examples, we are able to detect the inference flaws and return a positive divergence estimate after a few iterations. The estimated divergence also does converge to zero for more iterations where we expect inference to be near exact.

### 2. [C2ST] Our paper is not simply applying the classifier two-sample test (C2ST) to SBI calibration. We generalize the classifier two-sample test.

Indeed, using a classifier to perform a two-sample test is not a new idea. We did not intend to imply this but we acknowledge that the paper as submitted does not make this as clear as it should. That said, the classifier two-sample test is not directly applicable to simulation based calibration (SBC) due to four barriers:
1. In the past of SBC, it was not clear what space to run the classifier on, i.e. should it include parameters $\theta$, data $y$, or the joint, or include likelihood values, nor to interpret results as a divergence between the true and inferred posteriors.
2. More importantly, the classifier two-sample test only works with IID examples (there are two distributions P and Q, and we have IID observations from P and Q respectively.)  The SBI joint draws are not IID because the simulation table contains shared y. Think about our binary classification scheme, the example from class 0 is $(\theta, y)$, and the examples from class 1 is $(\tilde \theta_1,  y), … (\tilde \theta_M,  y)$; they all share a same y component.  
3. When the inference is done MCMC, there is additional auto-correlation in $\tilde y$, which further violates the C2ST requirement.
4. For two-sample tests, C2ST only creates binary labels and may perform poorly with imbalanced classification. Using naive binary classification in SBI will face a 1:S imbalancement.

Our paper solves these issues by developing a general label mapping framework:
- For (1), our paper formulates the SBC problem into a sample-based joint-space discriminative task for the first time. We develop the relevant theory for the SBI calibration and prove its relevance to posterior divergence. We also make use of additional information like log-likelihoods in a unified framework.
- For (2), our theorem 1 and 9 extend the traditional IID two sample tests, and allows a sequential sampling. 
- For (3) and (4), we extend the straightforward binary classification to a general framework “the label mapping” (line 132). The traditional binary C2ST is now a special case of our framework (example). In contrast, the multiclass classifier (example 4 in Section 2) is always balanced and allows the examination of auto-correlated samples (equation 12).  

### 3. [balanced classifier] The naive binary classifier can suffer from imbalanced labels, but our novel multiclass classifier always has perfect label balances.
Our novelly-developed multiclass classifier (example 4 in Section 2) is always perfectly *balanced*: each simulation run creates (M+1) examples from (M+1) classes; one from each. See Page 3 for an illustration table. Our theorem 3 proves that the multiclass classification divergence converges to the meaningful KL divergence between p and q and derives the convergence rate, under the limit $M \to \infty$, the simulation in which the naive binary classification is infinitely imbalanced.

That being said, the binary classifier might be more intuitive to the users than the multiclass one, and we still support the use of binary classifiers in our calibration. In our experiment, we find that even under highly imbalanced binary labeling (M=1000), with appropriate sample reweighting (Theorem 8 in Appendix), the binary classifier is still able to detect the inference flaws and output accurate divergence estimates.  


Meta-review: The reviewers unanimously recommend acceptance (7-6-7-7). Several minor issues have been raised and discussed during the author-reviewer discussion period. The authors are encouraged to take this feedback into account in the final version of the paper.
Decision: Accept (poster)