  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, coverage,
systematic uncertainties or a goodness-of-fit measure 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 autoencoders (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 to
calculate coverage probabilities without numerical integration for specific
"base-ordered" contours that are unique to normalizing flows. Furthermore we
show how systematic uncertainties can be included via effective marginalization
during training. The proposed extended supervised training incorporates (1)
coverage calculation, (2) systematics and (3) a goodness-of-fit measure in a
single machine-learning model. There are in principle no constraints on the
shape of the involved distributions, in fact the machinery works with complex
multi-modal distributions defined on product spaces like $\mathbb{R}^n \times
\mathbb{S}^m$. The coverage calculation, however, requires care in its
interpretation when the distributions are too degenerate. We see great
potential for exploiting this per-event information in event selections or for
fast astronomical alerts which require uncertainty guarantees.
