An end-to-end framework for uncertainty-aware validation of generative models

Introduction

Many scientific analyses compare observed data with a reference model in search of significant deviations, which may signal new phenomena. The reliability of the reference model is therefore critical: any mismodeling can be mistaken for a genuine signal and lead to a false discovery. Such reference models are often built from high-fidelity simulations that require substantial computing resources. In particle physics, for example, the background from known processes is modeled by simulating the full chain of hard process, parton showering, and detector response. The demands on the reference model are particularly high in model-agnostic searches, for example anomaly detection [1], where any deviation from it across the full feature space may be reported as a signal.

Generative models have emerged as a way to significantly reduce the cost of such expensive simulations, acting as fast surrogates (or even replacements) [2, 3, 4, 5, 6, 7] or as learned likelihoods in simulation-based inference [8, 9]. Like any fitted model, however, generative models carry an inherent error due to the finite size of the training set. Generated events are moreover not independent draws from the true density: they all come from a single model (or ensemble of models) fitted to a finite sample, whose fluctuations are imprinted in the model. No amount of generated data can therefore add information beyond the training sample [10]. In practice, any trained model will eventually fail a statistical test if it is checked against enough data. This happens even when the model is, in principle, capable of describing the data perfectly (Sec. II). Whether this failure matters depends on the use case in the analysis workflow [11]. A few representative examples are listed in Figure 1. Errors can be tolerable when the result does not impact the downstream task; this includes examples when the model’s output is further validated by exact downstream checks, or when only qualitative fidelity is required, as is the case in exploratory studies or pretraining. But when the fitted

density serves as the reference model of a statistical test, as a background estimate, a null hypothesis, or a calibration input, its error is indistinguishable from a signal. To be used reliably in this role, a generative model must therefore come with uncertainties, and these must enter the statistical test itself.

For models with few parameters, estimating uncertainties follows standard procedures. However, for modern density estimators, particularly machine learning driven generative models, with high-dimensional inputs and large numbers of non-identifiable parameters, the problem becomes technically challenging, and attributing robust uncertainties to them remains an active area of research.

To address this challenge, we propose SIGMA (Score-Induced Gaussian Mixture Approximation), a generative model whose uncertainties are tractable and can be propagated into a statistical test. To validate it, we further introduce an uncertainty-aware goodness-of-fit test that incorporates the derived uncertainties.

SIGMA models the density with a Gaussian mixture model and describes its uncertainties in the tangent space spanned by its score functions around the best fit of the density. Unlike for neural networks, this simple Ansatz makes the asymptotic maximum-likelihood covariance of the parameters tractable. The test then propagates these uncertainties by promoting the conventional point-null Neyman–Pearson likelihood-ratio test [12, 13, 14] (also referred to as simple hypothesis) to a composite one [15], in which the nuisance parameters encoding the uncertainties are profiled, as is common practice in LHC searches [16]. Unlike in searches, the nuisance parameters are not systematic uncertainties of the experiment, but the statistical uncertainties of the fitted generative model itself. Because they enter the test in the same standard form, conventional systematic uncertainties can be profiled alongside them, so that a single composite test accounts for both the error of the fit and the systematics of the analysis. The test outcome is therefore determined by the compatibility of the data with the model within its uncertainties, rather than driven by the fit error.

Calibrated with pseudo-experiments drawn from the composite-null, the test retains its nominal false-positive rate even when the number of test events exceeds the training sample size, provided that the uncertainties faithfully cover the error of the fit.

The paper is organized as follows. In Sec. II, we motivate the problem with a controlled one-dimensional example. Section III reviews related work. In Sec. IV, we describe the components of our algorithm, from the density model to the uncertainty quantification and the design of the goodness-of-fit test. § 2.2 applies the method to a controlled two-dimensional toy study, where the true density is known, and to the QCD background in a learned embedding of particle jets, a setting relevant for anomaly detection searches. We conclude in § 2.3 and outline future research directions.

In App. A we further discuss the intrinsic limits of generative models and standard goodness-of-fit tests. Appendix B explores a related approach to quantifying and propagating model uncertainties with ensembles of neural-network density estimators, and finds that the resulting uncertainties are insufficient for the composite-null goodness-of-fit test to close on large observed datasets. Appendix § 2.1.2 contains technical details of the experiments and studies the effects of varying a hyperparameter, Appendix § 2.1.3 derives the asymptotic distribution of the composite-null test statistic under the null hypothesis, and Appendix § 2.9 studies the choice of $\sigma$, another hyperparameter in the goodness-of-fit test.

Powerful tests reject every generative model

A generative model with fitted parameters $\hat{\theta}$ only ever provides an approximation $p_{\hat{\theta}}(x)$ of the true data distribution $p_0(x)$. Following Refs. [17, 18], the discrepancy has two classes of sources, both epistemic, i.e. uncertainties about the model rather than about the intrinsic randomness of the data:

i) The model bias, a systematic uncertainty: the model may be unable to represent the true distribution exactly, however well it is trained. This bias depends on the choice of architecture and does not shrink with the size of the training sample. It is distinct from the systematic uncertainties usually considered in physics analyses, such as imperfect knowledge of the detector or of the theory, which concern how well the training data, typically a simulation, describe reality.

ii) The statistical uncertainty of the fit: the finite training sample limits how precisely the parameters can be determined (data variability), and training may converge to different solutions depending on its initialization (optimization variability). Both vanish in the limit of infinite data and ideal optimization.

Our work addresses the statistical uncertainty ii). Even in the best-case scenario in which the model family contains the ground truth, the finite size of the training dataset places a fundamental limit on the accuracy of the fitted density. This limit has a direct and often underappreciated consequence for the validation of generative models via goodness-of-fit tests.

To understand the effect, we consider the following one-dimensional experiment: samples are drawn from a standard normal distribution, $p_0 = \mathcal{N}(0, 1)$, and the generative model is chosen to be perfectly well specified, $p_\theta= \mathcal{N}(\mu, \sigma)$ with trainable parameters $\theta= (\mu, \sigma)$. The parameters are estimated by maximum likelihood optimization on $N_{\mathrm{train}}$ training points, and the fitted model is then compared with independent test samples of increasing size $N_{\mathrm{test}}$ by means of a goodness-of-fit test.

The fitted model, indistinguishable from the truth at small test sample sizes, is eventually rejected with high significance once $N_{\mathrm{test}}$ grows past $N_{\mathrm{train}}$. This degradation is not specific to the test choice: any consistent goodness-of-fit test eventually detects the finite-statistics error between $p_{\hat{\theta}}(x)$ and the true density $p_0(x)$, so the $p$-values decrease with $N_{\mathrm{test}}$ regardless of the test used. The gray markers in Figure 2 show this trend for the Gaussian example and the point-null likelihood-ratio test that will be defined in Sec. IV (see Eq. (17)) calibrated with pseudo-experiments drawn from the fitted model.

This behavior is not a pathology of the example but a generic and quantifiable consequence of information limits [10] (see App. A for a derivation). How fast the $p$-value curve degrades depends both on the sensitivity of the goodness-of-fit test employed and on the quality of the model design, i.e. how well the model design matches the properties of the true distribution. Crucially, however, the degradation itself cannot be avoided: no amount of architectural or training improvement makes a density estimated from $N_{\mathrm{train}}$ samples exact.

How, then, can generative models be used reliably for science? It is worth being precise about what this “intrinsic” failure means. Standard goodness-of-fit methods test a point-null hypothesis: the fitted model is treated as an exact description of the data. Since the fitted model differs from the truth with probability one, rejection at large $N_{\mathrm{test}}$ is guaranteed by construction. However, this is not the question of scientific interest. The relevant question is whether the data are compatible with the generative model within its statistical uncertainty, the spread of fits that would result from different training samples and different runs of the optimization. This is a composite-null hypothesis. A test that ignores this uncertainty conflates genuine mismodeling with the unavoidable statistical fluctuations of the model fit, and its outcome is driven by the relation between $N_{\mathrm{test}}$ and $N_{\mathrm{train}}$ rather than by the quality of the model. On the other hand, the bias introduced by an incorrect model is one that uncertainties cannot solve, and needs to be addressed at the level of model design, e.g. by replacing parametric background functions that fail as the dataset grows with more flexible Gaussian processes [19]. The statistical uncertainty of the fit, however, can be propagated into the goodness-of-fit test itself, promoting the point-null to a composite one. With the uncertainties included, the median $p$-value remains close to its nominal value of 0.5 over the full range of test sizes we probe (blue markers in Figure 2).

Related work

a. Generative models and uncertainty quantification. Uncertainties on generative models have been estimated with Bayesian approaches [20, 21] and Monte Carlo dropout [22, 23]. While both return a spread of predictions, they do not come in a form that can be written down as a parametrized density function. This makes them incompatible with composite likelihood-ratio tests, which require a parametrized family of null densities.

The frequentist ensembling in Ref. [24] provides one such parametrization, describing the uncertainty through the weights of a fixed ensemble of trained models. Ensembling can be performed for any architecture, but the resultant uncertainties span only the coefficients of the trained models. As a result, the uncertainties cover the true density only under the strong assumption that it lies within the span of the trained ensemble members. Ensembling is also used in experimental searches, where the spread of independently trained models provides the uncertainty of a generative background model in the final fit [25]. Our SIGMA approach follows a similar logic but spans the full parameter space of the generative model, weakening the well-specification assumption, at the price of requiring a tractable model likelihood.

b. Data amplification and the statistical value of generated samples. A closely related line of work asks whether samples drawn from a generative model carry statistical precision beyond that of the training sample, quantified through an amplification factor [26, 27, 28, 29]. Refs. [28, 29] quantify the uncertainty of the learned densities, but do not propagate any uncertainties to the point-null hypothesis. Ref. [10] establishes that the generated data cannot contain more information than the data on which the model was trained, and Ref. [30] derives an information-theoretic bound on the achievable gain, posing a theoretical limit to amplification.

We address a prerequisite for trustworthy amplification: whether the fitted density, together with its estimated uncertainty, is compatible with independent data from the true distribution.

c. Validating generative models. Previous work typically validates the learned density through a low-dimensional summary test statistic: quantile occupancies [26], projected Jensen–Shannon divergences of high-level observables [27], the per-bin coverage of Bayesian uncertainties [28], the Kolmogorov asymptotics of a one-dimensional statistic [29], or classifier-based scores [31, 32], including the full distribution of classifier weights [33]. Since all of these tests treat the fitted density as exact, they eventually reject even a model that is as good as its training statistics allow, once $N_{\mathrm{test}}$ is large enough to resolve the error of the fit (Sec. II).

Two-sample tests for validating generative models have been compared in Refs. [34, 35]. The latter finds the Neyman–Pearson test to be overall the most robust, outperforming marginal-based comparisons, at a higher computational cost. The Neyman–Pearson test with a learned alternative hypothesis [12, 13, 14, 36] compares the hypotheses in the full feature space, assumes little about the position or shape of a deviation, and returns a $p$value that is calibrated under the assumption that the reference model is exact. While Refs. [37, 38] applied the Neyman–Pearson test to the validation of generative models in high-dimensional settings, these studies did not account for their uncertainties. Ref. [15] extends the Neyman–Pearson test to include nuisance parameters in the standard LHC form [16]. We adopt this structure, but in our case the nuisance parameters are not hand-chosen systematic uncertainties, but the estimation uncertainty of the generative model.

d. Gaussian mixture models. Gaussian mixture models are a classical tool of statistics [39]. Their maximum-likelihood fitting via expectation–maximization (EM) [40, 41], the computation of the observed information matrix [42], and the asymptotic theory of mixture likelihoods [43, 44] are all well established. So are the ways mixtures depart from textbook regularity: components can be degenerate or weakly identified, and the model is identifiable only up to relabeling [45]. These issues motivate the regularization choices in Sec. IV. We claim no novelty on this machinery. Our contribution is to use the tractable score and Fisher information of a fitted mixture to extend the point-null goodness-of-fit test to a composite one, so that the uncertainty in the fit enters the test itself.

e. Goodness-of-fit with estimated parameters. Standard goodness-of-fit tests such as Kolmogorov–Smirnov [46] assume a fully specified null. Testing a null hypothesis that was itself fitted to data is a classical problem in statistics [47, 48], and the standard remedy is the parametric bootstrap [49, 50, 51]. Particle physics knows the same process as calibration with pseudo-experiments, or toy Monte Carlo [52, 53]. Our bootstrap prescription of Sec. IV D mirrors this classical procedure, combined with a composite-null test in which the fit uncertainty enters as an explicit nuisance parametrization.

Method

The method, summarized in Figure 3, has three parts: a generative model with an unconstrained parametrization (Sec. IV A), a standardized description of its estimation uncertainty (Sec. IV B), and a goodness-of-fit test in which the model may move within that uncertainty (Sec. IV C). The calibration of the goodness-of-fit test is discussed in Sec. IV D.

Generative model

The model family should be chosen to satisfy two conditions. First, we find that the uncertainty must be propagated on all of the model’s parameters, since all of them fluctuate with the training sample. Appendix B shows that approaches propagating it on only a subset, such as the mixture coefficients of an ensemble of trained normalizing flows, fail. Second, every point of parameter space must define a valid density, so that Gaussian fluctuations of the parameters never leave the model family. Gaussian mixture models (GMMs) [44] with a suitable parametrization satisfy both. For $K$ components in $D$ dimensions,

where the weights $w_k$ are given by softmax logits and the covariances, $\boldsymbol{\Sigma}_k$, by Cholesky factors with a log-parametrized diagonal,

The logits $\ell_k$, means $\boldsymbol{\mu}_k$, log-diagonals $s_{k,i}$, and off-diagonal entries $[L_k]_{ij}$ ($i>j$) are all unconstrained real numbers, and fixing one of the $\ell_k$ using the $\sum_k w_k=1$ constraint removes the invariance of the softmax under a common shift of the logits. The parameter vector $\theta$ therefore ranges over all of $\mathbb{R}^{P}$, with

The parameter count grows as $P\sim KD^2/2$, and the number of components needed to resolve a complex density also tends to grow with dimension. In high dimensions, Gaussian mixtures thus become heavily over-parametrized [54], and neural density estimators such as normalizing flows have largely replaced them [2, 55]. In low dimensions, however, they remain expressive [56] and offer what neural estimators lack: a moderate number of parameters, identifiable up to a relabeling of the components [45], with analytic scores from which the uncertainty of the fit can be computed. This makes them well suited to the low-dimensional learned embeddings considered in § 2.2. An alternative, strongly related approach with an ensemble of independently trained flows, inspired by Ref. [24], is presented in App. B.

The GMM is fitted by maximum likelihood to a training sample $\{\boldsymbol{x}_i\}_{i=1}^{N_{\mathrm{train}}}$ drawn from the unknown true density $p_0$,

using the expectation–maximization (EM) algorithm [40], the standard tool for maximum-likelihood estimation of mixtures [44]. Starting from a $k$-means initialization, EM alternates between computing the posterior probability that each training point belongs to each component (E step) and re-estimating the weights, means, and covariances in closed form with the points weighted by those probabilities (M step); each iteration increases the likelihood, and the algorithm converges to a stationary point. The GMM likelihood is highly multimodal, so the resulting $\hat{\theta}$ is in general a local maximum [43].

For small displacements $\delta\theta$ around the fit,

where

is the score function. Once the fit is done, $\psi$ is a fixed function of $x$ (see App. § 2.1.2 for its computation). The test of Sec. IV C is built in this tangent space.

Estimation uncertainty

We test the fitted model on an independent test sample $\mathcal{D}=\{x_i\}_{i=1}^{N_{\mathrm{test}}}\sim p_0$. The standard approach treats the fit as a point-null hypothesis, $H_0^{\mathrm{point}}:\mathcal{D}\sim p_{\hat{\theta}}$. But $\hat{\theta}$ is itself a random variable, which scatters around the best parameter value by $O(N_{\mathrm{train}}^{-1/2})$. The modified hypothesis, a composite-null hypothesis, $H_0^{\mathrm{comp}}$, asks whether the test data are consistent with $p_{\theta}$ for some $\theta$ close to $\hat{\theta}$ within the range of values $\theta$ could take due to random

variation in the training sample. The test statistic of Sec. IV C encodes these fluctuations through a Gaussian constraint in standardized coordinates, constructed below, and Sec. IV D describes two ways of sampling $H_0^{\mathrm{comp}}$ to calibrate it. Both hypotheses are tested against an alternative hypothesis $H_1$ which extends the null family by adding additional components and thus provides the flexibility needed to detect possible deviations.

All expectations under the fitted model are estimated with a reference set $R_0=\{x_s\}_{s=1}^{S}\sim p_{\hat{\theta}}$, with $S\gg N_{\mathrm{test}}$. It is used for the Fisher information, for the normalization of the densities in the test, and to choose the additional components of the alternative hypothesis.

If the model is correctly specified, $\hat{\theta}$ is asymptotically Gaussian with covariance $\widehat{F}^{-1}/N_{\mathrm{train}}$ [43], where $F$ is the per-event Fisher information. We estimate $F$ on the reference set from the scores of Eq. (6),

and take $\widehat{\operatorname{Cov}}(\hat{\theta})=\widehat{F}^{+}/N_{\mathrm{train}}$, where $(\cdot)^{+}$ is a pseudo-inverse with a small ridge.² For simplicity of presentation, we write estimates from $R_0$ as averages over the whole dataset, but in practice each estimate uses a subset of $R_0$ for computational efficiency (see Appendix § 2.1.2).

We diagonalize $\widehat{\operatorname{Cov}}(\widehat{\theta}) = V\Lambda V^{\top}$, with eigenvalues $\lambda_1 \ge\cdots\ge\lambda_P \ge0$, and discard the directions with $\lambda_m > \Lambda_{\max}$. These are directions the training data barely constrain [44], along which the fit can move far from $\widehat{\theta}$ and the expansion of Eq. (5) breaks down. Collecting the $M \le P$ retained eigenvectors and eigenvalues in $V_M$ and $\Lambda_M$, we define standardized coordinates $v \in\mathbb{R}^{M}$ through

which measure each parameter displacement in units of its own uncertainty. The resulting coordinates are uncorrelated and have unit variance asymptotically, $v \sim\mathcal{N}(0,I_M)$, so that a displacement of $\lVert v\rVert= 1$ corresponds to a one-sigma shift of the corresponding parameter. The corresponding standardized score features,

satisfy

where the first identity is the vanishing mean of the score and the second follows from $\widehat{F} = V_M\Lambda_M^{-1}V_M^{\top}/N_{\mathrm{train}}$ on the retained subspace.

Uncertainty-aware goodness-of-fit test

We compare data and model with a profile[2] likelihood ratio between two exponential families built on $p_{\widehat{\theta}}$.

a. Null hypothesis. To propagate the fit uncertainty into the test, the null hypothesis must be promoted from a single fixed density $p_{\widehat{\theta}}$ to a family of densities that accounts for the uncertainties; we construct this family by tilting $p_{\widehat{\theta}}$ along the uncertainty-normalized score directions. An exponential tilt is used instead of $p_{\widehat{\theta}+Wv}$ because its log-likelihood is linear in $v$, which makes the test a convex problem. The null family reads:

with the normalization $Z$ defined in Eq. (13) below. By Eq. (5), we have that $g_v = p_{\widehat{\theta}+Wv}$ to first order in $v$, so the tilt parameters are the standardized parameter displacements of Eq. (8), and $H_0^{\mathrm{comp}}$ corresponds to $v \sim\mathcal{N}(0,I_M)$.

b. Alternative hypothesis. The alternative adds additional degrees of freedom via tilts along a fixed dictionary of $J$ localized kernels, following the New Physics Learning Machine (NPLM) approach of Refs. [12, 59],

so that $g_v = h_{v,0}$. The exponential form keeps $h_{v,\alpha}$ positive for all $(v,\alpha)$. The features $\widetilde{k}_j$ are Gaussian kernels of a common width $\sigma$, $k_j(x) = \exp(-\lVert x-c_j\rVert^2/2\sigma^2)$, standardized to zero mean and unit variance on $R_0$, with centers placed by $k$-means clustering of $R_0$ so that they tile the support of the model (see Appendix § 2.1.2). The normalization in Eq. (13) is estimated on the reference set,

The dictionary should tile the support of the model densely enough that no localized departure falls between kernels; we use $J = 500$. Beyond this point additional kernels add little: neighboring kernels of a common width overlap strongly, so the dictionary is highly collinear, and the ridge penalty in Eq. (15) suppresses the redundant combinations. The effective number of degrees of freedom of the test is therefore set jointly by the kernel width $\sigma$ and the ridge rather than by $J$ itself (see Appendix § 2.1.3), and our results are stable under moderate changes of $J$.

The width $\sigma$ sets the length scale on which the augmented model Eq. (12) is sensitive to departures from the null family, and its choice is bounded on both sides. From above, to ensure power, $\sigma$ should not exceed the scale of the discrepancies the test is meant to detect, which for the smooth densities considered here is set by the length scale of the data. From below, to ensure calibration, $\sigma$ is bounded by the validity of the calibration procedures to be described in Section IV D (this point is discussed further in Appendix § 2.9). To strike a balance between these considerations, we take $\sigma$ to be set by the characteristic scale of the data (alternatively, tests at several widths can be combined [60]); see Appendix § 2.1.2 for details.

c. Test statistic. On the test sample, the log-likelihood ratio of $h_{v,\alpha}$ to $p_{\widehat{\theta}}$, together with the Gaussian constraint on $v$ and a ridge penalty on $\alpha$, is

stead given by the sandwich formula [57]. We do not use it, since the null distribution is computed assuming the null hypothesis, under which the two estimates coincide. The sandwich then differs from $\widehat{F}^{+}/N_{\mathrm{train}}$ only when the null is false. This choice is also consistent with the parametric bootstrap of Section IV D, which intrinsically assumes well-specification.

Treating the $v$ as nuisance parameters with a Gaussian constraint [16], the test statistic is the profile likelihood ratio

The denominator lets the model move within its uncertainty to accommodate the data, and the numerator additionally switches on the localized distortions. The statistic $t_{\mathrm{comp}}$ therefore measures only the incompatibility that the model’s uncertainty cannot explain. As an uncertainty-blind baseline, the point-null statistic fixes $v=0$,

with no subtraction needed since $\ell(0,0)=0$. Because $\ell(v,0)$ is the restriction of $\ell(v,\alpha)$ to $\alpha=0$, both statistics are nonnegative. For numerical safety, the maximizations over $\alpha$ are also restricted to a box that bounds the log-density distortion of each kernel (see App. § 2.1.2 for details). Every maximization above is a convex problem. In Eq. (15), the data term is linear in $(v,\alpha)$, the penalties are concave quadratics, and

is a log-sum-exp of linear functions, and hence convex; the box constraint defines a convex set. The objective $\ell$ is therefore strictly concave, and each maximization has a unique optimum, which a quasi-Newton method reaches promptly.

d. Additional nuisance parameters. The construction above is the standard profile-likelihood treatment of constrained nuisance parameters [16], with $v$ encoding the statistical uncertainty of the fit. Nothing restricts the nuisance vector to this role: experimental and theoretical systematic uncertainties, parametrized as additional constrained directions in Eq. (15), can be profiled jointly with $v$, so that a single composite test accounts for all sources of uncertainty on the reference model. This holds provided the systematic variations admit a density parametrization, as in standard template-morphing approaches. We leave this combination to future work and focus here on the fit uncertainty alone, which is the component that existing validation practice omits.

Calibration

The null distribution of $t_{\mathrm{comp}}$ is obtained from $N_B$ pseudo-experiments that sample $H_0^{\mathrm{comp}}$. For each, the statistic $t_b$ is computed with the full procedure above, including the same penalties and constraints, and the $p$-value is the rank of the observed statistic,

We consider two ways of generating the pseudo-experiments. The statistic is the same in both; only its calibration differs.

a. Gaussian prescription. If $\hat{\theta}$ is Gaussian with the covariance of Sec. IV B, $H_0^{\mathrm{comp}}$ can be sampled directly,

For each pseudo-experiment, a value of $\theta_b$ is chosen, and pseudo-data are then drawn from the GMM at $\theta_b$. The test statistic $t_b$ is then computed on this pseudo-data in the same way as it is for the observed data. This procedure suffices to calibrate the observed statistic if the asymptotic approximation is reliable and the Gaussian distribution of $\theta_b$ reproduces the sampling distribution of $\hat{\theta}$.⁴

b. Bootstrap prescription. For a GMM with many parameters and a multimodal likelihood, the Gaussian description of $\hat{\theta}$ can fail: expectation maximization (EM) may converge to different local maxima for different training samples [43], and the scatter of the fit need not be Gaussian, nor small enough for Eq. (5) to hold. The parametric bootstrap relaxes this assumption by directly simulating the fitting procedure. Treating $p_{\hat{\theta}}$ as if it were the true density, we repeatedly draw training samples from it and refit; the spread of the refits $\hat{\theta}_b$ around $\hat{\theta}$ then serves as a proxy for the (unobservable) spread of $\hat{\theta}$ around the truth. The necessary assumption is that the spread from the fitting procedure is similar whether the data come from $p_0$ or from the nearby $p_{\hat{\theta}}$, which is weaker than the asymptotic conditions required for the Gaussian prescription. However, this assumption can fail even for a well-specified model, and in particular we find in App. § 2.9 that the bootstrap calibration breaks down if the kernel width of Sec. IV C is chosen too narrow compared to the characteristic scale of the data.

In practice, for each pseudo-experiment, we draw a training sample of size $N_{\mathrm{train}}$ from $p_{\hat{\theta}}$ and refit the full GMM with EM from a fresh $k$-means initialization, obtaining $\hat{\theta}_b$. We rebuild the entire construction at $\hat{\theta}_b$ exactly as for the nominal fit (i.e. the reference set, the Fisher information, the standardization, and the kernel dictionary), draw a test sample of size $N_{\mathrm{test}}$ from $p_{\hat{\theta}_b}$, and compute $t_b$.

These two prescriptions for generating the calibration distribution are, in a sense, mirror images of each other. The Gaussian prescription holds the fit fixed at $\hat{\theta}$ and moves the data-generating model to $\theta_b$, while the bootstrap prescription moves the fit to $\hat{\theta}_b$ and holds the data-generating model fixed at $\hat{\theta}$. Both, then, should reproduce the displacement between fit and truth, but the bootstrap does not assume that this displacement is small or Gaussian, and it captures the optimization variability. We will find that it remains calibrated to larger $N_{\mathrm{test}}$ (see § 2.2) and is our nominal prescription. The Gaussian prescription requires no refits and provides an explicit parametrization of the model uncertainty that can be easily propagated to the density and to downstream applications, and it remains useful where its assumptions hold. The two prescriptions also differ in their exposure to the truncation choice: $\Lambda_{\max}$ enters the Gaussian prescription twice, through the profiled directions in the test and through the fluctuations injected in the pseudo-experiments, whereas the bootstrap refits are constructed independently of it and dependence appears only in the test. As such, scanning $\Lambda_{\max}$, we find in App. § 2.1.2 that the bootstrap calibration is robust against the choice of $\Lambda_{\max}$, while the Gaussian prescription shows a stronger dependence.

In terms of the sources of Sec. II, the Gaussian prescription captures the data variability around a single optimum, while the bootstrap, which refits the model from a fresh initialization, also captures the optimization variability. This is why, in practice, we find that the bootstrap prescription yields uncertainty bands that are comparable to or larger than those of the Gaussian prescription. Neither, however, accounts for the model bias.

For the point-null baseline, the pseudo-experiments are simply drawn from $p_{\hat{\theta}}$ with the construction held fixed.

c. Asymptotic $\chi^2$ limit. In the quadratic approximation to $\ell$ and without the ridge penalty $\lambda_{\alpha}$ on the $\alpha$, $t_{\mathrm{comp}}$ follows a $\chi^2_J$ distribution under the Gaussian composite-null, as shown in App. § 2.1.3. The same holds for $t_{\mathrm{point}}$ under the point null, and the ridge lowers the effective number of degrees of freedom below $J$. We will see in § 2.2 that this approximate $\chi^2_J$ prediction holds (up to a decrease in the effective degrees of freedom due to $\lambda_{\alpha}$) in the regime where the asymptotics are reliable and the Gaussian and bootstrap prescriptions are in good agreement, but that it breaks down outside of this regime (as it should, since the derivation of this $\chi^2_J$ distribution assumes Gaussianity). As such, our results for test calibration are always evaluated using pseudo-experiments.

It is worth noting that the test described in Sections IV C and IV D is independent of how the density is estimated and of how its uncertainty is obtained, and it works provided this uncertainty is reliably estimated and can be propagated into the composite test.

Numerical experiments

In this section, we apply the method presented in Sec. IV to two problems: a two-dimensional toy model with known analytic ground truth and a four-dimensional problem from particle physics.

Toy model

We first validate the pipeline on a controlled two-dimensional toy problem, where the ground truth is known analytically and can be resampled at will. The target factorizes into two independent features, each chosen to stress a different aspect of the fit:

• Feature 1: a bimodal Gaussian mixture, $p(x_1) = \frac{1}{2}\mathcal{N}(x_1; -0.70, 0.12^2) + \frac{1}{2}\mathcal{N}(x_1; -0.30, 0.12^2)$, whose modes are separated by more than three times their common width, with analytic mean $\mu_1^* = -0.50$.

• Feature 2: a skew-normal distribution with location $\xi= 1.0$, scale $\omega= 0.75$ and shape $\alpha= 8.0$, i.e. strong positive skew, with analytic mean $\mu_2^* = \xi+ \omega\delta\sqrt{2/\pi}$ and $\delta= \alpha/\sqrt{1+\alpha^2}$.

The first feature probes the ability of the model to resolve well-separated modes, the second its ability to reproduce a long, asymmetric tail.

We draw $N_{\mathrm{train}} = 10^5$ training points and fit the Gaussian mixture model of Eqs. (1) and (2) for four choices of the number of components, $K = 16, 32, 64, 128$. We then estimate the model uncertainties following Sec. IV B and assess the quality of the fit with the point-null test of Eq. (17), which treats the fitted density as exact, and the composite test of Eq. (16), which accounts for the uncertainty of the fit. For the composite test we report both our default bootstrap calibration as well as the Gaussian prescription of Sec. IV D. Because the ground truth is available in closed form, a fresh test sample is drawn for every replica, so the closure test is not subject to reuse of validation data across toys.

Figure 4 shows the marginal distributions of the two features together with a $K = 64$ fit. The fitted model (black solid line) is compared with the empirical distribution of a held-out sample the same size as the training set ($\times$ markers), with the data-to-fit ratio shown in the bottom panels. The 68% bands are built under the Gaussian prescription (light blue) and under the bootstrap prescription (light red). The Gaussian and bootstrap bands are overall comparable, with the bootstrap band slightly larger in regions where the Gaussian approximation underestimates the spread of the fit around the truth.

Figure 5 shows the median $p$-value for the point null as well as the Gaussian and bootstrap composite nulls as a function of the test-set size $N_{\mathrm{test}}$ for the four values of $K$. For each curve, a median $p$-value of 0.5 would indicate correct calibration, since under the null hypothesis the $p$-value is uniformly distributed by construction. We perform 20 independent model fits for each value of $K$, each evaluated on 100 test datasets per value of $N_{\mathrm{test}}$. The point-null and Gaussian composite statistics are calibrated with 250 pseudo-experiments per fit, while the bootstrap calibration uses 100 refits in total per fit, with the same refits reused across all values of $N_{\mathrm{test}}$.

The point-null test is calibrated only in the small-sample regime. Its median $p$-value is compatible with 0.5 for $N_{\mathrm{test}} \le N_{\mathrm{train}}/10$, starts degrading beyond that point, and collapses to zero by $N_{\mathrm{test}} = N_{\mathrm{train}}$ for every value of $K$ we consider. The behavior is essentially independent of the number of mixture components $K$, as the test fails because it neglects the uncertainty of the fit, regardless of the flexibility of the model.

Promoting the test to the composite-null, with the model uncertainties estimated under the Gaussian prescription, largely restores calibration at intermediate values of $N_{\mathrm{test}}$. The median $p$-value stays compatible with 0.5 up to about $N_{\mathrm{test}} = N_{\mathrm{train}}/5$, but calibration still fails dramatically at $N_{\mathrm{test}} = 2N_{\mathrm{train}}$, with a median $p$-value less than 0.2 for all $K$.

The bootstrap calibration prescription fully restores closure, remaining calibrated up to $N_{\mathrm{test}} = 2N_{\mathrm{train}}$ for all $K$. This is consistent with the bootstrap capturing the optimization variability and the non-Gaussian scatter of the fit, which the Gaussian description misses. These results demonstrate that the point-null hypothesis is insufficient. We find that closure is dramatically improved through the use of a composite-null hypothesis, and that maintaining closure for $N_{\mathrm{test}} \gtrsim N_{\mathrm{train}}$ requires bootstrap (rather than Gaussian) calibration. Figure 6 illustrates this directly for $K = 64$ at $N_{\mathrm{test}} = N_{\mathrm{train}}$: the observed point-null test statistic lies almost entirely above its null distribution, while the composite-null statistic is well reproduced by its bootstrap null.

Embedded particle jets

Recent advances in representation learning have made it possible to construct powerful neural embeddings of highly complex physics data. Among other applications, these embeddings have enabled anomaly detection in a latent space [61, 62, 63, 64, 65, 66]. Anomaly detection in a latent space, however, requires a robust understanding of the density of the reference model. The method proposed in this paper offers a path to estimate a background data distribution and its uncertainties with statistical guarantees.

As a demonstration, we emulate a realistic physics analysis aiming to identify anomalous jet decays. We use the JetClass dataset [67], a large-scale public benchmark for jet tagging. The JetClass dataset provides 100M simulated jets for training, 5M for validation, and 20M for testing. Each jet is represented as a particle cloud of charged and neutral hadrons, electrons, muons, and photons and labeled by the initiating parton and decay process, with ten equally represented classes. We restrict the analysis to purely hadronic final states and consider five classes: decay of top quarks, W and Z bosons, Higgs bosons and QCD multijet production.

Following Ref. [65], the embedding is trained with a SupCon (contrastive) objective [68], which extends self-supervised contrastive learning [69] to labeled data: all jets sharing a class label are treated as positives, and jets from different classes as negatives. For training the embedding, we restrict the JetClass training and validation set to 10% of the available statistics, corresponding to 5M and 250k jets, respectively, equally distributed among the five considered classes. The resulting representation is shown in Figure 7.

Anomaly detection via two-sample testing [12, 13, 14] requires either a large reference (simulated) sample or an equivalent of the underlying probability density function in the embedding space. Since large simulated samples can be prohibitively expensive to generate, we apply the SIGMA method as described in Sec. IV to model the QCD background in the embedding space.⁵ We choose QCD multijet production because it is a dominant background whose simulated statistics are hard to match to the large number of events expected in LHC data. We apply our pipeline to the full JetClass QCD multijet test partition with $N_{\mathrm{test}}^{\mathrm{TOT}} = 2\mathrm{M}$ jets, which we split into a training pool and a testing pool each consisting of 1M jets. The dataset can be found on Zenodo [70].

Each SIGMA training and $p$-value evaluation is performed by sampling from these finite pools, so we report closure results only up to $N_{\mathrm{test}} = N_{\mathrm{train}} = 10^5$ events to mitigate correlations between different trainings and draws of pseudo-observed data.

Figure 8 shows the marginal distributions of the four latent features of our embedding for the QCD sample, together with a $K = 256$ fit. As for the two-dimensional problem, we compare the fitted model (black solid line) with the empirical distribution of a same-size held-out sample ($N = 10^5$) marked by the “×” symbol. The data-to-fit ratio is shown in the bottom panels. The 68% bands are built under the Gaussian prescription (light blue) and under the bootstrap prescription (light red).

Figure 9 shows the median $p$-value for the goodness-of-fit test using the point-null, the Gaussian composite and the bootstrap composite as a function of $N_{\mathrm{test}}$ for $K = 128, 256, 512$. For each value of $K$, we perform 10 independent model fits, each evaluated on 100 test datasets per value of $N_{\mathrm{test}}$. The point-null and Gaussian composite statistics are calibrated with 250 pseudo-experiments per fit; since the bootstrap calibration is considerably more expensive, it uses 100 refits in total per fit, with the same refits reused across all values of $N_{\mathrm{test}}$.

The results are qualitatively similar to those in 2D. At all values of $K$, the point-null test starts degrading before $N_{\mathrm{test}} \sim N_{\mathrm{train}}/10$. The composite test with the Gaussian prescription improves the calibration, enabling closure up to $N_{\mathrm{test}} \sim N_{\mathrm{train}}/10$ for $K \ge256$. The bootstrap prescription further improves the results allowing closure up to $N_{\mathrm{test}} = N_{\mathrm{train}}/5$ for $K \ge256$, with calibration degrading noticeably (but not completely) for larger $N_{\mathrm{test}}$.

The $K = 128$ SIGMA model performs visibly worse than either of the two larger models, and the difference between the Gaussian and bootstrap null prescriptions is substantially smaller than for larger $K$. We interpret this as underparameterization-induced model bias in the sense of Sec. II, where a mixture model with only 128 components is insufficiently expressive to fit the data. This model misspecification is not quantified in our statistical uncertainties, so sufficiently large $N_{\mathrm{test}}$ necessarily rejects this model.

By contrast, the calibration is evidently saturating with $K$ for large $N_{\mathrm{test}}$ for the $K = 256, 512$ models, with the larger difference between the bootstrap and Gaussian prescriptions indicating a breakdown in the asymptotics. This suggests that the finite $N_{\mathrm{train}} = 10^5$, rather than model capacity, is likely the primary barrier to achieving a more performant SIGMA in the context of this case study, as could be expected from fitting a model with several thousand parameters. However, even with these limitations, the bootstrap composite null is calibrated to values of $N_{\mathrm{test}}$ an order of magnitude larger than the point-null.

To characterize the extent to which these results are specific to SIGMA models, we also test parametrized ensembles of normalizing flows,⁶ and the results can be found in App. B 4. In short, we find that the SIGMA models achieve substantially better calibration properties than even ensembles of $\mathcal{O}(10^{2})$ flows, while being dramatically cheaper to train and evaluate.

Conclusion

All generative models trained on finite samples are inexact approximations of their training distribution, and a sufficiently powerful goodness-of-fit test will eventually detect deviations. This is a consequence of finite information, not a failure of design. A density estimated from $N_{\mathrm{train}}$ samples carries an inherent estimation error, so even a perfectly specified model is rejected by a powerful goodness-of-fit test. Validating a generative model against a point-null therefore measures the size of the test sample rather than the quality of the model, and no architectural or training improvement can change this.

An approximate model can nonetheless be useful, provided its estimation error is known. In this work, we make this statement quantitative. By linearizing a Gaussian mixture model around its best fit, we obtain the covariance of its parameters from its analytic score functions and Fisher information, and the goodness-of-fit test becomes a convex problem. We find that promoting the null from a point hypothesis to a composite hypothesis, together with a parametric bootstrap calibration procedure, is able to calibrate the goodness-of-fit test on datasets an order of magnitude larger than the point-null equivalent.

The success of SIGMA is in contrast with neural generative models trained through stochastic gradient descent, which typically do not converge to a well-defined optimum and which have non-invertible information matrices due to parameter degeneracies. We find in App. B that even a parametrized ensemble of normalizing flows is only able to capture part of the uncertainty, and cannot match the performance of SIGMA in our low-dimensional case studies.

Uncertainty quantification for generative models bears directly on how generated samples can be used. A generative model cannot create information that was not in its training data [10]. What it can do is produce arbitrarily many samples whose deviation from the truth is bounded by stated uncertainties. Generated samples can therefore be used beyond the training statistics only if these uncertainties are carried along into the downstream analysis, e.g. by using a composite-null hypothesis in the context of a goodness-of-fit or hypothesis test. For statistical inference, the true figure of merit of a generative surrogate is therefore the precision to which it can be trusted.

composite-null prescriptions, as the bootstrap-null would require retraining the whole ensemble of flows in each pseudo-experiment, which would be prohibitively expensive.

Within this framework, SIGMA still has clear limitations. Because the number of parameters of a GMM grows quickly with dimension, the method is suited to low-dimensional spaces rather than directly to high-dimensional ones. This further motivates the compression of high-dimensional data into a small number of summary statistics, such as through learned embeddings [61, 62, 63, 64, 65, 66].

Furthermore, SIGMA accounts for the statistical variability of the fit due to the finite training sample and to the stochasticity of optimization, but it is still subject to model bias: a model family that cannot represent the data will still be rejected by a sufficiently powerful test. The need to reduce model bias makes the choice of $K$ a trade-off. Increasing the number of components reduces mismodeling but it also adds parameters, leading to larger uncertainties and potentially undermining the calibration procedures. The Gaussian calibration manifestly requires the asymptotic regime where $N_{\mathrm{train}} \gg P$, and the bootstrap calibration relies on this indirectly through the assumption that the estimator spread is similar for fits to draws from $p_{0}$ and from $p_{\hat{\theta}}$.

We currently choose $K$, as well as the truncation hyperparameter $\Lambda_{\max}$, empirically. These model selection choices, as well as verifying that the model can be trusted beyond its training statistics, require closure tests on held-out samples at least as large as the desired inference dataset. This is costly in exactly the simulation-scarce settings where generative models are most useful, and the necessity of the bootstrap for the most trustworthy calibration adds the cost of refitting the model for every pseudo-experiment.

Each of these limitations offers promising avenues for progress. Both the restriction to low-dimensionality and the model bias considerations motivate exploring distinct model classes. In particular, the test of Sec. IV C is agnostic to the model class: it requires only tractable score vectors, which automatic differentiation supplies for any neural density estimator, and a faithful estimate of the uncertainty of the fit. Our tangent-space construction is then the exact, low-dimensional counterpart of the linearized Laplace approximations developed for deep networks [71, 72], and that literature’s low-rank and subspace curvature estimates offer a natural route to extending the composite-null test to larger models; the open problem is promoting such uncertainties from predictive heuristics to calibrated ingredients of a goodness-of-fit test.

On the other hand, the need for empirical tuning of $K$ and $\Lambda_{\max}$ could be eliminated with a data-driven procedure for hyperparameter selection. Such a procedure could involve the use of cross-validation to check closure on subsets of held-out simulation data.

Moreover, it would be interesting to consider applications with systematic uncertainties. Because the fit uncertainty enters as a standard constrained nuisance vector, the Gaussian calibration prescription composes naturally with the existing treatment of systematic uncertainties [15, 73]: additional nuisance directions are profiled jointly, and a single composite test accounts for all sources of uncertainty on the reference model. In contrast, composing our nominal bootstrap calibration prescription (which does not yield straightforwardly parametrized uncertainties) with additional nuisance parameters may be more challenging.

The immediate application of our method is background estimation for anomaly detection in low-dimensional learned embeddings. There, a validated model could enter the final test in closed form, together with its uncertainties, as the composite-null hypothesis. More broadly, as generative models increasingly substitute for simulation across various scientific fields, uncertainty-aware validation is necessary for these surrogates to serve as instruments of trustworthy statistical inference.

The code required to reproduce our experiments and implement SIGMA can be found here: https://github.com/gbadarac/UQ-GoF-GenModels/tree/paper.

Acknowledgments

We are especially grateful to Thea Årrestad for her supervision, guidance and support in the development of this work. We also thank Davide Valsecchi for discussions on uncertainty quantification for generative models and for his detailed feedback on the manuscript. We further thank Prasanth Shyamsundar for conversations on the information content of generated samples and the limits of data amplification. The authors acknowledge the financial support of the National Science Foundation under Cooperative Agreement PHY-2019786 (NSF AI Institute for Artificial Intelligence and Fundamental Interactions, http://iaifi.org/), and the National Science Foundation AI Expand program (OAC-2334265). C.R. further acknowledges support under a Postdoc.Mobility fellowship from the Swiss National Science Foundation (SNSF) (Project No. P500PT 222340). Claude (Anthropic) assisted with manuscript editing and code development. All scientific ideas and experiments are the original work of the authors.

Additional discussion on the limits of generative models

In this appendix, we quantify the statistical uncertainty of the fit introduced in Sec. II in terms of the Kullback–Leibler (KL) divergence between the true and the fitted densities, and derive the test sample size at which it becomes detectable by a goodness-of-fit test. We consider a parametric family $p_\theta$ with $\theta\in\mathbb{R}^{P}$ and assume the best-case scenario of Sec. II: the model is well specified, $p_0 = p_{\theta_0}$ for some interior point $\theta_0$ with finite, non-singular Fisher information $F(\theta_0)$, and the optimizer returns the global maximum-likelihood estimate $\hat{\theta}$ computed on $N_{\mathrm{train}}$ independent samples of $p_0$ (i.e. we isolate the data variability; any optimization variability would only add to the result below).

Expanding $\theta\mapsto\operatorname{KL}(p_0 \mathbin{\|} p_\theta)$ around $\theta_0$, the zeroth-order term vanishes because $p_0 = p_{\theta_0}$, and the gradient vanishes because $\theta_0$ minimizes the divergence over the family. The curvature at the minimum is the Fisher information: since the $\theta$-dependence of the divergence sits entirely in the cross-entropy $-\mathbb{E}_{p_0}[\log p_\theta]$, its Hessian is

To leading order, therefore,

the divergence of the fitted density is a quadratic form in the estimation error, with the information matrix as its metric. The maximum-likelihood estimate asymptotically fluctuates around the true value as $\hat{\theta}-\theta_0 \sim\mathcal{N}(0,F(\theta_0)^{-1}/N_{\mathrm{train}})$, and for a Gaussian vector $\mathbf{u}$ one has $\mathbb{E}[\mathbf{u}^{\mathsf{T}}A\mathbf{u}] = \operatorname{tr}(A\Sigma)$. Averaging Eq. (22) over the training sample then gives, for large enough $N_{\mathrm{train}}$,

This is the data variability, expressed in terms of the KL divergence. Note that the Fisher information cancels against its own inverse: every identifiable direction in parameter space contributes $1/(2N_{\mathrm{train}})$ nats to the expected divergence, independently of how tightly the data constrain it. A poorly measured direction fluctuates more, but the divergence is correspondingly flat along it, and the two effects compensate. Only the number of effective directions, $P$, survives. Equivalently, $2N_{\mathrm{train}}\operatorname{KL}(p_0 \mathbin{\|} p_{\hat{\theta}})$ is asymptotically the Wald statistic $N_{\mathrm{train}}(\hat{\theta}-\theta_0)^{\mathsf{T}}F(\hat{\theta}-\theta_0)$, which, like the likelihood-ratio statistic of Wilks’ theorem, is $\chi^2_P$-distributed with mean $P$. For the Gaussian example of Figure 2, $P=2$, and indeed $\operatorname{KL}\simeq\hat{\mu}^2/2 + (\hat{\sigma}-1)^2$ with $\operatorname{Var}(\hat{\mu})=1/N_{\mathrm{train}}$ and $\operatorname{Var}(\hat{\sigma})=1/(2N_{\mathrm{train}})$, so that $\mathbb{E}[\operatorname{KL}]\simeq1/N_{\mathrm{train}}$.

A goodness-of-fit test performed with $N_{\mathrm{test}}$ samples becomes sensitive to a discrepancy of size $D$ (in KL) once $N_{\mathrm{test}}D \gtrsim O(1)$. This is the Neyman–Pearson benchmark: for the most powerful test of $p_{\hat{\theta}}$ against the specific alternative $p_0$, the expected log-likelihood ratio per test sample is $\mathbb{E}_{p_0}[\log(p_0/p_{\hat{\theta}})] = \operatorname{KL}(p_0 \mathbin{\|} p_{\hat{\theta}}) = D$, so the evidence accumulated over the test set, $N_{\mathrm{test}}D$, must exceed its $O(1)$ null fluctuations. Substituting Eq. (23), which represents the discrepancy that the most powerful test would detect, this threshold is crossed once

This argument should be understood as establishing a sensitivity ceiling for goodness-of-fit tests, describing when the training-set discrepancy becomes detectable in principle, not necessarily when a specific test statistic used in practice (such as the test statistics described in Sec. IV C) rejects. The actual rejection point of a test depends on its power against this particular alternative and on its empirical calibration, and can in general occur later, and with a different $O(1)$ prefactor, than this bound suggests: a statistic built on $J$ degrees of freedom, for instance, resolves the same $D$ only once $N_{\rm test}D \gtrsim\sqrt{J/2}$, which rescales Eq. (24) by a constant without affecting the scaling in $N_{\rm train}$ and $P$.

Two qualifications delimit the validity of Eq. (23) itself. First, it is the leading term of an asymptotic expansion: it requires $N_{\rm train}$ large enough for the quadratic approximation (A2) to hold across the sampling distribution of $\hat{\theta}$, which fails when $P$ is not small compared to $N_{\rm train}$ or when some directions become non-identifiable (e.g. degenerate mixture components); in that regime the effective number of parameters is smaller than $P$. Second, under misspecification the same expansion around the KL projection $\theta^*$ of $p_0$ onto the family yields a variance term $\operatorname{tr}(A^{-1}B)/(2N_{\rm train})$, with $A$ the expected curvature and $B$ the score covariance ($A=B=F$ for a correct model), on top of the bias floor $\operatorname{KL}(p_0 \mathbin{\|} p_{\theta^*})$, which does not shrink with $N_{\rm train}$. The floor eventually dominates the discrepancy, and the rejection point saturates at $N_{\rm test}\sim1/\operatorname{KL}(p_0 \mathbin{\|} p_{\theta^*})$ independently of further growth in $N_{\rm train}$: this is the bias that uncertainties cannot absorb and that must be addressed at the level of model design.

These considerations sharpen the usage criteria of Figure 1 [11]. When the fitted density enters an inference as the reference hypothesis, its estimation error may be neglected only while it lies below the sensitivity of that inference, and the validation test used to establish this must be at least as powerful as the downstream analysis (see the discussion of Sec. III); otherwise, mismodeling that validation cannot see may still surface downstream as a false positive. Equation (24) bounds how long such safety can last for any analysis: beyond $N_{\rm test}\gtrsim N_{\rm train}/P$, even a perfectly specified model is rejected when tested, and the uncertainties of the fit must be incorporated as a composite null, as described in Sec. IV C.

Normalizing flows and the limits of ensemble uncertainties

Normalizing flows [55] are among the most widely used and flexible neural density estimators in high-energy physics, both as surrogates for expensive simulation [2, 31] and as background models [4, 5, 6, 7, 25], including density estimation in learned embedding spaces for anomaly detection searches [66]. Unlike the Gaussian mixtures of the main text, a flow has no tractable parameter covariance, so the score-based, full-parameter-space construction of Sec. IV is unavailable.

In this appendix, we consider the same two numerical experiments from § 2.2, a two-dimensional toy with a bimodal Gaussian and a skew-normal feature and a four-dimensional latent space for jet representations, and we ask the following two questions: (1) does the epistemic uncertainty of a normalizing flow ensemble admit a parametrization that a likelihood-ratio test can profile over? and (2) with this uncertainty propagated in the null, does a normalizing flow ensemble pass a composite goodness-of-fit test such as the one presented in the main text?

Ref. [74] provides such a parametrization through the maximum-likelihood $w_i f_i$ ensembling of Ref. [24], in which independently trained flows are frozen and combined into a weighted mixture whose weights carry the uncertainty. This appendix extends that work in two directions: to more complex target distributions, and to a goodness-of-fit test built on the same construction as the main text.

We describe how we train the normalizing flow models in Sec. B1, how we estimate uncertainties with the $w_i f_i$ ensembling in Sec. B2, and how we incorporate the uncertainties in the goodness-of-fit test in Sec. B3. The results of our numerical experiments, presented in Sec. B4, show that a normalizing-flow ensemble admits such a weight parametrization, and that propagating it restores closure up to and beyond $N_{\rm test}=N_{\rm train}$ for the largest ensemble ($M=128$) in the two-dimensional toy, but not in the four-dimensional embedding, where the propagated weight uncertainty is negligible compared to the flow’s misspecification: the composite test collapses onto the point-null and the model is rejected well before $N_{\rm test}=N_{\rm train}$.

Density estimation with normalizing flows

A normalizing flow learns an invertible, differentiable map $T_\phi$ from a simple base distribution to the target,

and is trained with a maximum likelihood objective. We use nflows [75] with masked rational-quadratic autoregressive layers, reverse permutations between layers, dropout 0.2, and a standard normal base distribution. The architecture is chosen separately for each dataset by a grid search on held-out likelihood: four layers (8 bins, 4 blocks, hidden width 64) for the two-dimensional toy, and a larger flow with six layers (10 bins, 8 blocks, hidden width 128) for the four-dimensional jet embedding. The training settings are the same for both datasets and were kept fixed in both grid searches: Adam [76] with batch size 512 and an initial learning rate of $10^{-5}$, cosine-annealed to $10^{-6}$ over 1000 epochs, with early stopping (patience 10 epochs) on a 20% validation split of each member’s training sample. In both cases training ends by early stopping, typically after $\sim100$ epochs, well before the 1000-epoch limit.

For each of the two problems (the two-dimensional toy and the four-dimensional jet embedding), we train a pool of 512 flows, each an independent run on its own bootstrap replica of the $N_{\mathrm{train}} = 10^5$ training set, so the spread across members reflects both data and optimization variability. An ensemble of size $M$ is built by drawing $M$ members uniformly at random from this pool, and we show $M = 1, 32, 128$. Both the fitted density and its propagated weight uncertainty of Eq. (29) depend on which members happen to be drawn, so a single draw gives a noisy estimate of the uncertainty. To remove this dependence, we repeat the draw ten times, obtaining ten independent ensemble choices of size $M$, each fit independently (Sec. B 2). In the marginal figures below the fitted model is the median over the ten choices and the band is the per-choice weight uncertainty averaged over them. This same per-choice uncertainty is what the goodness-of-fit test profiles over: it enters the null as the Gaussian prior $\pi(\boldsymbol{u}) = \mathcal{N}(\hat{\boldsymbol{u}}, \Sigma_u)$ of Sec. B 3, so the band shown here is exactly the uncertainty propagated into the test. Each of the ten choices is then tested and calibrated independently, and the ten results are pooled into the reported $p$-value and its band (Sec. B 4).

Uncertainty quantification via $w_i f_i$ ensembling

Following [24], we freeze the trained flows $f_i(x)$ as fixed basis functions and place the uncertainty on a small set of mixture weights $w_i$. The ensemble density is

where the sum-to-one constraint preserves normalization and leaves $M-1$ free weights $\boldsymbol{u} = (w_1,\ldots,w_{M-1})$. The underlying assumption is that the true density lies in the span of $\{f_i\}$. The weights are signed ($w_i \ge0$ is not imposed because allowing cancellations between the networks improves expressivity) and are fit by minimizing the negative log-likelihood on the training sample,

which is convex in $\boldsymbol{u}$ on the feasible set $\{\boldsymbol{u}: g(x_n)>0\ \forall n\}$. We find $\hat{\boldsymbol{u}}$ by damped Newton steps with a feasibility line search. The resulting estimator $\hat{g}(x)=g(x;\hat{\boldsymbol{u}})$ is the best estimate of the ensemble density.

Treating $\hat{\boldsymbol{u}}$ as an M-estimator [77], its covariance is given by the sandwich formula [57]

Here $\mathcal{L}$ is the mean negative log-likelihood of Eq. (27), $A$ its Hessian at the optimum, and $B$ the covariance of the per-event scores $\boldsymbol{s}_n$, whose components are the weight-difference sensitivities $s_{n,j}=-(f_j(x_n)-f_M(x_n))/g(x_n)$. We denote this $(M-1)\times(M-1)$ covariance by $\Sigma_u=\operatorname{Cov}(\hat{\boldsymbol{u}})$. Propagating it to the density through the same sensitivities gives the pointwise predictive variance

i.e. an uncertainty band $\hat{g}(x)\pm\sigma_g(x)$ that supplies the pointwise epistemic uncertainty later propagated into the goodness-of-fit test of Sec. B 3.⁷

Equation (B4) is the exact analogue of the sandwich covariance discussed in the main text,⁸ but restricted to the $M-1$ mixture weights: it captures only re-weightings $\sum_i\delta w_i f_i(x)$ of the fixed members, not deformations of any individual flow. As shown in Sec. B 4, this is the reason why the composite test built on it approaches the uniform null $p$-value only in the 2D toy, and fails to improve on the point-null in the 4D embedding.

Uncertainty-aware goodness-of-fit test

We test the ensemble density with the composite-null likelihood-ratio construction of the main text (Sec. IV C), now applied with the mixture weights as the profiled parameters.

a. Null. The null is the family of re-weighted mixtures $g$ of Eq. (26), with $\boldsymbol{u}$ constrained by the Gaussian prior centered on the fit,

The trainable parameters $\boldsymbol{u}$ are profiled under this prior when the statistic is evaluated.

b. Alternative. The alternative enriches the null by a multiplicative exponential tilt along a fixed dictionary of localized kernels, as in the main text [59],

with Gaussian kernels $k_j(x) = (2\pi\sigma^2)^{-D/2}\exp(-\lVert x-c_j\rVert^2/2\sigma^2)$. The kernel dictionary is similar to that of the main text (Sec. IV C), using $J$ Gaussian kernels of a single width $\sigma$ (values in Sec. B 4), but in this case centered on the first $J$ points of the sample under test. Since the $J=500$ kernels are each about as wide as the spread of the data and overlap heavily, so every region populated by the data is covered whichever way the centers are placed, and we expect the choice of placement to have no significant effect on the test. Setting $\alpha=0$ recovers the null ($h=g$), so the null is nested in the alternative. The exponential tilt preserves the sign of $g_u$, and $h$ is normalized up to the Monte Carlo estimate $\widehat{Z}$ of Eq. (33).

As with the test in the main text, two regularizations act on the tilt: an $L_2$ ridge $\frac{1}{2}\lambda\lVert\alpha\rVert^2$ and a box $\lvert\alpha_j\rvert\le A$ on the tilt coefficients. The normalization (B8) is estimated by importance sampling from the equal-weight member mixture $q(x)=\frac{1}{M}\sum_m f_m(x)$,

with $y_r\sim q$ and $R=N_{\rm ref}=N_{\rm test}$ reference points. The reference weights are self-normalized so that $\widehat{Z}(u,0)=1$.

c. Test statistic. The test statistic follows the construction of the main text (Sec. IV C), with the mixture weights $u$ in place of the whitened coordinates $v$ and the prior $\pi(u)$ in place of the Gaussian constraint on $v$. Following the standard LHC nuisance-parameter prescription [16], we define the penalized log-likelihood

The composite statistic is the profile likelihood ratio

where the second term lets the ensemble move within its weight uncertainty to best fit the data, and the first additionally switches on the localized kernels, so that only the discrepancy the weights cannot absorb contributes. As in the main text, the point-null statistic fixes the weights at their best fit,

d. Calibration. The training data (used to fit the flows and the weights $\widehat{u}$) are kept strictly separate from the sample $\mathcal{D}$ on which $t$ is evaluated. The null distribution of $t$ is built from $B$ pseudo-experiments. We draw posterior-predictive toys: each toy draws a weight vector $u^{(b)}\sim\pi$ and then a dataset $\mathcal{D}^{(b)}\sim g_{u^{(b)}}$, so that the null itself carries the weight uncertainty. Because the weight fit of Eq. (27) is convex, unlike the non-convex mixture fit of the main text, whose toys must be re-fit from a cold start, re-fitting each toy would recover a $u^{(b)}$ in the vicinity of $\widehat{u}$ regardless of initialization. We then expect drawing $u^{(b)}\sim\pi$ to perform similarly, so we do not re-fit the data-generating weights per toy and thus only show the Gaussian calibration. A full bootstrap calibration which refits the whole ensemble of flows in every pseudo-experiment could outperform this procedure, but would be computationally prohibitive.

Because the weights are signed, $g_u$ is not a probability mixture and cannot be sampled member-by-member. We instead draw toys by sampling-importance-resampling (SIR) from the positive proposal $q$. The statistic is then evaluated on each toy exactly as on the data, re-maximizing over $u$ and $\alpha$ as in Eq. (35), giving $\{t^{(b)}\}$, and the $p$-value of the observed statistic is the empirical tail probability $p=(\#\{t^{(b)}\ge t(\mathcal{D})\}+1)/(B+1)$. For the point-null statistic, toys are drawn from $g_{\widehat{u}}$ without sampling the weights. Because the tilt is regularized (ridge and box) and its kernel dictionary is strongly collinear, neither statistic follows the asymptotic Wilks $\chi^2$, whose number of degrees of freedom would equal the number of tilt coefficients; the effective number of degrees of freedom is much smaller. We therefore do not rely on Wilks’ theorem, but calibrate the statistics empirically from the pseudo-experiments, against which all $p$-values in Sec. B 4 are evaluated.

Results

We apply the construction of Secs. B 1–B 3 to the two-dimensional toy of § 2.2 A ($N_{\rm train}=10^5$), for ensembles of $M$ flows up to $M=128$. The tests use the kernel dictionary of the main text ($J=500$, $\sigma=0.4$), ridge $\lambda=1$ and box $A=3$. Each null is calibrated on $B=100$ pseudo-experiments, and the $p$-value is estimated from 100 draws of $N_{\rm test}$ events from the truth. We then repeat the construction on the four-dimensional jet embedding of § 2.2 B, with the same $N_{\rm train}$, ensemble sizes and test settings, and the four-dimensional kernel width of the main text, $\sigma=0.53$.

a. Density quality. Figure 10 shows the fitted two-dimensional marginals for $M=1,32,128$. A single flow ($M=1$), which carries no parametrized uncertainty, reproduces the bulk of the target but visibly mismodels the peaks and the valley of the bimodal feature 1. As $M$ increases the ensemble median tightens onto the target (most visibly the bimodal valley of feature 1) and a weight-uncertainty band develops, narrow in the bulk and widening in the low-statistics tails. In the four-dimensional embedding (Figure 11) the picture is different: the fitted marginals show structured deviations from the held-out data, most visibly around the sharp peaks, while the weight-uncertainty band remains negligibly small at every $M$. The propagated uncertainty therefore does not cover the discrepancy between the ensemble and the data, which anticipates the goodness-of-fit results below.

b. Goodness-of-fit. Figure 12 shows the median test $p$-value versus $N_{\rm test}$ for $M=32,128$. As expected, the test sensitivity increases with $N_{\mathrm{test}}$: the $p$-value stays near nominal while the test sample is too small to resolve the finite-$N_{\mathrm{train}}$ error of the fit, and falls once it does. The single flow degrades earliest, by $N_{\mathrm{test}}\sim10^4$, as anticipated in Sec. II: a flow’s large parameter count and imperfect inductive bias push the rejection onset to $N_{\mathrm{test}}\ll N_{\mathrm{train}}$. A single flow is therefore not a valid reference model at these test sizes. Ensembling and propagating the weight uncertainty push this onset outward, and the benefit grows with $M$. A moderate ensemble ($M=32$) delays the degradation to $N_{\mathrm{test}}\sim2$–$5\times10^4$ but still reaches the floor $1/(B+1)$ by $N_{\mathrm{test}}=N_{\mathrm{train}}=10^5$ with little difference between the composite and point-null (median $p\approx0.09$ vs $0.07$ at $N_{\mathrm{test}}=5\times10^4$). This happens although the $M=32$ ensemble closely reproduces the two-dimensional marginals (Figure 10): the test probes the joint density, where the residual error of the fit is resolved. Only at $M=128$ does propagating the weight uncertainty make a clear difference: the composite test stays near nominal, $p\approx0.43$, whereas the point-null falls to $p\approx0.16$ at $N_{\mathrm{test}}=5\times10^4$ and to the floor by $N_{\mathrm{test}}=N_{\mathrm{train}}$. With $M=128$ the composite test gives $p\approx0.43$ out to $N_{\mathrm{test}}=2\times10^5$, i.e. it does not reject the model even in the genuine oversampling regime $N_{\mathrm{test}}>N_{\mathrm{train}}$.

In four dimensions the picture changes (Figure 13). The single flow degrades by $N_{\mathrm{test}}\sim2\times10^3$ and the ensembles by $N_{\mathrm{test}}\sim10^4$. There is no improvement from propagating the weight uncertainty: the composite $p$-value is inseparable from the point-null within the pseudo-experiment variations at every $M$ and every $N_{\mathrm{test}}$, with no separation even at $M=128$. All constructions reach the floor $1/(B+1)$ by $N_{\mathrm{test}}\sim2\times10^4$, well below $N_{\mathrm{train}}=10^5$, and even $M=128$ is fully rejected.

c. Discussion. Since every null is calibrated with pseudo-experiments drawn from the composite-null itself, the loss of $p$-value at large $N_{\mathrm{test}}$ reflects a failure of the ensemble uncertainty to cover the error of the fit, not an artifact of the statistic. Ensembling the flows and propagating the weight uncertainty push this rejection outward in the 2D toy study, and with $M=128$ the model survives up to $N_{\mathrm{test}}\simeq N_{\mathrm{train}}$ and into the oversampling regime $N_{\mathrm{test}}>N_{\mathrm{train}}$ probed here, up to $N_{\mathrm{test}}=2N_{\mathrm{train}}$. In the more complex 4D problem, the construction fails to close: ensembling still improves, but the weight covariance $\Sigma_u$ spans only re-weightings $\sum_i\delta w_i f_i(x)$ of the frozen members, and in the embedding the resulting uncertainty is negligible compared to the discrepancy between the ensemble and the data. With an uncertainty this small, the profiled weights can barely move, the composite null effectively reduces to the point-null, and the two statistics coincide: even an ensemble of $M=128$ flows is rejected well before $N_{\mathrm{test}}=N_{\mathrm{train}}$. This is the limitation anticipated in Sec. B2: re-weighting frozen flows can make up for the error of a simple target, but not for that of a more complex one.

Implementation details

The nominal fit and every bootstrap refit use EM with a $k$-means initialization (specifically, the $k$-means++ initialization implemented in the GaussianMixture class in scikit-learn [78]). The scores of Eq. (6) can be computed analytically or by central finite differences,

which costs $2P$ evaluations of the log-density on each sample. The normalization of the tilted model for the alternative hypothesis, as well as the standardization of the features, are calculated using reference samples drawn from a pool of size $2\times10^6$ from the best-fit SIGMA model. $10N_{\mathrm{test}}$ samples are used to compute the normalization to ensure that its uncertainty is subdominant, and $10^6=10N_{\mathrm{train}}$ samples are used to compute the Fisher matrix. The feature standardization is performed with the whole pool.

The $J=500$ kernel centers $c_j$ are obtained by $k$-means clustering of a subsample of $R_0$. The width $\sigma$ of the kernels is taken to be uniform and isotropic, with $\sigma=0.4$ in the two-dimensional case study and $\sigma=0.53$ in the four-dimensional case study.

The kernel amplitudes carry the ridge penalty $\lambda_\alpha\lvert\alpha\rvert^2$ of Eq. (15) and are restricted to

which bounds the log-density distortion of each kernel by $A$. We use $A=10$, i.e., local distortions within a factor $e^{\pm10}$. Both regularizations are applied identically to the data and to every pseudo-experiment. $\lambda_\alpha$ is taken to be $1$ throughout.

The $\Lambda_{\max}$ hyperparameter is taken to be $0.008$ in all of the studies in the main body, and Figure 14 shows the effects of changing $\Lambda_{\max}$ in the 2D example with $K=64$. As anticipated in the main body, the performance of the Gaussian-calibrated null is strongly dependent on the choice of $\Lambda_{\max}$, since its sampling from the null hypothesis is directly dependent on $\Lambda_{\max}$, while the bootstrap-calibrated null is more robust, as $\Lambda_{\max}$ only affects the definition of the test statistic and not the sampling of the null.

Null distribution of the test statistic

We show that, to second order in $(\nu,\alpha)$ and for $\lambda_\alpha=0$, $\operatorname{tcomp}$ follows a $\chi^2_J$ distribution under the Gaussian composite null. Let $f=(\Phi,\tilde{k})\in\mathbb{R}^{M+J}$ denote the full feature vector, which has zero mean under $p_{\hat\theta}$ and, by Eq. (10), covariance

with $C_{\alpha\alpha}=\operatorname{Cov}[\tilde{k}]$ and $C_{\alpha v}=C_{v\alpha}^{\mathsf{T}}=\operatorname{Cov}[\tilde{k},\Phi]$. Write $N\equiv N_{\mathrm{test}}$, $a\equiv N_{\mathrm{test}}/N_{\mathrm{train}}$, $u=(v,\alpha)$, and $s=\sum_i f(x_i)=(s_v,s_\alpha)$ for the feature sums over the test sample. Expanding $\log\tilde{Z}(u)=\frac{1}{2}u^{\mathsf{T}}\Sigma u+O(\lVert u\rVert^3)$ gives

where the identity in the upper block comes from the Gaussian constraint. The numerator and denominator maximizations give $s^{\mathsf{T}}H^{-1}s/2$ and $\lVert s_v\rVert^2/2(1+a)$, respectively, and the block-inverse identity yields

Here $r$ is the vector of kernel sums with their projection on the score sums removed, and $S$ is the Schur complement of the $v$-block of $H$. Under the Gaussian composite null, the test sample is drawn from $p_{\hat{\theta}+Wz}\simeq g_z$ with $z\sim\mathcal{N}(0,I_M)$. The tilt shifts the mean of $f$ by $\Sigma_{\cdot v}z=(z/N_{\mathrm{train}},C_{\alpha v}z)$, so that, including the spread of $z$,

Substituting gives $\operatorname{Cov}(r)=S$. Since $r$ is asymptotically Gaussian, $t_{\mathrm{comp}}\sim\chi^2_J$, or $\chi^2$ with $\operatorname{rank}(S)$ degrees

of freedom if the kernels are linearly dependent on $R_0$. The result is independent of $\alpha$ and of the correlations $C_{\alpha v}$. The $M$ standardized directions drop out because the constraint on $v$ matches the distribution of $z$.

A ridge $\lambda_\alpha> 0$ adds $2\lambda_\alpha I_J$ to the $\alpha$-block of $H$, which leaves $r$ unchanged and replaces $S$ by $S + 2\lambda_\alpha I_J$ in the quadratic form. The statistic is then a weighted sum of $J$ independent $\chi^2_1$ variables with weights $\sigma_j/(\sigma_j + 2\lambda_\alpha) < 1$, where $\sigma_j$ are the eigenvalues of $S$. For $t_{\mathrm{point}}$, with $v \equiv0$ and data drawn from $p_{\hat{\theta}}$, the same argument gives $t_{\mathrm{point}} = s_\alpha^{\mathsf{T}}(NC_{\alpha\alpha})^{-1}s_\alpha\sim\chi^2_J$.

Figure Figure 15 examines the empirical null distributions for the four-dimensional particle-jets problem of § 2.2 B with $K = 256$ and $N_{\mathrm{train}} = 10^5$, for the three null prescriptions of Sec. IV D and increasing test-set size. At $N_{\mathrm{test}} = 1\mathrm{k}$, all three empirical nulls agree with each other and with the observed statistic, and a $\chi^2$ fitted to each null returns $\hat{\nu} \approx93$, consistent with the predicted effective $\nu\approx82$ and the presence of additional noise due to the finite sample used to estimate $Z$. In the small-sample regime, the asymptotic approximation holds and the test cannot resolve the finite-sample uncertainty, and the three calibrations are interchangeable.

The tests separate as $N_{\mathrm{test}}$ grows. At $N_{\mathrm{test}} = 20\mathrm{k}$, the two composite nulls still nearly reproduce the observed distribution (with bootstrap somewhat outperforming Gaussian), while the point-null distribution has already shifted below that of the observed test statistic.

At $N_{\mathrm{test}} = 10^5 = N_{\mathrm{train}}$, all three null distributions fall below the observed test statistic distribution (consistent with the $p$-values in Figure 9), but the point-null and Gaussian prescription fail dramatically while the bootstrap calibration remains fairly close to the observed distribution.

Effects of kernel width

In this appendix, we study the effects of the choice of kernel width $\sigma$ on the behavior of the goodness-of-fit test. As discussed in the main body, we find that choosing $\sigma$ too small relative to the length scale of the data results in spurious calibration failure, establishing a tradeoff where $\sigma$ must be chosen large enough for the test to be reliable, and small enough for it to be sensitive.

In Figure 16, we show results of the closure test for the two-dimensional example for $K = 16, 32$ with all parameters chosen as in the main body, except $\sigma= 0.1$ instead of $0.4$. As with the larger kernels, the point null is insufficient and is rejected for all but the smallest values of $N_{\mathrm{test}}$; the Gaussian null outperforms the point null, but is also miscalibrated for $N_{\mathrm{test}} \ge10\mathrm{k}$.

On the other hand, the bootstrap-calibrated null deviates from closure in the opposite way: the median $p$-value is greater than $0.5$ at large $N_{\mathrm{test}}$, indicating that the distribution of the test statistic is shifted rightward in calibration pseudo-experiments relative to the distribution in observed data. In other words, the test finds more discrepancy, on average, with the best-fit trained model in pseudo-experiments (where that model is guaranteed to be well-specified) than in the observed data (where it is not), so this deviation from calibration is clearly spurious.

We ascribe this to an interplay between the structure of the SIGMA model, which is a GMM, and the test. When trained on a finite training set without strong regularization, GMMs generally overfit by placing narrow components in regions of feature space with disproportionately many events. This overfitting is not pathological in our context, as it is part of the variance we aim to capture with uncertainties (as opposed to model bias, which we cannot capture; strong regularization would induce such bias).

The fine-grained structure induced by overfitting is then generically present in our trained models but not in

the data-generating process, so an alternative hypothesis comprised of small-length-scale Gaussian functions will systematically fit a dataset drawn from a trained SIGMA model better than it would observed data. As such, we do not recommend using SIGMA models in conjunction with a goodness-of-fit test sensitive to length scales substantially smaller than the natural length scales of the data, hence our choice of $\sigma$ in the main body to be comparable to the data standard deviation.

References