Unifying supervised learning and VAEs - coverage, systematics and goodness-of-fit in normalizing-flow based neural network models for astro-particle reconstructions

Neural-network based predictions of event properties in astro-particle physics are getting more and more common. However, in many cases the result is just utilized as a point prediction. Statistical uncertainties and coverage (1), systematic uncertainties (2) or a goodness-of-fit measure (3) are often not calculated. Here we describe a certain choice of training and network architecture that allows to incorporate all these properties into a single network model. We show that a KL-divergence objective of the joint distribution of data and labels allows to unify supervised learning and variational autoen-coders (VAEs) under one umbrella of stochastic variational inference. The unification motivates an extended supervised learning scheme which allows to calculate a goodness-of-fit p-value for the neural network model. Conditional normalizing flows amortized with a neural network are crucial in this construction. We discuss how they allow to rigorously define coverage for posteriors defined jointly on a product space, e.g. R n × S m , which encompasses posteriors over directions. Finally, systematic uncertainties are naturally included in the variational viewpoint. The proposed extended supervised training with amortized normalizing flows incorporates (1) coverage calculation, (2) systematics and (3) a goodness-of-fit measure in a single machine-learning model. There are no constraints on the shape of the involved distributions (e.g. Gaussianity) for these properties to hold, in fact it works with complex multi-modal distributions defined on product spaces like R n ×S m . We see great potential for exploiting this per-event information in event selections or for fast astronomical alerts which require uncertainty guarantees.


Introduction
Deep neural networks have achieved great results over the last couple of years.While the breakthroughs were initially made in industrial applications, for example in image processing [1], in recent years their application in fundamental science has become more and more ubiquitous.A typical application of deep learning in experimental high-energy physics is concerned with the reconstruction of particle interactions [2].These include discrete quantities like particle type [3] or continuous quantities like direction, position and energy [4].Traditionally, reconstructions of continuous quantities are performed by parametrized likelihood fits [5] which allow to calculate confidence intervals with standard Frequentist or Bayesian methods [6].Neural networks on the contrary are often used to produce point estimates [7], and there is no universal agreed-upon notion how to calculate confidence intervals.Often they are just ignored, the result is registered 1 arXiv:2008.05825v5[cs.LG] 14 Jan 2024 as an observable and used in a binned likelihood analysis [4].Per-event uncertainties are not necessarily required for this use-case.However, there are many situations where precise uncertainty quantification is important.An example are per-event reconstructions in high-energy neutrino telescopes like IceCube [8] or gravitationalwave observatories like LIGO [9] in the context of multi-messenger astronomy [10].These experiments send out alerts to the astronomical community to perform follow-up observations with other experiments.On the one hand, these alerts are time-critical and should be sent out to the public as fast as possible.On the other hand, the uncertainties of the inferred direction must be precise and in particular not biased.While the time-critical aspect of this use-case is in favor of neural networks, obtaining unbiased uncertainties can be a challenge.A naive solution could employ Bayesian neural networks [11] or approximations like certain dropout-variations [12] or ensemble-methods [13].However, these methods model a posterior over network weights, not over the actual physics parameters.An alternative is to parametrize the likelihood function with a neural network and perform a standard likelihoodbased Frequentist or Bayesian analysis.Recently, this was done for gravitational-wave signals with flow-based networks [14].Such methods, however, inherit the disadvantages of likelihood-based inference.They can be fall into local optima during an optimization routine or have potentially long running times when Markov Chain Monte Carlo [15] (MCMC) is used to obtain samples.An alternative that recently has gotten popular is likelihood-free inference with neural networks, in which the posterior is directly learned from data.This has been applied to gravitational wave posteriors modeling the posterior as a parametrized Gaussian or mixture of Gaussians [16] and going beyond to use autoregressive normalizing flows for more complex shapes [17].This paper is about the same line of thinking, where the posterior or more generally parts of the joint distribution of data and labels are learnt using normalizing flows.
We first discuss that the training process in supervised learning can be recast as variational inference of the true posterior distribution over labels, where the variational approximation is parametrized by a neural network.The derivation involves the Kullback-Leiber (KL)-divergence [18] of the joint distribution of data and labels.This is a known derivation and the final loss function is sometimes called "conditional maximumlikelihood objective" in the machine learning literature [19].While the loss indeed represents a likelihood with respect to the neural network parameters, we emphasize it is more useful to think of it as an approximation of the posterior over the latent variables.Standard supervised learning usually uses the mean-squared-error (MSE) loss function, which corresponds to a standard-normal posterior as an approximation of the true posterior.This can often be a bad approximation.We compare it to the slightly more complex case of a Gaussian with a single covariance parameter and in particular approximations based on normalizing flows [20], which in principle have arbitrary approximation capabilities.Since all Gaussian approximations, including the basic MSE loss, can be thought of as special cases of normalizing flows, so-called affine flows, it seems sensible to view all supervised loss functions from this angle.While normalizing-flow base distributions can be arbitrary, it turns out that there are advantages of starting with a standard normal as a base distribution that go beyond the argument of simple evaluation.The standard normal automatically allows for straight-forward coverage tests, which are discussed in Section 6.
Usually, variational inference with neural networks is employed in unsupervised learning in the context of variational autoencoders (VAEs).As we will show (Section 2.2), one can derive the Evidence Lower Bound (ELBO) [21] of a variational autoencoder from the same joint KLdivergence as the supervised "maximum likelihood objective".This derivation not only explicitly indicates how supervised learning and VAE training is connected, but it also sheds some light on the interpretation of VAEs.Importantly, we use this connection to motivate a mixed training scheme which we call "extended supervised" learning (Section 2.3).This allows to calculate Bayesian goodness-of-fit values (Section 8) for the trained model.Figure 1 indicates an overview of the resulting picture.Traditional variational inference (Fig. 1 a) is typically discussed on a per-event level or in the context of neural networks as approximating the posterior over the network parameters in Bayesian neural networks [11].However, one is really interested in the posterior over physics parameters, ideally for all events simultaneously.Figure 1 b) illustrates that supervised learning, extended supervised learning and unsupervised VAEs can be interpreted as stochastic variational inference using the KL-divergence of the joint distribution of data, observed and unobserved labels (latent variables).Because only sub-parts of the joint distribution are effectively parametrized in these approaches, one obtains "explicit" variational inference solutions for the posterior distributions, in the sense that the conditional structure of the posterior is explicitly modeled and thereby different for each datapoint.This is sometimes also called "local" variational inference [21] in the literature.The unified viewpoint in combination with amortized conditional normalizing flows naturally leads to answers to the following questions: 1. How can we do coverage tests with neural networks on complex base distributions, including distributions over directions? 2. How are systematic uncertainties incorporated in the training process? 3. How can we do goodness-of-fit checks on neural-network predictions?
The first half of the paper is concerned with the unified viewpoint and derives the various loss functions from the joint KL-divergence.The second half then answers the three questions above.

Monte Carlo estimates and the joint KL-divergence
Physics experiments perform measurements on the final outcomes of a causal chain of events.These measurements are inherently probabilistic due to noise and the randomness from particlephysics interactions, and the measured observable follows a specific probability distribution.The probability distribution over possible measurement outcomes is called the likelihood function when viewed as a function of its parameters.It is a central object in both Frequentist and Bayesian statistical analyses to perform parameter estimation [6].Its shape is entirely determined by the laws of nature in combination with the detector response.However, because the laws can be convoluted and the experiment can be very complex, it is usually not possible to write down an explicit analytic expression.A common practice is to estimate the likelihood function from Monte Carlo simulations where all these complex effects are considered [22].This estimated likelihood function is then used to perform inference or calculate confidence intervals.
Neural networks allow to skip this estimation step completely, because the Monte Carlo samples themselves are drawn from the true joint distribution P t (x, z o ) of observations x and parameters of interest z o .Let us call the corresponding true data generating function P t (x; z o ) which is the true probability distribution of the measured data x given the properties z o .Here z o stands for recorded or observed properties in the simulation, for example the position of a particle interaction.Connected to this data generating function, there exists a true posterior distribution P t (z o ; x) and a true prior distribution P t (z o ).The true joint distribution follows the distribution which includes the detector response and selection effects inherent to the measurement -it also includes artificial selection effects.For example, if the generating function of the direction of injected particles is uniform, the actually recorded Monte Carlo events will generally not be uniform due to detector effects.The implicitly contained true prior P t (z o ) is this non-uniform distribution over z o of the actually registered events, not the uniform generating distribution.
Since a Monte Carlo simulation draws samples x i , z o,i from the joint distribution P t (x, z o ), we can always evaluate any expectation value under the true joint probability distribution as where i indexes the N samples.The following sections make use of a specific choice for f , namely f = ln Pt(x,zo) q(x,zo) , which yields an expectation value that equals the KL-divergence between two distributions P t (x, z o ) and q(x, z o ).The KL-divergence [18] is a natural quantity to measure the distance Fig. 1: Variational inference (VI) examples for simulation data x i and labels z i indexed by i = 1 . . .N , which comprise the whole dataset of size N .The exception is the single event example in (a), which has a single datum x j as input.
between two probability distributions.In particular, if q ϕ is a parametrized probability distribution, the KL-divergence defines a loss function over ϕ that achieves its minimum when q ϕ is equal to P t .It is therefore often used in variational methods which perform inference via optimization [21].It turns out that the joint KL-divergence can be used to derive the loss functions in both supervised learning and unsupervised VAEs and thereby unifies them as two connected approaches to variational inference in slightly different circumstances.

Supervised learning
The KL-divergence of the true joint distribution P t (z o , x) and an approximation q ϕ (z o , x) can be written as

=
x zo

=
x zo The distributions involving P t can not be evaluated analytically, but as discussed above Monte Carlo simulations yield samples from P t and hence provide sample-based estimates of the integrals.We only parametrize the conditional distribution over labels, q ϕ (z o ; x) with ϕ, and leave the distribution q(x) unparametrized as it is typically not of interest.Taking the sample estimate yields the following update rule over ϕ to minimize the KL-divergence objective and thereby minimize the supervised loss function L supervised (ϕ): = arg min The approximation of the marginal likelihood, q(x i ), can be dropped as a constant part with respect to changes in ϕ.Additionally we can drop the true posterior evaluations P t (z o,i ; x i ) which are also constant with respect to ϕ.If the distribution q ϕ is actually parametrized by a network whose parameters are ϕ, the derivation shows that minimizing the KL-divergence between the true posterior and an approximation given by a neural network is equivalent to standard neural network training where the goal is to minimize negative log-probability over labels.This has been discussed previously [19], but it is usually described as a maximum likelihood objective with respect to the network parameters.Here, we emphasize that it is really more useful to think of q ϕ (z o ; x i ) as a parametrized posterior over labels z o given the data x i , not as a conditional likelihood function with respect to the parameters ϕ.In unsupervised learning the labels z o,i are not available so the precise above objective does not work.In the following we show that a simple replacement of the posterior KL-divergence in eq.7 with the respective reverse KL-divergence leads to a tractable solution via the reparametrization trick [23] which contains the evidence lower bound (ELBO) and thereby the VAE objective.

Unsupervised learning: Variational autoencoders
The following derivation is motivated by an experimental physicists' point of view when there is access to a Monte Carlo simulation.Here, one way to think of latent variables in unsupervised variational autoencoders is to imagine them as unobserved labels that are not recorded in the simulation.The implied direct comparison, and renaming, of "observed" into "unobserved" labels is justified if the latent variables have positive mutual information with the data.Mutual information is a non-linear generalization of correlation [24].In events of a static particle physics detector with a single particle interaction, for example, the only properties that can have positive mutual information with the data are properties of the particle interaction.These properties can be in principle be labeled.Typical labels are the position or direction of a particle at the interaction point. 1he aim of the following exercise is to start with the same KL-divergence as for the supervised loss, but with unobserved labels z u instead of observed labels z o , and deduce which modifications have to be made in order to obtain the ELBO loss of a variational autoencoder as a partial tractable objective.It is not the aim to deduce strict new numerical results, but to indicate how supervised learning and unsupervised learning are precisely connected in the variational viewpoint.
Starting with the same joint KL-divergence, but now using unobserved labels z u , we can write = E x [D KL (P t (z u ; x); q ϕ (z u ; x)] which results in two intractable terms.The second term D KL,M,const can for further discussion be ignored since q(x) is typically not parametrized and therefore irrelevant in optimization schemes.
In the first term, the outer integral over x is explicitly kept, because we do have samples from x and so we can in principle evaluate expectation values over x.The inner part, however, involves an expectation over the intractable KL-divergence of the conditional distribution P t (z u ; x) with q ϕ (z u ; x) (orange).In order to proceed to some extent, we could replace the KL-divergence within the expectation value with any generalized divergence measure D(P t ; q), for example f-divergences [25], as long as it shares the property that it has a minimum when P t (z u ; x) = q ϕ (z u ; x).However, there is only a subset that will straightforwardly lead to the ELBO loss of the VAE.Among those, a natural choice is the reverse KL-divergence, D rev.KL (q ϕ (z u ; x); P t (z u ; x)), which we choose for simplicity.This leads to the following surrogate loss term: As long as q ϕ has arbitrary approximation capabilities, the surrogate term has the property that the minimum over the parameters ϕ is equal to the minimum of the original KL-divergence arg min ϕ D KL,joint(x,zo) (P t ; q ϕ ) = arg min ϕ Surrogate KL (P t ; q ϕ ), (15) because both KL-divergence (orange) and reverse KL-divergence (blue) are equal to zero when the two involved distributions are equal.So far this is just a theoretical exercise, as P t is inaccessible.The surrogate loss, however, contains the ELBO loss.In order to see this, we can rewrite the surrogate loss as which results in the negative ELBO loss, two constant terms, and a residual term R(ϕ, θ, ψ).The important step in this reformulation is the introduction of an auxiliary distribution p θ/ψ (z u , x) = p θ (x; z u ) • p ψ (z u ) from eq. 18 to eq. 19.The negative entropy term in eq.18 by itself is tractable, but behaves divergent when minimized over ϕ, since the negative entropy will tend to infinity.Therefore, the introduction of p θ/ψ (z u , x) is a necessary ingredient to obtain anything that has a chance for non-trivial behavior.We also indicate by the subscript ψ that the prior p ψ (z u ) can be parametrized by ψ, even though in many applications it is just taken to be a fixed standard normal distribution.With the introduction of p θ/ψ (x, z u ), the surrogate KL-divergence then splits into two terms in eq.19.
The first term is the KL-divergence between P t (x) • q(z u ; x) and p θ/ψ (x, z u ), which is equal to the ELBO loss when the constant term E x [lnP t (x)] is pulled out of the integral.This joint KL-divergence between P t (x) • q ϕ (z u ; x) and p θ/ψ (x, z u ) has been used before in the context of the InfoVAE [26] or more general VAE architectures with additional constraints [27] and is by itself another non-standard starting point to derive the ELBO loss.Here it arises as a product in the derivation which is connected to the supervised loss derivation via the KL-divergence of the joint distribution.It is important to note that we have to use this more complicated construction, compared to just start with this simpler KL divergence, in order to see the similarity to the supervised loss derivation and then be equipped with a canonical way to derive the extended supervised case in Section 2.3.
The second term is a residual term R(ϕ, θ, ψ).This term is not tractable, because any samples drawn from q ϕ (z u ; x) can not be evaluated by the inaccessible density P t (z u ; x).However, one can deduce that after an ELBO optimization with solution θ * and ϕ * and flexible enough density parametrizations, the residual term is always bounded from below by zero because it is equal to a proper KL-divergence: This construction therefore makes it explicit that ELBO optimization alone does not necessarily lead to a joint density that matches the true density because R can be greater than zero -a fact that is lost in the standard ELBO derivation [28] based on the marginal likelihood and Jensen's inequality and therefore also in many papers since the original VAE paper [23].Furthermore, beside the assumption that we could in principle observe the latent variables but choose not to, hence the term "unobserved", in general we do not know the exact values.Because of the symmetry of the KL-divergence under diffeomorphisms, this allows for extra functional freedom of the distributions and the involved mapping.Instead of complete freedom, however, it was pointed out in [29] that the determined final distribution matches the true one up to certain class of transformations A, which they call "identifiability" up to A. In particular, when all terms in the ELBO depend on extra observed input, the prior has a certain structure and the data x are Gaussian observations, A turns out to be a global scaling and permutation of latent dimensions. 2 Regarding the above derivation we can write R(ϕ * , θ * , ψ * ) ∼ A 0 to denote this situation, which means the residual term is zero within the identifiability class A. Extra constraints on the mutual information between data and labels [26], the total correlation of latent variables [30], a better prior parametrization [31] or extra conditional dependencies [29] are often used to improve the ELBO solution.From the above derivation these are all well motivated, as all of those use extra constraints besides the ELBO and therefore have the potential to reduce R(ϕ * , θ * , ψ * ) within a given identifiability class A whose properties will depend on the constraints.There are also approaches that change the relative strength of the data PDF term with either a data-PDF or latent-PDF prefactor [32] [26].These can be motivated with a balancing of the often very different dimensionality between the two PDFs.
For the discussion in the following, we will set this relative scaling to unity without loss of generality.We can further form the sample approximation of the ELBO and add the above mentioned constraints to form a loss function arg min θ,ϕ,ψ = arg min where the term C(θ, ϕ, ψ) indicates the constraint.For a total correlation constraint [30], for example, we would have C(θ, ϕ, ψ) = γ • D KL (q ϕ (z u ); j q j,ϕ (z u )), where the index j describes the different marginal distributions of q ϕ .There is also a tunable prefactor γ to balance the different loss terms.As discussed before, if flexible enough density estimators are used and the final parameter solution is denoted by θ * , ϕ * and ψ * , it follows that R(ϕ * , θ * , ψ * ) ∼ A 0, the ELBO saturates, and Surrogate KL (P t ; q ϕ * ) ∼ A 0 and D KL,joint(x,zo) (P t ; q ϕ * ) ∼ A 0.
In our view, a few advantages arise from this derivation of the variational autoencoder.
1.The derivation is connected to the KLdivergence derivation of supervised learning.2. It explicitly shows that there are three joint distributions involved.The first is P t (z u ; x) • P t (x), the true underlying joint distribution.
The second is q ϕ (z u ; x) • P t (x), a distribution where the conditional distribution over the latent variables is exchanged for a tractable approximation q ϕ (z u ; x).The third is another tractable approximation p θ (x; z u )•p ψ (z u ) which is parametrized in the opposite causal direction.In the literature, the true distribution P t is typically simply denoted by p, going back to the original ELBO marginal likelihood derivation [28] or VAE publication [23].This can be confusing, as p (i.e.P t ) and p θ are often used interchangeably.3. Extra constraints [27] [29] where the first term involves an intractable KLdivergence similar to the VAE and the second term is the supervised loss function which we abbreviate by S. Next we construct a surrogate term similar to the VAE derivation in order to obtain a tractable objective.
Surrogate KL,1 In this term we again replace the KL-divergence (orange) by a reverse KL-divergence (blue).In contrast to the VAE, we can define a second surrogate term as where the expectation is taken with respect to the parametrized distribution q φ (z o ; x) (turquoise) that is determined in the supervised part S. The tilde indicates that gradients with respect to φ are not propagated through, which is used to completely decouple the supervised from the unsupervised part during learning.This second surrogate loss will be useful for consistent goodness-of-fit procedures (Section 8 = arg min since the supervised part leads to q φ (z o ; x) ≈ P t (z o ; x) and the surrogate losses have the same global minimum as the joint KL-divergence, similar to the VAE derivation.
Next, we reformulate the second surrogate loss as and we end up with a supervised term S and similar terms to the unsupervised derivation.The residual term R(θ, ϕ, φ, ψ) again is equal to a KLdivergence up to a certain identfiability class after training.There is a slight difference, in that the prior in the ELBO is now defined over the joint space (z u , z o ) instead of just z u alone.The same derivation would work with the first surrogate loss by replacing q φ (z o ; x) with P t (z o ; x) everywhere.In particular, if in addition to using the first surrogate loss, the joint prior is split up , every part of the resulting ELBO objective is conditioned on z o as an extra parameter.This is an important factor for identifiability guarantees, as outlined in [29], and might be interesting to study on its own.
We use the second surrogate loss here to later have self-consistent goodness-of-fit calculations (see Section 8) at all times during training.In the following, we discuss two types of loss functions that can be formed using the previous results.

Extended supervised loss
The first is a loss definition that can be used to perform what we call extended supervised training, and will be important for the calculation of a goodness-of-fit as described in Section 8.It is a sample-based application of the tractable parts in eq.34: the ELBO-like term and the supervised part S(φ): = arg min We again add a constraint term C(θ, ϕ, φ, ψ) similar to the VAE loss to potentially improve the identifiability of the unsupervised dimensions z u .The true observed labels are denoted by z o,i , samples from q φ (z o ; x) are denoted by z o,i,φ .The symbol ∼ indicates that the gradient is not propagated in order to decouple the supervised part during training.If the other surrogate term Surrogate KL,1 (eq.27) had been used, this decoupling would have happened automatically.Because the supervised training is effectively decoupled, one can also choose to first train the supervised part L supervised , and only later train the rest.An extended supervised training can therefore always be started with an already finished supervised model and can be viewed as an add-on to it.

Semi-supervised learning
The extended supervised loss can also be adapted for semi-supervised learning.In semi-supervised learning, parts of the training data have labels, and parts are unlabeled.In the derivation of the VAE loss we argued that latent variables that share mutual information with the data can in principle be labeled.In semi-supervised learning, this assumption is automatically implied -parts of the data are not labeled, but could be in principle, if the data comes from a Monte Carlo simulation.Taking the variational viewpoint of the previous sections, a naive solution might be to use the supervised loss for data with labels, the VAE-loss for data without labels, and use the same distribution q ϕ in both losses to share parameters.A more natural solution, however, is to use the extended supervised loss instead of the pure supervised loss.
= arg min For unlabeled data, we take the loss of the variational autoencoder over the combined space z comb = {z o , z u } and treat the combined variable as unsupervised.The parametrization does not change for data with or without labels, only the sampling strategy differs.Therefore, the Ansatz seems to be an elegant and natural way to unify both types of data.

Toy Monte Carlo
Several toy Monte Carlo datasets have been created to perform empirical tests in the following sections.They are designed to mimic electromagnetic showers of electron-neutrino interactions [33] in a Cherenkov neutrino detector like IceCube [8] in a 2-D setting.In reality, such showers consist of charged particles that emit Cherenkov light in a coherent light front at the Cherenkov angle that changes its shape as detection modules are further away from the interaction point [34].For ice, which is the detection medium in IceCube, the light front becomes nearly isotropic for large distances.In general, depending on the position and orientation of the shower with respect to a detection module, the shape of the resulting arrival time distribution of photons varies.In this toy simulation the arrival time distribution is parametrized by a gamma distribution.In addition, the expected number of detected photons falls off exponentially with distance and depends on the orientation of the shower with respect to the module.These PDFs are parametrized such that they qualitatively mimic the photon PDF behavior of the real IceCube detector [8].At detection, the actual number of observed photons follows a Poisson distribution.Fig. 2 illustrates the first two datasets that are used in the next sections.
Fig. 2 (a) shows a detector that consists of a single module with two example neutrinos A and B. The two bars next to the photon collection module indicate the logarithmic expected mean of the number of observed photons from each neutrino interaction.Fig. 2 (b) shows the two corresponding photon arrival distributions.For event A, the data distribution is more concentrated than for event B, because the particle interacts closer to the detection unit.Also the number of observed photons is larger.Fig. 2 (c) shows the configuration of a second set of Monte Carlo simulations for a larger detector and a simulated threshold condition that at least 5 photons are observed for an event to be recorded.Also depicted are two associated example events (C and D) together with their expected number of photons in the various photon collection modules.
A summary of all datasets is given in table 1. Dataset 1 and 2 have two intrinsic degrees of freedom, the position of each neutrino interaction.The deposited energy and thereby emitted photons is always the same.Dataset three has an additional degree of freedom by also randomizing the direction.Dataset four also simulates a falling energy spectrum between 1 and 100 GeV.Datasets five to seven involve a larger detector and also four degrees of freedom.The last dataset, which contains track-like topologies, effectively emulates relativistic muon tracks by putting multiple eletromagnetic showers along a track that moves with the speed of light and distributing the energy in equal parts among those losses.More information on how these more complex datasets are used is given in Section 8.
The detection process can be described by an inhomogeneous tempo-spatial Poisson point process [35].The spatial part is restricted to the position of the collection modules, while the temporal detection can happen at all times.The corresponding likelihood function for observed labels z o is an extended likelihood [36], L(z o ), which we can write as an explicit evaluation of the data generating distribution D(k, t; z o , x).For a total number of modules N , this can be written as where λ j is the expectation value of a Poisson distribution for module j, k j the detected number of photons in module j, and p j (t i ) the probability distribution of photon arrival time t i in module j.The Poisson mean λ j and the shape of p j (t i ) depend on the event parameters z o and module positions x j .

The importance of flow-based models
All derivations so far assumed non-specific parametrizations of conditional probability density functions, i.e. posteriors and data generating PDFs, with neural networks.A general way to parametrize a probability density function with a neural network is via conditional amortized normalizing flows [37] [20].Normalizing flows are defined using a flow-defining function ρ ⃗ F (ẑ) whose parameters we define as ⃗ F .This function can be used as z o = ρ ⃗ F (ẑ) to transform a base random variable ẑ, usually following a standard normal distribution, to the desired target random variable z o .For normalizing flows to work, this function has to be invertible and differentiable, i.e. it has to be a diffeomorphic.If these properties are satisfied, then denoting the probability density of the base by p b (ẑ) and the target by q(z o ) we can calculate the log-probability of the target where is the Jacobian of ρ ⃗ F with respect to ẑ. Flows can be composed of multiple other flows and the resulting log-probability simply involves a sum over all log-determinants.In general p b (ẑ) can be arbitrary, but for simplicity and the possibility to calculate coverage (see Section 6) we use the   2 or of 400 modules as depicted in Fig. 11 standard normal distribution.One can also generate samples from q(z o ) by first sampling ẑi from p b (ẑ) and then transform the samples via where the samples now depend on the flow parameters ⃗ F .This makes the samples differentiable, known as the reparametrization trick [23], and is a key feature in the variational autoencoder and extended supervised losses (Section 2.3).

Conditional normalizing flows
Standard normalizing flows only describe PDFs without conditional dependencies, but we would like them to describe conditional PDFs like the posterior distribution.There are in principle multiple ways how this can be achieved.One way is to extend specific normalizing flows designed for high-dimensional image data, like NICE [38] and GLOW [39].These type of flows contain neural-network conditioners, typically MLPs, as part of their flow definition, and one can add a data representation as an additional input to these conditioners.This has been done before, e.g. in [40].There are another class of normalizing flows which are parameter-efficient Euclidean normalizing flows without coupling layers like "radial flows" [37] or "Gaussianization flows" [41] that have been shown in low Euclidean dimension (D ⪅ 20) to have a good performance on density estimation.This is the dimensionality regime we are interested in.These types of flows also allow a natural way to create a conditional normalizing flow by just predicting all the flow parameters by a neural network, which is the strategy we follow in this paper.The way to describe such a conditional normalizing flow q(z o ; x) with conditional input x is to make the transformation ρ ⃗ F dependent on x via a non-linear neural network transformation, as indicated in Fig. 3 (a).In general, a neural network with parameters ϕ that takes x as an input predicts the parameters ⃗ F , which in turn defines ρ ϕ (ẑ; x) ≡ ρ ⃗ F ϕ (x) (ẑ).The log-probability of a conditional normalizing flow PDF then looks like Compared to a standard normalizing flow (eq.40) which has flow parameters ⃗ F , the free parameters of such a conditional normalizing flow are actually the neural network parameters ϕ of the network used to encode the data x.As a specific example, Fig. 3 (b) shows an affine flow where a neural network predicts a mean vector μ and a width σ when the base distribution is a standard normal distribution.The resulting probability distribution q ϕ (z o ; x) of this flow describes a symmetric Gaussian distribution with mean μ and standard deviation σ.This is a common choice in some regression problems in high-energy neutrino physics [42] and used later in some comparisons.
The most common choice in supervised learning is to fix σ = 1, which results in a Gaussian PDF with unit variance which corresponds to the Mean-Squared-Error (MSE) loss function.Generic conditional normalizing flows with a standard normal distribution as base distribution therefore naturally generalize the MSE-loss.

Normalizing flows on spheres and tensor products of manifolds
Normalizing flows can be defined on manifolds like 1-spheres (S 1 ) or 2-spheres (S 2 ).In physics, 2-spheres are in particular interesting because directions in space are naturally defined on the 2-sphere, in particular if one wants to avoid issues in the polar regions which come up when looking at the zenith and azimuth seperately.Manifolds of dimension n are always embedded in an Euclidean embedding space R n+1 .In general, as described in [43], the log-determinant factors from eq. 40 and eq.41 of manifold flows differs from Euclidean flows as ln det(J Eucl. ) where the Jacobian J Emb. is calculated as if the transformation acts in embedding coordinates and an additional orthonormal projection matrix E projects into the tangent space at the transformation coordinates.If the Jacobian corresponds to an Euclidean transformation, the formula reduces to the standard case.In the end, manifold normalizing flows work similar to Euclidean normalizing flows, in that one can calculate the log-probability and sample efficiently from the distribution.The base distribution p 0 also has to be defined on the manifold, and for spheres a typical choice is the uniform distribution [43].In section 6, we argue that it is actually useful to have a fixed transformation from the Gaussian distribution in the Euclidean plane to the uniform distribution on the sphere as an additional step to facilitate coverage calculation.
In order to describe distributions that are defined simultaneously over the direction and position of an interaction, we can define a normalizing flow over a product space of manifolds.In the example of a 2-dimensional position and 1-dimensional direction, which appears in the toy simulation in Section 3, this space would be R 2 × S 1 , i.e. 2-dimensional Euclidean space and a 1-sphere.We use an autoregressive structure similar to [44] to connect the PDFs over the different manifolds and create a joint PDF on the tensor product space.In order to still use a single combined Gaussian base distribution, we again employ some transformations from the plane to the sphere for all manifold sub-parts as described in Section 6.This allows to calculate coverage also for such tensor product distributions.

Architecture and training
In all further comparisons in this paper, we split the parameters of flow-based posteriors into two parts.The first part consists of an encoder with gated recurrent units (GRUs) [45] and a subsequent aggregate MLP to encode the data into an internal representation h.The GRU reads in a time-ordered sequence of photons, where each photon is characterized by its detected position and time.The second part is a multilayer perceptron with 2 layers that further maps that internal representation to the respective flow parameters F .The process is illustrated in Fig. 3 (c).The generative model used in section 8 uses label and latent variables as an input to an MLP that maps to the flow parameters.
For the training of conditional normalizing flows in supervised learning, for example in section 5, we evaluate the log-probability of the labels for a given input 41 and minimize the negative log probability as defined in eq. 9 in batches using stochastic gradient descent.
For extended supervised training, which is used for the goodness-of-fit calculation in section 8, we additionally learn a generative model for the data and an additional posterior over a latent space as defined in the loss function 36.The additional latent space posterior and generative model are trained using the extra ELBO term.
We train the supervised loss and the ELBO at the same time in batches using stochastic gradient descent, but make sure that gradients are strictly separated, which corresponds to the "second surrogate loss" described in section 2.3.
At the end of training we adopt stochastic weight averaging (SWA) [46] in all cases.We found this to reduce the fluctuations and find a more stable solution.More details on the architecture and on the training procedure is given in appendix A.

Computational efficiency
One of the advantages of conditional normalizing flows compared to classical likelihood approaches is its computational efficiency.Astronomical alerts sent by the IceCube detector, for example, often take hours from the time the neutrino enters the detector until the alert is sent to the community as an astronomical telegram (ATEL).The reason are time-intensive profile likelihood scans that are needed for uncertainty contours.Spherical conditional normalizing flows, on the other hand, can produce a full-sky scan in seconds using multi-resolution HEALPIX [47] grids because they can both sample and evaluate the PDF.In a first step, samples drawn from the normalizing flow define the grid by guiding where pixelation has to be finer, and where it can be coarser.In a second step the PDF is then evaluated with the multi-resolution grid found in step one.With the other properties described in this paper, like systematics inclusion and coverage guarantees, this makes it an appealing alternative to likelihood scans that take hours.The only time consuming aspect is the training, which can take a few weeks on a single GPU for more complex models.This is of course not really an issue, as the model only needs to be trained once and can then readily be used for inference.

An example application with conditional Euclidean flows
In the following we study three Euclidean conditional normalizing flows in a simple example to infer the position of a neutrino interaction in the toy Monte Carlo.With the true likelihood of the true data generating function (eq.39) and N (ẑ; 0, 1)   the inherent prior distribution we can construct a true posterior distribution.We assume a flat prior for simplicity.We can then compare the true posterior with the posterior approximations obtained using various normalizing flows after training.Fig. 4 shows such a comparison for the example events from section 3. It compares a flexible Gaussianization flow [41] with an affine   flow (see Fig. 3 for details).For dataset 1, which consists of a single photon collection module, the resulting posterior is highly non-Gaussian.The flexible Gaussianization flow can approximate the posterior much better than an affine flow with a single width parameter.For dataset 2 the events generally contain more observed photons and the posteriors become more Gaussian-like as dictated by the Bernstein-von-Mises theorem [48].Both flows have more comparable shapes in this example.To assess performance more generally, Fig. 5 shows the posterior approximation quality versus the number of parameters of the second part of the flow.The number of parameters of the first GRU encoder part is held fixed and not included in the quantity displayed on the x-axis.The second part, which consists of an MLP with two layers, is varied in its hidden dimension.For the Gaussianization flow, additionally the number of flow parameters is varied, which specifically for Gaussianization flows leads to more internal layers and more kernel basis elements (see appendix A for details).In general, the flexible Gaussianization flow is able to reach better density estimation by having a lower KL divergence to the true distribution than the simpler affine flow.As can be seen in Fig. 5 b), the MSE loss with a standard normal posterior is roughly the same scale as the other posterior approximations in dataset 2, and the corresponding KL-divergence to the true posterior only a little worse than for an affine or Gaussianization flow.This happens by chance here, since true posteriors have various shapes and scales and no reason to match a standard normal distribution.For dataset 1 (Fig. 5 a)), on the other hand, the standard normal using the MSE-loss is much worse in terms of KLdivergence to the true posterior and therefore not shown in the figure.The better approximation performance of the Gaussianization flow has to do with the possibility to adjust the potential complexity of the flow itself, i.e. increase the size of F .For an affine flow, on the other hand, F is fixed to be a mean μ and width σ, and one can only increase the complexity by increasing the MLP hidden dimension.In all cases, the final performance saturates, where more parameters do not help performance.
The results suggest that flow approximation capabilities are not the bottleneck in this toy study.As mentioned earlier and indicated in Fig. 4 a), for a single detection module the posterior shapes are highly non-Gaussian.Yet the Gaussianization flow has a good approximation quality, i.e. a KL-divergence close to zero (Fig. 5 a).In dataset 2, which consists of four detection modules, the posterior shapes tend to be more Gaussian due to the larger amount of photons per event (Fig. 4 b), so a close posterior approximation should be easier to achieve than for dataset 1 when using the same normalizing flow.The opposite is the case, however, as seen in the slightly higher KL-divergence offset from zero (Fig. 5 b).These results suggest that either more training data is required for this more complex dataset or the encoding architecture needs to be improved -or a combination of both.We come back to this at the very end of the paper.

Coverage of the neural network posterior
In Frequentist or Bayesian analyses it is important to know how well confidence or credible intervals cover the true values.In a Frequentist analysis, in the limit of large amounts of data, Wilks' Theorem [49] implies that the quantity λ = −2 • (lnL(θ t ) − lnL(θ 0 )) is χ 2 ν distributed with ν corresponding to the dimensionality of θ when θ 0 is the value that maximizes the likelihood function and θ t the true value that generated the data.It is connected to the fact that likelihood scans around the optimum have an approximately Gaussian shape in the large-data limit.A similar statement appears in a Bayesian analysis from the Bernstein-von-Mises theorem [48], which implies that for large amounts of data the posterior p(θ; x) becomes Gaussian.For any n-dimensional multivariate Gaussian the quantity (x − µ)C −1 (x − µ) is χ 2 n distributed [50], which again implies the quantity λ Bayes = −2 • (lnp(θ t ; x) − lnp(θ 0 ; x)) is also χ 2 distributed in the Gaussian limit of p(θ; x).Both of these coverage calculations require the assumption of the large-data limit, but are very efficient because they do not require a large numerical effort.With normalizing flows that have a Gaussian distribution as base distribution, we can use a similar methodology to calculate coverage probabilities for contours from a target distribution of any shape and dimension without resorting to numerical integration.The resulting coverage results are obtained for specific unique credible intervals that are defined by the base distribution.This is explained in the following.Figure 6 a) illustrates that central simply connected3 credible intervals at the base distribution transform to simply connected intervals in the target space because of the diffeomorphism connecting the two spaces.The normalizing-flow relation (eq.40) ensures that the interval in the target space still covers the same probability mass.This interval will be called "base-ordered" in the following, as it is the unique interval obtained using an ordering principle in the base PDF.If no ordering principle is applied, there are infinitely many ways to find an interval that contains a given probability mass.A similar issue arises in Frequentist confidence intervals, where a unique interval can be obtained by a likelihood ordering principle [51].Importantly, the base-ordered interval is always simply connected and in general different from the unique interval obtained by ordering the probability mass directly in the target space, as indicated in the right plot in Fig. 6.As such, for a given normalizing-flow there are always two unique credible intervals: a "base-ordered" interval that is transformed into the target space, and a "target-ordered" interval that is directly constructed in target space.The construction of the target-ordered interval requires numerical integration, in particular to calculate coverage probabilities for these intervals.In contrast, the base-ordered interval can be analytically calculated in base space utilizing the Gaussian base distribution, and transformed via the normalizing flow mapping interval to the target space.This in principle works for any dimension or target shape.More importantly, coverage probabilities can also be analytically calculated by utilizing a statistical relationship between the Gaussian distribution and the χ 2 distribution.This is indicated in Fig. 7 a).If the samples z o,i follow the target q ϕ (z o ), the corresponding samples ẑi at the base must follow a standard normal distribution, and this implies the quantity . This fact can be used to calculate coverage probabilities of the base-ordered contours, which by probability conservation (eq.40) is also valid for the transformed base-ordered intervals in target space.The same idea works distributions on manifolds like spheres.In [52] the authors discuss how to define a flexible normalizing flow on a sphere via stereographic projection of a normalizing flow in the plane.A more stable alternative turns out to be directly to start with a base distribution on the sphere and parametrize a flexible flow which is intrinsic to the manifold [43].We can combine both of these ideas to define a flexible flow on the sphere that also allows to define coverage.The methodology is illustrated in Fig. 7 b) and consists of multiple sub-flows.
The base distribution is again a standard normal distribution in R n .It is transformed to another distribution in R n that itself corresponds to the flat distribution on the n-sphere once it is stereographically projected.Then follows the stereographic projection, and finally the intrinsic flow on the sphere.The first two flows from the standard normal in R n to the flat distribution on the sphere can be combined without actually invoking a Riemannian manifold flow which normally would involve more complicted Jacobian factors4 (see eq. 42).To do so, the distribution in the plane is first written in spherical coordinates.The radial coordinate of the Gaussian distribution in R n is then transformed to θ d−1 on the sphere, while the remaining angular coordinates (θ 1 , . . ., θ d−2 , ϕ) stay unchanged.The end result of these transformations turns out to be (see appendix B for a derivation) for the special cases of the 1-sphere 2-sphere, respectively.For higher-dimensional spheres, no analytical function for the radial part can be written down, but individual sub-parts of the flow are analytical.Evaluations and inverses can be obtained using bisection and potentially Newton iterations if derivatives are desired.An issue that potentially arises with base-ordered credible intervals, in particular those defined on the sphere, is illustrated in Fig. 6 b) and c).
Because the number of possible transformations for a given flow function class is often large, there can be many different ways to generate similar target-ordered contours, which correspond to different relative base-ordered contour alignments.In Fig. 6 b) a spherical normalizing flow based on the recursive flow in [43] is trained to describe a shape of the letter "M" on the sphere.After training, the 10% and 30% base-ordered contours enclose regions of low PDF values, and are misaligned with the respective target-ordered contours.In the extreme case they can even be on the opposite side of the sphere.The exact numerical amount of over-and under-coverage for the two contour types is in general always slightly different, but in such a situation, a coverage calculation that would indicate over-coverage with respect to the target-ordered contours can simultaneously indicate under-coverage with respect to the baseordered contours.Only when coverage is exact, i.e. the contours enclose the true values exactly as predicted, the two contour definitions agree.
In practice however, one often has slight over and under coverage, and one wants to ensure that both definitions at least agree on the sign and that their "center of gravity" overlaps.In order to achieve this, one can add an extra regularization term to the loss function during training.In Fig. 6 c), the same normalizing flow is trained on the same target shape, but now with an added spherical contour regularization term R c Here x mean is the mean of the target PDF, x base,mid the base mean in target space, and x i,base,50 a point i on the 50% base contour transformed to target space.The first term, a dot product of the target mean with the center of Illustration of a stable normalizing flow q Ψ on the d-sphere starting with a d-dimensional standard normal distribution to allow for exact coverage calculation.The flows ρ 1 and ρ 2 involve only the radial coordinate due to rotational symmetry.The second flow ρ 2 is equal to the stereographic projection with a hyper plane that splits the sphere into two hemispheres.This is different from the visualization for illustrative purposes.contours, adjusts the center of the base-ordered contours to be aligned with the mean of the distribution.It avoids the situation where the base-ordered contour is aligned antipodal to the target-ordered contour.The second term, the variance of dot products of the target mean with N points on the 50% baseordered contour, adjusts the 50% base contour in target space to symmetrically surround the mean of the target distribution.In principle, several regularization terms are possible, but the above terms work reasonably well.In Euclidean space, dot products would be replaced by Euclidean distances.If the PDF is highly degenerate and the target-ordered contours consist of disjoint sets, for example if the PDF consists of two narrow peaks that are widely apart, even regularized base-ordered contours might not be reliable in the sense of the same sign as the target-ordered contours -for certain probabilities base-ordered contours might indicate overcoverage while the corresponding target-ordered contours might indicate undercoverage, or vice versa.In these situations one either has to resort to numerical techniques to evaluate target-ordered contours, since those are usually the ones of interest, or directly work with the base-ordered contours.For visualization purposes, an issue might arise with base-ordered contours when the PDF has more than 2 dimensions.In this case one typically uses 1-d and 2-d marginal views of the PDF, and base-orderend marginal contours are not trivial to compute, while target-ordered contours can simply be constructed using samples, as is for example done in Fig. 11 b) for a 4-d PDF.It is however possible to construct base-ordered contours for conditional sub-parts in an autoregressive PDF.It might be possible to create a resulting marginal base-ordered contour with the correct coverage via an appropriate combination of such conditional base-ordered contours, which could then be included in visualizations of marginal distributions.We leave this for future work.
We can test coverage in a coverage plot as indicated in Fig. 8 which shows expected versus actual coverage probabilities of base-ordered contours for a 3-d posterior over position and direction and the corresponding validation loss curve.The posterior in this example is trained on dataset 3, in which the position, direction and energy of neutrino interactions are different from event to event.We use a factorized posterior q ϕ (x p , y p , θ; x) = q ϕ1 (x p , y p ; x) • q ϕ2 (θ; x p , y p , x), which allows to learn the spherical part separately from the Euclidean part in a stable manner via an autoregressive structure [44].While the Euclidean part again uses a 2-d Gaussianization flow, the spherical part employs the previously described strategy to reach a flat distribution on the 1-sphere (eq.43) and afterwards uses convex Moebius transformations [43] parametrized by a neural network as an intrinsic flow.Coverage is calculated using the base-ordered contours.This is an example for coverage of a joint posterior defined on R 2 × S 1 , for which a numerical calculation of target-ordered coverage is already computationally too time-consuming.The base-ordered coverage calculation, however, would work analogously for higher dimensions and on more than two manifolds without a noticable computational increase.
We can observe interesting coverage behavior during training, as depicted in Fig. 8. Once the first phase of training is over, in which the loss function decreases rapidly, good coverage seems to be already reached, even though the overall training has not finished yet.The second phase, in which the optimization is much slower has been called "random diffusion phase" by [53].The authors argue the network learns a better compression of the input data in this stage.In the context of normalizing flows, it seems the first phase brings the learned posteriors in line with the labels to reach proper coverage.In the second diffusion phase, the posterior regions shrink as much as they can while maintaining coverage.This is an empirical observation, and we leave more in-depth studies for future work.

Systematic Uncertainties
A standard practice to incorporate systematics in both Frequentist and Bayesian analyses is to first assume they can be parametrized by a parameter ν which follows a statistical distribution, for example a Gaussian with a known mean and width, e.g.p(ν; µ ν , σ ν ).Such a parametrization does not make sense in a Frequentist sense, but should rather be understood subjectively as the ignorance about the true value of the systematic parameter ν in question.In a Frequentist analysis, usually every systematic is then included as a log-probability penalty term in the likelihood function with respective nuisance parameters in a profile likelihood approach [54].In a Bayesian analysis the systematic uncertainties influence the joint distribution of x and θ, i.e. p(x, θ) → p(x, θ; ν).A marginalized joint density can then be obtained as P t,M (x, θ) = p(x, θ; ν)p(ν)dν.Because it applies to the joint distribution, it also applies to any term in Bayes' theorem, P t (θ; x, ν) = P t (x; θ, ν)P t (θ; ν) and one can for example obtain the marginalized posterior as which now includes the extra uncertainty about the systematics parameter ν.
In the previous sections we discussed how supervised learning and VAEs approximate the underlying distributions of interest via variational inference.A Monte Carlo simulation that first draws systematic parameters ν and then records events with labels θ and data x according to the detector response will automatically produce samples from P t,M (θ, x) or any true conditional distribution of interest.These are the distributions that the neural network will approximate.Incorporating systematics is therefore possible at no change of the underlying procedure, as long as the simulation one is using for training additionally samples from the systematic distributions.No explicit form of the PDF of the systematic parameters has to be known, only sampling is required.For supervised learning, for example, one can just replace any term P t by P t,M in eq.7 and the neural network based approximation q ϕ will approximate the marginalized true posterior P t,M (θ; x).It is already common practice in modern high-energy physics experiments to generate Monte Carlo simulations with potentially marginalized systematic parameters [55], so there is often no computational overhead in actually including those uncertainties.
To illustrate the effect of marginalization over systematic uncertainties, Fig. 9 shows actual versus expected coverage of the contours of a posterior over the logarithmic energy log 10 (E ν /GeV).The systematic uncertainty in this context is the expected number of observed photons given a certain energy deposition, i.e. the light yield, which is a common systematic uncertainty in neutrino detectors.In this example, it is modeled as a flat prior on the light yield with varying width with up to 50% relative uncertainty.For a given prior width, the normalizing flow is trained on dataset 4, which  includes 4 degrees of freedom -position, direction and energy -and for each event the per-event light yield is additionally drawn with a random draw from the prior.As shown in Fig. 9 a), the wider the systematic prior that is used in the simulation the higher is the over-coverage of the true values.Fig. 9 b) shows the different PDFs for an example event from the test dataset.As expected, the PDFs get wider and have larger entropy for a larger width of the light-yield prior that is used during training.In general, the result from base-ordered and target-ordered contours is pretty similar, which is reflected in the PDF shapes which are unimodal and pretty close to Gaussian.For 30% and 50% systematic uncertainty though, the PDFs become slightly non-Gaussian and the respective target-and base-ordered coverage probabilities start to diverge, as expected.

Goodness-of-Fit of the neural network model
The coverage calculation is only a meaningful information for data that is similar to the training data.If new input data is given to a fully trained neural network, for example real data instead of a Monte Carlo simulation, it is desirable to test if this new data can be described by the neural network model.If this is the case then the coverage results might apply 5 .Traditionally, a test between a model and data is done via a goodness-of-fit test.In the following we describe how the extended supervised training procedure described in Section 2.3 allows to calculate a goodness-of-fit via posterior predictive checks.The structure of the extended supervised model is shown in Fig. 10.On top of a posterior q ϕ , we have an additional posterior q φ over unobserved Fig. 10: Generic structure of the extended supervised loss for goodness-of-fit tests and semi-supervised applications based on four normalizing flows.The respective base space Gaussian random variables ẑ are denoted by circles and the respective flow-defining functions by ρ.The arrows from the variational distributions q to the flow-defining functions ρ indicate sampling, where the samples serve as an input to ρ.
latent variables, and a prior p ψ and p θ as a generative model.All of those are (conditional) normalizing flows, and are crucial to calculate for posterior predictive simulations.Posterior predictive checks are a standard methodology in Bayesian analysis to do goodness-of-fit tests of the resulting posterior distribution (see [56], chapter 6.3).In the extended supervised scheme, a posterior-predictive p-value, sometimes also called "Bayesian p-value" [56], can be defined via x,z I T (x,z)>T (x obs ,z) p θ (x; z)q ϕ (z; x obs )dxdz, (48) and numerically calculated using samples from p θ (x; z) • q ϕ (z; x obs ), i.e. from the posterior predictive distribution.The quantities T (x, z) are called test quantities, and I T (x,z)>T (x obs ,z) the indicator function, where the integral effectively counts the fraction of samples where T (x i , z i ) > T (x obs , z i ) for an infinite number of samples x i , z i .A comparison of a given p-value with the ensemble of p-values from the training dataset gives a goodness-of-fit criterion.In contrast to Frequentist goodness-of-fit testing, the test quantities also depend on the parameter z, and the p-value distribution of the true model can be non-uniform and is usually more peaked around 0.5 [56].Only in certain limits, for example for infinitely precise posterior distributions, does the p-value distribution approach the uniform distribution [56].In our case, we additionally do not expect to reach a uniform distribution as the neural network model will only yield an approximation of the true model, in particular since we do not have access to the exact data-generating function, but learn an approximation of it in tandem with the posterior distribution.We therefore determine the null hypothesis p-value distribution from the training data distribution.We propose to use a normalized version of the logarithm of the data-generating PDF (eq.39) as a test quantity for the Poisson process data.The division by the number of detection modules N and number of total observed photons N tot in the whole detector makes different realizations of the data comparable 6 .For other data-generating PDFs, or a test quantity based on the posterior distribution, this division would probably not be necessary.Any test quantity based on the data-generating PDF or posterior can readily be calculated since they are part of the supervised training and modeled by normalizing flows.If those normalizing flows are bad approximations of the true underlying PDFs, the test quantity will be less constraining.Any improvement in approximation of the normalizing flows is therefore desired for a better overall goodness-of-fit procedure.
In order to illustrate the overall procedure of goodness-of-fit testing, the extended supervised training is performed on dataset 5 (see table 1 for information on all datasets) using the loss function described in eq.36.In principle, the posterior q ϕ could be a model trained in a purely supervised fashion as described in section 2.1, with its weights fixed, and only afterwards the additional normalizing flows are trained using the loss function in eq.36.However, in the following we train everything together from scratch, including the supervised term and the ELBO-like term.For the ELBO-like term, we sample one latent variable z u,i per observed input pair z o,i , x i using the reparametrization trick, and minimize the overall resulting loss function, including L supervised , with stochastic gradient descent.The training evolves similar to supervised training until the overall loss function does not reduce significantly anymore (see appendix A for training details).Dataset 5 contains shower-like interactions in a detector consisting of 400 modules, where the interactions are constrained to lie within a certain containment region (the black dashed region in Fig. 11).The photon arrival time PDFs p(t; z o , x) that appear in the data-generating function (eq.39) are modeled by a 1-dimensional conditional normalizing flow that takes as extra information the position of each module.Additionally, the mean count expectations of the Poisson distributions are predicted, and together with the normalizing flow, allow to describe the complete data-generating function.In the following, we apply this trained model also to events of dataset 6 (uncontained showers) and dataset 7 (track-like events).Figure 11 a) shows the interaction of three example events involved in dataset 5 and 7.The first two events are a low energy (blue) and high-energy (orange) shower event from dataset 5.The third event is a 20 m track-like event (green) from dataset 7 (see table 1 for a list of all datasets).All three events are within the containment volume, as defined by dataset 5 (contained showers), which is also indicated as the black bounding box.The expected value of photon counts per module is shown as vertical bars next to each module.Figure 11 b) shows a comparison of the final posterior distributions of the two shower events over the four degrees of freedom after training, .i.e. a visualization of the learned joint distribution over z o and z u , q ϕ (z o ; x) • q φ (z u ; z o , x).The energy is not used as a label during training, so the latent variable from the unsupervised part of the extended supervised training can be identified with the energy degree of freedom.The high-energy shower (orange) produces more photons, and allows to better inform about the true event properties than the lower energy shower.
The p-value distribution of the test set of dataset 5 after training, calculated using eq.48, is shown in Fig. 12 as the orange curve.In addition, two other p-value distributions are depicted, where the final trained model is applied to events from dataset 6 (shower-like uncontained) and dataset 7 (track-like).As can be seen, the p-value distribution for dataset 5 is not totally uniform, possibly because the individual PDFs have not all fully converged during training.In general, full convergence of all normalizing flows is indicated when the p-value distribution of the training dataset is fairly symmetric.Towards the beginning of the training process, the p-value distribution of the training dataset is still highly asymmetric.The p-value distribution of dataset 6, uncontained shower events, has a similar shape as contained showers, except it also contains a peak for p-values close to 1.The events that pile up close to 1 are events that lie outside the containment region, which are not part of the training dataset and therefore get such a bad goodness-of-fit.For dataset 7, which contains track-like events, there is a skew in the distribution which extends over the whole p-value range and most events have a p-value close to 1. Figure 13 shows the effect of applying a cut that retains events with a p-value smaller than 0.99.The training dataset is almost unchanged, while the cut removes most events outside the containment region for uncontained shower-like events.For track-like events, additionally events    after cut Fig. 13: Effect of the p-value cut (p val < 0.99) on the spatial distribution (upper row vs. center row) and the distribution of the number of detected photons (lower row) for datasets 4 to 6.The center row also indicates the percentage of retained events for the respective dataset.
inside the detector are substantially reduced.The distribution of the number of photons per event in the lower row indicate that the events that remain for the track-like dataset are low energy events that have an intrinsically low information content and are with the given model not statistically separable to a low energy shower-like event.Therefore, a simple energy cut can remove those events and one would end up with a strong separation, in this example a complete background reduction above a certain energy threshold.An alternative is a cut on high-entropy events, which is possible with normalizing flows, since they allow to calculate the per-event differential entropy to arbitrary precision.This can also be seen in Fig. 11 b).The higher energy shower has a visibly lower overall entropy.This example is obtained with toy simulations that only qualitatively model a situation commonly encountered in neutrino telescopes, in particular with a 2-d compared to a 3-d setting.However, with 400 detection modules and the inclusion of relevant effects like scattering, absorption and typical Cherenkov emission features, the data should broadly mimic a real world experiment like IceCube DeepCore [57] both in statistical and computational aspects.The result illustrates that it is in principle possible to train an extended supervised model on a signal dataset, apply the resulting model to any data event, and based on the resulting p-value decide to keep the event or exclude it from further analysis.The remaining events should be similar to events from the training sample, within uncertainties given by the training statistics and any included systematics.A further second cut on an energy variable or the entropy should then be enough to yield a statistically optimal event selection for the given model.If the machine learning model is improved, for example by a better data encoding scheme or a more powerful normalizing flow, the p-value separation and thereby the final event selection should also improve.Besides such an application for data selection, it might also be interesting to study goodness-of-fit results for concrete data to Monte Carlo comparisons or for data anomaly detection.

Related Work
We discussed that standard supervised training can be viewed as variational inference, and variational inference can be viewed as an approximate likelihood-free inference approach.An active topic in likelihood-free inference are neural sequential inference procedures [58].These procedures typically require guided re-simulations for a single event, and eventually should converge to the true posterior.From our perspective, simple supervised learning with normalizing flows, in combination with coverage and goodness-of-fit guarantees as described earlier, might be more suited for certain applications in high-energy physics or astro-particle physics where processing time is valuable -for example when posteriors for large amounts of data are desired.Also single astronomical multi-messenger alerts [10] can benefit.An initial posterior that can be trusted due to the discussed coverage and goodness-of-fit checks could be send to the community quickly.One can later always decide to refine the result with an iterative likelihood-free procedure like automatic posterior transformation (APT) [59] which would naturally blend in due to its use of normalizing flows.For gravitational waves, posterior inference with autoregressive flows has been discussed recently in [17].In this line of work, coverage is checked with 1-dimensional marginal distributions which seems more computationally demanding than our approach.Systematics or goodness-of-fit tests are not discussed in that paper.
Another related line of work discusses inference of physics posteriors specifically based on coupling layers [60], originally dubbed conditional invertible neural networks (c-INNs), more recently using a conditional variant [40] of GLOW [39] where conditionality is added to the coupling layer inputs.For Euclidean normalizing flows in low dimensions, which is what we are interested in, there exist expressive normalizing flows without coupling layers like Gaussianization flows [41].These normalizing flows, similar to certain manifold flows, can be made conditional by predicting their whole parameter set by a neural network, which is what we use in this paper.A further difference to c-INNs is that we interpret the base space not as a latent space, but purely as an auxiliary space that can be utilized for coverage calculations.
Normalizing flows have also been used to model the likelihood function directly [14].Pure likelihood modeling inherits some problems of the standard likelihood machinery: a time consuming optimization or MCMC procedure.However, in our view it is useful to model the likelihood as the data generating PDF in the proposed "extended supervised training" (Section 2.3) in combination with an associated posterior.
A central part of our derivation is based on the joint KL-divergence.The importance of the joint KL-divergence has been emphasized in other approximate Bayesian inference, for example to unify expectation propagation in Gaussian processes [61] and as mentioned earlier for the derivation of the supervised loss for neural networks [19].In the latter example, however, the loss function is not called the log-posterior, but conditional maximum likelihood objective with respect to the neural network parameters.While both views are correct, the posterior interpretation is more fruitful in our opinion, in particular in combination with normalizing flows which shows that supervised learning can behave as rather well-approximated variational inference.Another related recent work applied normalizing flows for a neutrino oscillation analysis [62].The authors do not use the joint KL-divergence in the derivation, which makes it less obvious how to include systematics, although they briefly mention one approach among two alternatives that is similar to our suggestion.The authors also argue for a vague agreement of the base distribution with a Gaussian after training, but do not discuss quantifiable coverage tests.
Likelihood approaches typically only allow to define a rigorous goodness-of-fit for special random variables.However, realistic applications typically involve unbinned likelihoods, for which no rigorous goodness-of-fit can be defined [63].Furthermore, coverage and systematics checks are often computationally challenging.For traditional machine learning based regression, and even for other methods with normalizing flows, as pointed out above, coverage or goodness-of-fit are either impossible or not straightforward to compute.
The proposed extended supervised training with conditional normalizing f ĺows combines 1) systematics 2) coverage and 3) a rigorous goodness-of-fit measure in one variational inference algorithm, without restrictions to the underlying PDFs.We are not aware of any other Ansatz in the literature that has these properties.

Conclusion
In the first part of the paper we have shown that the supervised learning loss, an extended supervised loss and the unsupervised ELBO-loss of variational autoencoders can all be derived under one unifying theme of variational inference.The inference in all these approaches can be traced back to the minimization of the KL-divergence between the true joint distribution of data and labels and an approximation for certain sub-parts of the joint distribution based on neural networks.
For supervised learning, the derivation of the loss function from the KL-divergence of the joint distribution is well-known but the connection to variational inference of the posterior distribution is typically not being made.From our perspective, however, this interpretation is very useful in high-energy physics, even when using an MSE loss function which corresponds to a standard normal posterior approximation.After all, a central task of a physics analysis is to calculate a posterior distribution, and the variational perspective shows that neural networks learn to predict exactly that.Conditional amortized normalizing flows with a standard normal base naturally generalize the MSE loss, and are crucial to make the variational inference perspective usable in practice.
A simple replacement of observed with unobserved labels in the KL-divergence allows to eventually derive the ELBO loss of variational autoencoders.In this derivation, the usual encoder q ϕ and decoder p θ are both manifestly auxiliary, and are explicitly separated from the true distribution P t which p θ and q ϕ both try to approximate within a given identifiability class of transformations.This is different from the usual derivation via the marginal likelihood, e.g. as discussed in the original VAE paper [23], where p θ is often used to denote both the true posterior and the recognition model, which obfuscates the fact that those are typically not the same and also do not have to be the same parametric class.
Conditional normalizing flows with a Gaussian base distribution have a second advantage.They can be used to calculate coverage without numerical integration for certain unique "base-ordered" contours of the normalizing-flow PDF which are defined via the base space.In contrast to "target-ordered" contours created in the target space by starting at the highest PDF values, which can consist of disjoint sets, the base-ordered contours are always simply connected in the topological sense.We further expanded the idea to calculate coverage for distributions on m-spheres S m , or more generally to posteriors defined autoregressively on tensor products like R n × S m .To enable coverage for such distributions we derived the transformation from the m-dimensional standard normal Gaussian to the flat distribution on the m-sphere, which for 1 and 2 dimensions is conveniently simple and invertible.A further intrinsic flow on the sphere can then describe a more complex distribution.This two-step procedure ensures stable training, since the transformation from the plane to the sphere is fixed and does not change during optimization.We also discussed a regularization term that can be added to the loss function in order to ensure that base-ordered and target-ordered contours align.This can be in particular important for spherical distributions.If the distributions are very non-Gaussian with widely separated peaks, the coverage probabilities for base-ordered and target-ordered contours can diverge -which is only a problem when over-or undercoverage is observed.In such a case extra care in the interpretation is necessary.
We demonstrate a coverage calculation for a joint 3-d posterior of position and direction (x, y, θ) for shower-like neutrino interactions in a 2-d toy Monte Carlo.Empirically, we observe that coverage of the base-ordered contours is already obtained once the training phase enters the random diffusion phase as defined by [53].Since the loss still decreases during this phase, only more slowly, and in the variational interpretation this implies shrinking posteriors, in this second training phase the posterior approximations get slowly better while coverage more or less stays intact.
Viewing all these models as coming from the joint KL-divergence makes it also natural to include systematic uncertainties in the final approximated posterior distributions.They can be included via sampling from systematics priors during the Monte Carlo generation process, which lets the normalizing flow approximate the marginalized true posterior distribution -a process one might call "effective marginalization".The training process of the model otherwise stays the same.We test this method with a dataset with a varying energy-scale systematic and an approximation to the energy posterior distribution.The resulting posterior regions are wider and are over-covering more for a higher systematic uncertainty, as expected.
Finally, we test the extended supervised learning scheme.It is related to semi-supervised learning and involves a hybrid of the supervised and VAE loss functions.In addition to the label posterior q φ , a latent posterior q ϕ , a prior p ψ , and a decoder p θ are learned either jointly or separately as an add-on to the supervised posterior q φ .The extra distribution q ϕ learns latent variables to describe information not captured by the posterior of known labels q ϕ , while the prior and decoder define a generative model.The combination of a posterior and a generative model allows to perform Bayesian predictive simulations, which in combination with a test quantity T (x, z), which we suggest to be the normalized logarithm of the learned data-generating PDF, allows to calculate a p-value which can be used as a goodness-of-fit criterion.Any data that are not similar to the training data within the given uncertainty of the statistical approximations of the neural networks lead to p-values close to 1 and can be excluded.We test this concept with a neural network model trained on shower-like neutrino interactions contained in the center of a hypothetical detector consisting of 400 detection modules.In comparison with shower-like interactions in the full detector, the model calculates a small p-value for events on the boundary of the detector.We also apply it on a dataset containing track-like energy depositions mimicking muon tracks, for which the separation is markedly visible.Such a goodness-of-fit selection strategy can potentially allow a very simple but effective data selection.
Let us mention two aspects that require further investigation.The first is related to the stability of the training algorithm and the final learned distributions.We used stochastic weight averaging to average over several iterations at the end of training to obtain a more stable posterior solution.This serves as an approximation to ensemble averaging, which is basically averaging under a variational posterior over network weights.Due to varying learning rates for different model complexities, slight fluctuations in the final loss values remained (see Fig. 5).Ideally, one would like a scheme where repeated training schedules lead to repeatable posterior regions which are indistinguishable within a certain margin which can be defined before training.This needs to be studied more systematically.
A second aspect is related to the data encoding.In all studies we have split the conditional posterior into a data encoding based on a GRU and aggregate MLP, and a second MLP which further maps that encoding summary to the normalizing-flow parameters.Our Monte Carlo tests show that at some level of flow complexity the performance bottleneck of the posterior approximation can in part stem from the data-encoding architecture.This particular encoding was only chosen for simplicity and serves as a proof-of-concept, but it is not an optimal encoding.In an application of the method, it is advised to use the best encoding mechanism available for the given experiment.Any improvement in data encoding not only translates to tighter posterior contours, but also to a more powerful goodness-of-fit calculation.

A Implementation details A.1 Training
All training procedures use ADAM [64], which gives much faster convergence than standard gradient descent (SGD) [65].For the experiments in section 5 we use a mostly fixed learning rate which adjusts itself to keep fluctuations of the validation loss at a certain level.For the affine flows, we had to fix the learning rate completely, because too large learning rates would sometimes lead to drastic jumps in the loss function.This might be a problem that only occurs for flows which are defined by less than a handful of parameters.For Gaussianization flows [41], similar behavior could be seen for single-layer flows with a single basis function.Once the number of flow parameters is increased beyond a certain level, it does not seem to be an issue anymore.For the experiments in section 6, 7 and 8 we start with a high learning rate and then use a decreasing learning rate scheduler that lowers the learning rate by fixed factors until it reaches a value of about 0.0001.Towards the end of training, once the loss reaches a more or less constant level, the final iterations are averaged using stochastic weight averaging (SWA) [46].As argued by the authors, SWA approximates model ensemble averages.Model averages, in turn, can be used to calculate the expectation value under the approximate weight posterior since SGD iterations itself can be interpreted as samples from a variational approximation of the weight posterior [66].

A.2 Posterior
The posterior parametrization is indicated in Fig. 3 (c).The first encoding part is similar for all experiments.A GRU encodes the input data sequentially, the output is aggregated and mapped by a single layer MLP to a representation h that summarizes the data.
The second step in the posterior is another MLP that maps the summary representation h to the respective flow parameters.We use two inner layers for this second MLP, and vary its number of dimensions for the various experiments.The parameter choices are described in Section A.4.

A.3 Data-generating PDF in goodness-of-fit
For the data-generating PDF for the extended supervised model used for the goodness-of-fit, we have a discrete Poisson part and a continuous part.We use a normalizing flow for the continuous part p(t; z o ) of the generating distribution (eq.39).The involved MLP that is used in the amortization also outputs a logarithmic mean expectation for photon counts to calculate the related Poisson distribution, so we model the discrete and continuous part in a joint prediction step.By simultaneously predicting the flow parameters and the Poisson mean, it is possible to model the correlation between the two.The input to the MLP are the event parameter labels z o and the module position in (x,y) coordinates.

A.4 Parameter choices
We choose a single-layer GRU with hidden dimension 10 and for the aggregation MLP choose intermediate dimension 15 and map to a 20-dimensional representation h for the first 3 datasets in the paper.This encoding part then has roughly 1000 parameters.For the larger detector used in datasets 5-7 the involved dimensions are a bit larger.For the comparison in Fig. 5 the dimensions in the two inner layers of the second MLP vary between 1 and 100.For the Gaussianization flow, additionally the number of flow layers varies between 1 and 5, and the number of basis elements per flow varies between 1 and 10, as the flow complexity increases.

B Transforming a Gaussian to the flat distribution on the sphere
The following calculation derives the transformation ρ tot from the d-dimensional Gaussian distribution to the flat distribution on the d-sphere (see Fig. 7 b)) which is discussed in Section 6.The flow can be split up as ρ tot = ρ 2 • ρ 1 .The first flow ρ 1 transforms the standard normal to a distribution p f that corresponds to the distribution in the plane that is the stereographic projection of the flat distribution on the sphere.We can derive p f by starting with the flat distribution on the sphere.
The flat distribution on the d-sphere is defined as [67] p flat,sphere = sin(θ 1 ) 1 • . . .f +1 , which is expressed here entirely as a connection to the radial coordinate in the plane via r f = x p,j .It therefore follows that or vice versa which is a relation that is indicated in Fig. 7 b).
Applying the flow ρ −1 2 to the flat distribution on the sphere we obtain where p f is now the corresponding PDF in the plane after the stereographic projection.We now need to find the flow ρ 1 (indicated in Fig. 7 b)) to transform a standard normal Gaussian distribution to p f or vice versa.Similar to the transformation from the sphere to plane, this transformation can be done entirely in the radial coordinate when the Gaussian distribution is written in spherical coordinates.The radial transformation can be calculated using the cumulative distribution functions of the radial PDFs.The two radial PDFs are where p r,g is the radial distribution of the d-dimensional standard normal distribution, also known as the χ-distribution, and p r,f the radial distribution of p f .The corresponding CDFs which map from R + to [0, 1] then follow to be8  where Γ(x, y) is the upper incomplete Gamma function and 2 F 1 is the Gauss hypergeometric function.The transformation ρ 1 then can be written with the corresponding CDFs as r f = ρ 1 (r g ) = CDF −1 r,f (CDF r,g (r g )), (64) which can in general not be written down analytically except for d = 1 and d = 2. Using bisection and Newton iterations it is possible to evaluate this expression and its inverse for higher d.The Newton iterations make the result differentiable, which is required if it is to be used in variational autoencoders.Finally, the transformation ρ tot = ρ 2 • ρ 1 can be written as ρ tot,2 (r g ) = arccos 1 − 2 • exp(−r 2 g /2) , (d = 2) (67) which has a simple and invertible structure for d = 1 and d = 2 after some manipulation with trigonometric identities, while higher dimensions again require bisection and Newton iterations.Because ρ tot defines a flow from the d-dimensional standard normal distribution to the flat distribution on the d-sphere, besides its usage in normalizing flows for coverage as describe in Section 6, it can be used as a non-standard way to generate uniform samples on the d-sphere.This can be done by first sampling from the d-dimensional standard normal distribution, switching to spherical coordinates, and finally transforming the radial coordinate to the last coordinate θ d−1 on the sphere using ρ tot , while keeping all other angles θ 1 , . . ., θ d−2 and ϕ as they are.

Fig. 2 :
photon collection module log.expected photon count [a.u.] MLP(c) Common data encoding for all experiments.

Fig. 3 :
Fig. 3: General conditional flow (a) and an affine conditional flow (b) parametrization of the approximate posterior in supervised learning.Choosing σ ϕ = 1 yields a shifted standard normal distribution which is the PDF used in the MSE loss.General normalizing flow parameters are denoted by ⃗ F and the parameters of the encoding neural network are denoted by ϕ.The common encoding scheme for all experiments is depicted in (c), which amortizes the parameters ⃗ F .

Fig. 4 :
Fig. 4: A comparison of posteriors of the position for the example events A and B from dataset 1 and events C and D from dataset 2. The normalizing flow posterior is shown together with a 68% contained probability mass in black.The 68% probability mass contour of the true posterior assuming a flat prior is shown in white.The true event positions are marked in red.The upper row shows the result for Gaussianization flows, the lower row for an affine flow (a Gaussian) with a single covariance parameter.
Performance vs. complexity for dataset 1.
Performance vs. complexity for dataset 2.

Fig. 5 :
Fig.5: Posterior approximation performance of a Gaussianization flow[41] (GF), an affine flow with a single variable width σ (affine), and an affine flow with σ = 1 (MSE).Along the x-axis the number of parameters and the potential flow complexity increases, while the encoding complexity is held fixed (see text).The respective upper plot shows the validation loss.The black dotted line shows the loss obtained from the true posterior.The respective lower plot shows the sample-based KL-divergence.

Fig. 6 :
Fig. 6: Illustration of base-ordered contours (red) and target-ordered contours (black) for a 1-d normalizing flow (a) and for a spherical normalizing flow (b)+(c) with different coverage probabilities.The spherical example illustrates that the base-ordered contour can be misaligned with the target-ordered contour.The spherical examples also show the inverse transformations of target-ordered contours in the base space.19

Fig. 7 :
Fig. 7: Coverage schematics to indicate the connection to the χ 2 distribution at the base distribution (a) and to indicate the discussed mapping strategy to calculate coverage for spherical distributions (b).

Fig. 8 :
Fig. 8: Illustration of the coverage behavior for a posterior calculated with dataset 3 (a).The coverage is calculated for different stages of training (iteration in the legend), including the end of training where the model is averaged using SWA.The validation loss curve is shown in (b), where the iterations in fig (a) are indicated by the shaded gray area.

Fig. 9 :
Fig. 9: Illustration of the effect of systematic uncertainty inclusion via effective marginalization.Each model is trained on data with a different flat prior on the overall light yield per energy deposition (in percentage).The test dataset has no light yield uncertainty included.The left figure (a) shows actual coverage versus expected coverage for base-ordered and target-ordered contours using the test dataset.The right figure shows posterior distributions for an example event from the test dataset.
(a) Example interactions of three events overlaid on the detector used in datasets 5 − 7. Detector modules are denoted as black dots.The relative expected photon count in each module is visualized as vertical bars with the respective color.
var. [a.u.](b) Visualization of the final posterior approximation (marginal PDF in 1-d and target-ordered contours of the marginal PDFs in 2-d) for the two shower-like events visualized in Fig. 11 a).The Monte Carlo truth is indicated with markers and vertical lines.

Fig. 11 :Fig. 12 :
Fig.11:A few events from the datasets for the goodness-of-fit example.

Table 1 :
Properties of the datasets used for various comparisons.The detector geometry always consists of 1 or 16 modules as depicted in Fig.