1.1. Bayesian deep learning

In theory, the uncertainty of $ w $ in an neural network can be derived directly from Bayes Theorem, which is the basis of Bayesian Deep Learning (MacKay, Reference MacKay1992; Neal, Reference Neal1996; Andrieu et al., Reference Andrieu, de Freitas, Doucet and Jordan2003). From Eq. (2), if we choose $ A=w $ and $ B $ as the training data, we know that the posterior uncertainty in $ w $ given a training dataset is determined by the prior uncertainty in $ w $ , and a likelihood that indicates how likely the training data are if the weight vector is equal to a certain $ w $ .

Since the function $ f(..) $ in deep learning is highly nonlinear, this likelihood, and hence the uncertainty pdf for the weight vector $ w $ , will have an extremely complex shape (e.g., many local minima and maxima). For low-dimensional systems, one can use methods such as Markov-Chain Monte Carlo (MCMC), Langevin sampling, and Hybrid or Hamiltonian MC to accurately describe the posterior pdf (Andrieu et al., Reference Andrieu, de Freitas, Doucet and Jordan2003; Van Leeuwen et al., Reference Van Leeuwen, Cheng and Reich2015), but these methods do not scale well for high-dimensional problems, as is often the case in deep learning applications. To alleviate this issue, approximations to the posterior pdf have been proposed, such as the Gaussian or Laplace approximation, or variational Bayes (approximated pdf by a finite set of pdfs of a certain shape, e.g., a Gaussian mixture). These methods need the computation of the Hessian, the second derivative of the loss function, which can be computationally expensive. Apart from the crude Gaussian approximation, all these methods remain very expensive in high-dimensional problems, which pose difficulties in practice (Blundell et al., Reference Blundell, Cornebise, Kavukcuoglu and Wierstra2015).

1.2. Bagging and bootstrapping

A popular machine-learning method for uncertainty quantification related to the weights $ w $ is bootstrapping or Bagging Breiman (Reference Breiman1996), which uses an ensemble of weight vectors to generate an ensemble of outputs. In this method, the weight vector ensemble can be generated in several ways, e.g., by randomly drawing initial weights and by randomly varying the order in which training batches are used during training of the neural network, or by giving each neural network a subset of the training data. When a new input vector is presented to the neural network, each trained weight vector in the resulting ensemble is used to generate one output value and interpreted to have the same probability. This, then, results in an ensemble of output values that is suggested to represent its uncertainty pdf. Unfortunately, while this method gives a spread in the output values, this spread does not represent a proper uncertainty estimate. The reason is that during training, the ensemble of initial weight vectors is transformed into an ensemble of trained weight vectors. Because this transformation is nonlinear, the pdf of the equally probable initial weight vectors is deformed into a pdf with a complicated shape determined by the likelihood of each weight vector to produce the training data (corresponding to the $ l\left(B|A\right) $ factor above). Hence, the trained weight vectors are not equally likely anymore, and neither are the output values based on the neural network’s that use these weights. The simple treatment of equally likely output values, as assumed in Bagging, is inappropriate.

This issue with the likelihood of the weight vectors is more troublesome than might be anticipated. In principle, the likelihood value of each trained weight vector, also called the importance weight, can be properly calculated and included to evaluate the ensemble of weight vectors. However, in practice, one weight vector typically fits the training data significantly better than the others, resulting in a much higher importance weight. This dominance effectively collapses the ensemble to a single neural network, eliminating the diversity needed for meaningful uncertainty estimation and yielding an unrealistically low (often near-zero) uncertainty. While this challenge is acknowledged in Bayesian Deep Learning, where methods have been developed to construct ensembles of approximately equally likely weight vectors, it is overlooked in Bagging and related approaches, essentially making the uncertainty quantification useless. We return to this critical issue in detail in Section 3.3.

1.3. Monte-Carlo Dropout

Another popular method for uncertainty quantification related to the weights $ w $ is Monte-Carlo (MC) Dropout (Gal and Ghahramani, Reference Gal and Ghahramani2016). In conventional dropout, a method used to avoid overfitting, the weight of each neuron is set to zero during the training phase with a probability $ P $ . After training, all weights are used, but their trained value is multiplied by $ P $ to generate one neural network. In Monte-Carlo Dropout, the same dropout procedure used during training is also used during prediction. In other words, the trained weights are also put randomly to zero with probability $ P $ during prediction. This prediction procedure is performed several times to generate an ensemble of weight vectors, leading to an ensemble of output values. This output ensemble is then suggested to provide an uncertainty estimate.

The issue with MC Dropout is that during training one realization of a trained weight vector is created. This weight vector is randomly perturbed by dropout to generate an ensemble of perturbed trained weight vectors. What should be the next step is to importance-weigh each of those perturbed weight vectors by their likelihood value, as discussed above for Bagging, but that step is omitted. Hence, similar to Bagging, the uncertainty quantification is inappropriate. Furthermore, it fails a fundamental requirement of Bayesian inference: the posterior uncertainty should decrease as more data becomes available. This behavior has been shown to be absent in MC Dropout (Osband, Reference Osband2016; Verdoja and Kyrki, Reference Verdoja and Kyrki2021, and references in these papers).

1.4. Estimating uncertainty distribution parameters directly

The above discussion focused on uncertainties in the weight vector $ w $ . A second category of uncertainty-quantification approaches tries to estimate the pdf of output uncertainty directly. Instead of predicting output $ z $ , the neural network predicts:

(4) $$ {p}_{\alpha }(z)={f}_{\alpha}\left(x,{w}_{\alpha}\right) $$

where $ {p}_{\alpha }(z) $ denotes the uncertainty pdf of output $ z $ with parameters $ \alpha $ , and the subscript $ \alpha $ on the weight vector $ w $ and on $ f $ denotes that the target of estimation are the distribution parameters $ \alpha $ , not $ z $ itself as in Eq. (1). The difference from the methods above is that the uncertainty pdf in $ z $ is not determined via an ensemble of weight vectors $ w $ , but instead, one weight vector is trained to provide the pdf parameters directly. Often, a Gaussian assumption is used, in which the mean and variance are estimated; however, other standard distributions can also be used (Nix and Weigend, Reference Nix and Weigend1994; Wang et al., Reference Wang, Li, Ouyang and Wang2016).

1.5. Evidential deep learning

Evidential deep learning for deep regression is an extension of the above methodology in which not only the parameters of the pdf of $ z $ are learned, but also the uncertainty in those parameters, e.g. Amini et al. (Reference Amini, Schwarting, Soleimany and Rus2020). Specifically, a Normal-Inverse-Gamma distribution is used to describe the uncertainty in the mean and variance parameters. In this way, part of the epistemic or structural uncertainty is addressed.

1.6. Quantile regression

While the two methods discussed above assume a certain shape for the uncertainty pdf, typically Gaussian, in quantile regression the pdf is allowed to have any shape, at least in principle. The neural network predicts:

(5) $$ {P}_z={f}_P\left(x,{w}_P\right) $$

where $ {P}_z $ denotes predefined quantiles of the uncertainty pdf of output $ z $ , and the subscript $ P $ on the weight vector $ w $ and on $ f $ denotes that the target of estimation iare these quantiles (Wang et al., Reference Wang, Li, Ouyang and Wang2016; Pfreundschuh et al., Reference Pfreundschuh, Eriksson and Duncan2018; Sønderby et al., Reference Sønderby, Espeholt and Heek2020). Quantile regression works well for a low-dimensional output space but becomes impractical in high-dimensional problems, because the number of output variables required to describe histograms increases dramatically with the dimension. To avoid the latter issue, this method is often used to estimate the marginal uncertainty pdfs instead of the full joint pdf.

1.7. Deep Ensembles

However, the methods in the second and third categories fail to account for the uncertainty in the weight vector $ {w}_p $ . This omission is partly corrected in so-called Deep Ensembles (Lakshminarayanan et al., Reference Lakshminarayanan, Pritzel and Blundell2017). Deep Ensembles generates an ensemble of pdf parameters, or, alternatively, an ensemble of histograms, by training multiple networks, each initialized with different untrained weight vectors and using distinct training batch orders. The resulting ensemble of pdf parameters is further averaged to represent the pdf of $ z $ , or the histograms themselves are averaged. Similar to Bagging and MC dropout, the simple average of the ensemble overlooks the fact that the trained weight vectors $ {w}_p $ , and thus the output ensemble members, are not equally likely.

As mentioned above, the literature has concentrated on quantifying uncertainty in the weight vector $ w $ or tried to estimate the uncertainty pdf of $ z $ directly. None of the machine-learning methods discussed above has considered the uncertainty in the training and testing data, or in the new input data. While uncertainties in the training, testing, and new input data might be small in some fields, they can be significant in the geosciences, where real-world data, such as direct measurements, model outputs, or reanalysis data, often contain substantial errors that cannot be ignored.

1.8. Conformal prediction

A completely different category is conformal prediction Balasubramanian et al. (Reference Balasubramanian, Ho and Vovk2014). It is based on the concept of dissimilarity, or nonconformity, of a new prediction compared to previous input data and their predictions from the neural network, so that the prediction interval size is larger for new “unusual” instances and smaller for “typical” instances. First, the neural network model is trained on a dataset. Then, a nonconformity score is defined by the user, which measures the extent to which a data point deviates from the model’s predictions. For instance, one can use the difference between the predicted value and the true value as a score. Given new input data, conformal prediction uses the nonconformity scores from the testing data to calculate a prediction interval. The idea is to find the smallest set or interval such that, with a predefined confidence level (like 90%), the true value falls inside. Hence, the method produces a confidence interval, not an uncertainty estimate as Bayesian methods do, and it ignores uncertainty in the weight vector. Furthermore, in the geosciences, the training, testing, and new input data all contain uncertainties, and especially including uncertainty in output training data is not straightforward with this methodology. In fact, one typically relies on Bayesian methods to do this, resulting in an interesting, but perhaps ad hoc mixture of Bayesian and significance testing methodologies.

1.9. Recent reviews on uncertainty quantification

To conclude, while many methods have been developed for quantifying uncertainty in the output of a deep-learning neural network, they are either incomplete or statistically inappropriate. Recent articles by Abdar et al. (Reference Abdar2021), Gawlikowski et al. (Reference Gawlikowski2023), and He et al. (Reference He, Jiang, Xiao, Xu and Li2025) provide comprehensive overviews of existing techniques, but they lack critical discussions of the statistical properties of these methodologies. As a result, understanding the true meaning of uncertainty estimates from these approaches remains challenging. The review by Cheng et al. (Reference Cheng, Quilodrán-Casas and Ouala2023) discusses the connection between data assimilation, machine learning, and uncertainty quantification, but the uncertainty quantification methods mentioned above are quoted without a critical review.

1.10. Our contribution

In this paper, we present a novel approach that accounts for all sources of uncertainty. We derive an expression for the exact pdf and provide an efficient methodology for practical computation. In Section 2, we outline how to systematically incorporate all uncertainty sources. In Section 3, we introduce a practical approach to calculate the relevant uncertainties and combine them effectively. Specifically, we explain how proposal densities, a standard method in statistics, can be used to address the likelihood issue associated with trained weight vectors $ w $ . In Section 4, we demonstrate the methodology on a simple toy problem and on an application to predict cloud process rates and offer a physical interpretation of the resulting uncertainty. Finally, we conclude with a discussion in the last section.

3.1. Source 1: Uncertainty in the new input

The influence of uncertainty in the new input vector $ x $ on the output $ z $ is the most straightforward to implement. One can draw $ {N}_x $ samples from the probability density function $ p\left({x}_t|x\right) $ ,which describes the uncertainty around the observed input vector $ x $ . Each of these sampled input vectors is then propagated through the neural network to generate corresponding samples of the output $ z $ . This Monte Carlo approach provides an empirical approximation of the distribution of $ z $ induced by the uncertainty in $ x $ . Mathematically, this means that we draw the samples from the uncertainty pdf, $ p\left({x}^t|x\right) $ , as:

(13) $$ p\left({x}^t|x\right)=\frac{1}{N_x}\sum \limits_{i=1}^{N_x}\delta \left({x}^t-{x}_i^t\right), $$

where $ {x}_i^t $ are the drawn samples centered around $ x $ . Using this expression, we can rewrite Eq. (12) as:

(14) $$ p\left(z|x,{\theta}_{tr},{\theta}_{te}\right)=\frac{1}{N_x}\sum \limits_{i=1}^{N_x}\int \int \int p\left(z|w,{x}_i^t,{\theta}_{te}^t\right)p\left(w|{\theta}_{tr}^t\right)p\left({\theta}_{tr}^t|{\theta}_{tr}\right)p\left({\theta}_{te}^t|{\theta}_{te}\right)\;d{\theta}_{tr}^t,d{\theta}_{te}^t dw. $$

3.2. Source 2: Uncertainty in the training and testing data

The uncertainty in the training and testing data can be treated in the same way as the uncertainty in the new input vector $ x $ . Similar to Eq. (13), we draw samples from $ p\left({\theta}_{tr}^t|{\theta}_{tr}\right) $ and $ p\left({\theta}_{te}^t|{\theta}_{te}\right) $ and write:

(15) $$ p\left({\theta}_{tr}^t|{\theta}_{tr}\right)=\frac{1}{N_1}\sum \limits_{j1=1}^{N_1}\delta \left({\theta}_{tr}^t-{\theta}_{tr,j1}^t\right) $$

and similar for $ p\left({\theta}_{te}^t|{\theta}_{te}\right) $ . This leads to:

(16) $$ p(z|x,{\theta}_{tr},{\theta}_{te})=\frac{1}{N_x}\sum \limits_{i=1}^{N_x}\frac{1}{N_{X1}}\sum \limits_{j1=1}^{N_{X1}}\frac{1}{N_{X2}}\sum \limits_{j2=1}^{N_{X2}}\int p(z|w,{x}_i^t,{\theta}_{te,j2}^t)p(w|{\theta}_{tr,j1}^t)\hskip0.2em dw, $$

where $ {N}_{X1} $ and $ {N}_{X2} $ are the number of samples drawn from the uncertainty pdf of the input–output pairs in the training and testing datasets, respectively. Note that we have $ {N}_{X1} $ samples for each training input–output pair, and $ {N}_{X2} $ samples for each testing input–output pair.

3.3. Source 3: uncertainty in the weights

3.3.1. Basic derivation

We now turn to the estimation of the weight distribution conditioned on the training data in Eq. (16), namely $ p\left(w|{\theta}_{tr,j1}^t\right) $ . As discussed in the Introduction, this pdf is inherently complex due to the highly nonlinear relationship between the neural network weights and the training data. Consequently, the exact form of this distribution is unknown, and direct sampling from it is not feasible.

Instead, we adopt a pragmatic approach inspired by standard practices in machine learning. In typical training routines, the weights are obtained through a minimization procedure (e.g., stochastic gradient descent), given a fixed neural network architecture and a predefined optimization algorithm. Under these conditions, the final weight vector is determined primarily by two stochastic components: the initialization of the weights and the order in which the training data batches are processed. By varying these components, we can generate multiple plausible realizations of the trained weights. These realizations serve as approximate samples from the unknown distribution $ p\left(w|{\theta}_{tr,j1}^t\right) $ , enabling a Monte Carlo approximation of the pdf. However, as explained below, considerable effort is needed to obtain a practical algorithm.

Let us first introduce the random batch order in our formalism and derive an expression for $ p\left(w|{\theta}_{tr,j1}^t\right) $ that takes this uncertainty into account. By denoting the batch order $ b $ , we can write:

(17) $$ p\left(w|{\theta}_{tr,j1}^t\right)=\int p\left(w,b|{\theta}_{tr,j1}^t\right)\; db=\int p\left(w|b,{\theta}_{tr,j1}^t\right)p\left(b|{\theta}_{tr,j1}^t\right)\; db=\int p\left(w|b,{\theta}_{tr,j1}^t\right)p(b)\; db, $$

where we used the fact that the batch order is independent of the training data. By further employing Bayes Theorem in the form of Eq. (3), and writing the input–output pair of $ {\theta}_{tr,j1}^t $ as $ \left({X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ , we find:

(18) $$ p\left(w|b,{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right)=\frac{l\left({Z}_{tr,j1}^t|w,b,{X}_{tr,j1}^t\right)}{p\left({Z}_{tr,j1}^t|b,{X}_{tr,j1}^t\right)}p\left(w|b,{X}_{tr,j1}^t\right). $$

We can evaluate $ p\left(w|b,{X}_{tr,j1}^t\right)=p(w) $ , since the input data $ {X}_{tr,j1}^t $ and the batch order $ b $ alone have no information on the weights $ w $ . Therefore, the conditional pdf is equal to the prior distribution over the weights.

The likelihood $ l\left({Z}_{tr,j1}^t|w,b,{X}_{tr,j1}^t\right) $ quantifies how well a particular weight vector $ w $ , when applied to the input data $ {X}_{tr,j1}^t $ in the batch order $ b $ , reproduces the corresponding output data $ {Z}_{tr,j1}^t $ of the training dataset. As explained in the Introduction, the denominator, $ p\left({Z}_{tr,j1}^t|b,{X}_{tr,j1}^t\right) $ serves as a normalization constant that ensures the posterior is a valid pdf. Since we are primarily interested in drawing samples from the posterior rather than evaluating it explicitly, this term does not need to be computed in practice.

Inserting Eq. (18) into Eq. (17) we find:

(19) $$ p\left(w|{\theta}_{tr,j1}^t\right)=\int \frac{l\left({Z}_{tr,j1}^t|w,b,{X}_{tr,j1}^t\right)}{p\left({Z}_{tr,j1}^t|b,{X}_{tr,j1}^t\right)}p(w)p(b)\; db $$

In the next section, we focus on determining the likelihood $ l\left({Z}_{tr,j1}^t|w,b,{X}_{tr,j1}^t\right) $ and show how to perform the integral.

3.3.2. The likelihood and the integral

To determine the likelihood, we use the relationship between the training pairs and the weight vectors as specified by the neural network:

(20) $$ {Z}_{tr,j1}^t=f\left({X}_{tr,j1}^t,w,b\right)+\unicode{x025B} \left({X}_{tr,j1}^t\right), $$

where $ \unicode{x025B} \left({X}_{tr,j1}^t\right) $ is the difference between the predicted $ f\left({X}_{tr,j1}^t,w,b\right) $ and the output of training $ {Z}_{tr,j1}^t $ . Because we already incorporated the uncertainty in the input and output training data above in Section 3.2, we can interpret this difference as the uncertainty in the neural network function itself. Rewriting this relation as $ f\left({X}_{tr,j1}^t,w,b\right)={Z}_{tr,j1}^t-\unicode{x025B} \left({X}_{tr,j1}^t\right) $ , and since $ {Z}_{tr,j1}^t $ is given, the uncertainty distribution of $ f\left({X}_{tr,j1}^t,w,b\right) $ is equal to that of $ \unicode{x025B} \left({X}_{tr,j1}^t\right) $ shifted by the constant vector $ {Z}_{tr,j1}^t $ .

As an example, assume the uncertainty in $ \unicode{x025B} \left({X}_{tr,j1}^t\right) $ to be Gaussian, with zero mean and error variance $ {\sigma}_f^2 $ for each element $ m $ of vector $ \unicode{x025B} \left({X}_{tr,j1}^t\right) $ . We then find the likelihood:

(21) $$ l\left({Z}_{tr,j1}^t|w,b,{X}_{tr,j1}^t\right)=A\exp \left[-\frac{1}{2}\sum \limits_{m=1}^{N_{tr}}\frac{{\left({Z}_{m, tr,j1}^t-f\left({X}_{m, tr,j1},w,b\right)\right)}^2}{\sigma_f^2}\right], $$

where $ A $ is a normalization factor, $ {N}_{tr} $ is the total number of the input–output pairs in the training dataset. The term under the exponent can be recognized as a typical machine-learning loss function.

We now turn to the calculation of the integral in Eq. (19). We can draw initial weight vector samples $ {w}_k $ from $ p(w) $ and batch order samples $ {b}_k $ from $ p(b) $ , which leads to:

(22) $$ p\left(w|{\theta}_{tr,j1}^t\right)=\frac{1}{N_w}\sum \limits_{k=1}^{N_w}\frac{l\left({Z}_{tr,j1}^t|{w}_k,{b}_k,{X}_{tr,j1}^t\right)}{p\left({Z}_{tr,j1}^t|{b}_k,{X}_{tr,j1}^t\right)}\delta \left(w-{w}_k\right), $$

where $ {N}_w $ is the total number of drawn weight vector samples. This expression shows that the posterior pdf of the weight vectors is represented by the initial weight vectors $ {w}_k $ , each with a so-called importance weight:

(23) $$ {\alpha}_{j1,k}=\frac{1}{N_w}\frac{l\left({Z}_{tr,j1}^t|{w}_k,{b}_k,{X}_{tr,j1}^t\right)}{p\left({Z}_{tr,j1}^t|{b}_k,{X}_{tr,j1}^t\right)} $$

This implies that any posterior moment of the weight vector $ w $ (e.g., the mean or covariance) must be computed using importance weights $ {\alpha}_{j1,k} $ . However, in practice, the values of $ {\alpha}_{j1,k} $ tend to vary dramatically between samples and, typically, one of them takes a value close to 1, while the others are nearly zero. This is because the likelihood contains an exponential term, as shown in Eq. (21). The sum in the brackets varies considerably with the weight vector $ {w}_k $ , and the exponent of that sum will enhance that variation substantially. As a result, the posterior distribution becomes sharply peaked, and only a single weight vector significantly contributes to the estimate of $ p\left(w|{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ . The effective ensemble size collapses to one, eliminating any representation of uncertainty. This issue is known as filter degeneracy in the particle filter literature (Doucet et al., Reference Doucet, de Freitas and Gordon2001; Van Leeuwen, Reference Van Leeuwen2019; Evensen et al., Reference Evensen, Vossepoel and van Leeuwen2022). It renders the use of Eqs. (22) and (23) impractical, as they fail to maintain a diverse and informative ensemble of weight vectors.

3.3.3. The proposal density

Unlike Eq. (22), where samples of the weight vector $ w $ are drawn from the prior distribution $ p(w) $ , we can instead draw samples from a distribution of trained weights. This approach is mathematically formalized by introducing a proposal density $ q\left(w|{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ which explicitly depends on the training data (Doucet et al., Reference Doucet, de Freitas and Gordon2001; Van Leeuwen, Reference Van Leeuwen2010; Ades and Van Leeuwen, Reference Ades and Van Leeuwen2012). To incorporate this, we multiply and divide the prior pdf $ p(w) $ by the proposal density $ q\left(w|{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ . This procedure will allow us to perform importance sampling using samples drawn from this proposal rather than from the $ p(w) $ . Using this idea in Eq. (19), we find:

(24) $$ {\displaystyle \begin{array}{l}p\left(w|{\theta}_{tr,j1}^t\right)=\int \frac{l\left({Z}_{tr,j1}^t|w,b,{X}_{tr,j1}^t\right)}{p\left({Z}_{tr,j1}^t|b,{X}_{tr,j1}^t\right)}p(w)p(b)\; db\\ {}=\int \frac{l\left({Z}_{tr,j1}^t|w,b,{X}_{tr,j1}^t\right)}{p\left({Z}_{tr,j1}^t|b,{X}_{tr,j1}^t\right)}\frac{p(w)}{q\left(w|b,{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right)}q\left(w|b,{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right)p(b)\; db\\ {}=\frac{1}{N_w}\sum \limits_{k=1}^{N_w}\frac{l\left({Z}_{tr,j1}^t|{w}_{j1,k},{b}_k,{X}_{tr,j1}^t\right)}{p\left({Z}_{tr,j1}^t|{b}_k,{X}_{tr,j1}^t\right)}\frac{p\left({w}_{j1,k}\right)}{q\left({w}_{j1,k}|{b}_k,{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right)}\delta \left(w-{w}_{j1,k}\right),\end{array}} $$

where the weight vectors $ {w}_{j1,k} $ are now drawn from the proposal density $ q\left(w|{b}_k,{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ . The conditioning on the training data means that the weight vectors have been trained on the training data. The details on how they are trained are described in the proposal density.

A standard choice in machine learning is to draw the weight vectors by minimizing a loss function, which is the same as maximizing the likelihood. This minimization is combined with a methodology to avoid overfitting. It is not easy to write down the proposal density in this case. No matter these details, from Eq. (24) we see that $ p\left(w|{\theta}_{tr,j1}^t\right) $ is now represented by a set of weights $ {w}_{j1,k} $ , each with importance weight given as:

(25) $$ {\beta}_{j1,k}=\frac{1}{N_w}\frac{l\left({Z}_{tr,j1}^t|{w}_{j1,k},{b}_k,{X}_{tr,j1}^t\right)}{p\left({Z}_{tr,j1}^t|{b}_k,{X}_{tr,j1}^t\right)}\frac{p\left({w}_{j1,k}\right)}{q\left({w}_{j1,k}|{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right)}. $$

The $ {\beta}_{j1,k} $ consist of three factors that depend on the value of the weight vector, the prior pdf, the proposal density, and the likelihood. The prior $ p(w) $ is typically a uniform distribution around zero, i.e., the probability $ p\left({w}_{j1,k}\right) $ is the same for any drawn weight vector $ {w}_{j1,k} $ . The proposal density $ q\left(w|{b}_k,{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ does depend on the trained weight vector. To see this, denote a prior weight vector sample as $ {u}_k $ and the corresponding trained weight vector as $ {w}_{j1,k} $ . The proposal density is related to the prior density via the standard transformation of variables formula from statistics, as:

(26) $$ q\left({w}_{j1,k}|{b}_k,{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right)={\left|\frac{\partial {w}_{j1,k}}{\partial {u}_k}\right|}_{b_k}p\left({u}_k\right), $$

where $ {\left|\partial {w}_{j1,k}/\partial {u}_k\right|}_{b_k} $ is the Jacobian of this transformation. This Jacobian is a complicated function and depends on $ {u}_k $ and $ {b}_k $ . Although different prior samples $ {u}_k $ will have the same value of $ p\left({u}_k\right) $ , the Jacobians are different for different $ {u}_k $ and $ {b}_k $ , and, consequently, the proposal densities of the trained weights will be different from each other. The likelihood values of the trained weights will also be different because some trained weight vectors will fit the training data better than others. This means that the importance weights $ {\beta}_{j1,k} $ will be different for each weight vector. It is well known from the particle filtering literature that if the training data set is large, the differences in the $ {\beta}_{j1,k} $ will be huge, again leading to one importance weight having a value close to 1, while all others have values very close to zero.

Since Bagging, Deep Ensembles, and MC dropout use the trained weights in this way, they are of little value for providing information on $ p\left(w|{\theta}_{tr,j1}^t\right) $ . However, we can construct a proposal density that does not have this problem, as explained below.

3.3.4. A new proposal density

To overcome the filter degeneracy issue, we design the proposal density $ q\left(w|{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ such that all importance weights $ {\beta}_{j1,k} $ are equal, or at least nearly so. One intuitive approach is to enforce a fixed likelihood value across all sampled weights. Specifically, we define a target likelihood value $ {l}_0 $ and seek intermediate weight vectors $ {\tilde{w}}_{j1,k} $ that satisfy the identity:

(27) $$ l\left({Z}_{tr,j1}^t|{\tilde{w}}_{j1,k},{b}_k,{X}_{tr,j1}^t\right)={l}_0. $$

To ensure that the $ {\beta}_{j1,k} $ are equal, all values $ q\left(\tilde{w}|{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ need to be equal as well. According to Eq. (26) this implies that the Jacobian of the transformation from the untrained weights $ {u}_k $ to the trained weights $ {\tilde{w}}_{j1,k} $ must not depend on $ k $ . However, variations in $ {u}_k $ that change $ {\tilde{w}}_{j1,k} $ such that $ l\left({Z}_{tr,j1}^t|{\tilde{w}}_{j1,k},{b}_k,{X}_{tr,j1}^t\right) $ increases or decreases are not allowed, because the likelihood has to be fixed at $ {l}_0 $ . This means that the matrix $ \left(\partial {\tilde{w}}_{j1,k}/\partial {u}_k\right) $ is not full rank, and its determinant vanishes:

(28) $$ \left|\frac{\partial {\tilde{w}}_{j1,k}}{\partial {u}_k}\right|=0 $$

resulting in a proposal density $ q\left(\tilde{w}|{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ that is not well defined. As a result, while this method is conceptually appealing, it is not mathematically viable for practical use.

To avoid this issue, we define the final weight vectors in the proposal density by adding a small perturbation of size $ \delta $ to each intermediate weight vector $ {\tilde{w}}_{j1,k} $ , as:

(29) $$ {w}_{j1,k}={\tilde{w}}_{j1,k}+\delta {r}_{j1,k}\left({u}_k\right), $$

where $ {r}_{j1,k}\left({u}_k\right) $ is a function of the untrained weight vector $ {u}_k $ . The Jacobian can now be evaluated via a Taylor-series expansion around $ \left|\partial {w}_{j1,k}/\partial {u}_k\right| $ , and, since $ \delta $ is small, we only need to retain up to the first term in the expansion, leading to:

(30) $$ \left|\frac{\partial {w}_{j1,k}}{\partial {u}_k}\right|=\left|\frac{\partial {\tilde{w}}_{j1,k}}{\partial {u}_k}+\delta \frac{\partial {r}_{j1,k}}{\partial {u}_k}\right|\approx \left|\frac{\partial {\tilde{w}}_{j1,k}}{\partial {u}_k}\right|+\delta Tr\left[ Adj\left(\frac{\partial {\tilde{w}}_{j1,k}}{\partial {u}_k}\right)\left(\frac{\partial {r}_{j1,k}}{\partial {u}_k}\right)\right]+O\left({\delta}^2\right), $$

where $ Tr $ denotes the trace, and $ Adj $ denotes the adjugate. If we now choose $ \left(\partial {r}_{j1,k}/\partial {u}_k\right) $ as the pseudoinverse $ \left(\partial {r}_{j1,k}/\partial {u}_k\right)= Adj{\left(\frac{\partial {\tilde{w}}_{j1,k}}{\partial {u}_k}\right)}^{\dagger } $ we find:

(31) $$ \left|\frac{\partial {w}_{j1,k}}{\partial {u}_k}\right|=\delta n+O\left({\delta}^2\right), $$

where $ n=\dim \left({w}_{j1,k}\right) $ , and we used that $ \left|\partial {\tilde{w}}_{j1,k}/\partial {u}_k\right|=0 $ . This shows that, up to order $ \delta $ , the Jacobian is a constant, independent of the weight vector $ {w}_{j1,k} $ . This, then, means that the proposal density $ q\left({w}_{j1,k}|{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ in Eq. (26) has the same value for each weight vector $ {w}_{j1,k} $ .

However, by adding a small perturbation to the trained weight vectors, the likelihood is no longer exactly equal to $ {l}_0 $ . To assess the impact of this perturbation, we perform a Taylor series expansion on the likelihood around the trained weights. This yields:

(32) $$ l\left({Z}_{tr,j1}^t|{w}_{j1,k},{b}_k,{X}_{tr,j1}^t\right)={l}_0+\delta {r}_{j1,k}{\nabla}_{w_{j1,k}}l+O\left({\delta}^2\right), $$

in which $ \delta {r}_{j1,k} $ is the pertubation. Since the weights $ {\tilde{w}}_{j1,k} $ have been obtained through training, the likelihood is already near its maximum, implying that the gradient $ {\nabla}_{w_{j1,k}}l $ is small. By choosing the perturbation magnitude $ \delta $ to be sufficiently small, the first-order term becomes negligible, and the perturbed likelihood remains very close to $ {l}_0 $ . This approximation restores the practical viability of the method.

The target value for the likelihood $ {l}_0 $ is determined as follows. First, we train the neural network with the unperturbed training data, resulting in a weight vector denoted by $ {w}_0 $ . Next, we use the $ {N}_{X1} $ perturbed versions of the training data and pass each of them through the neural network using this fixed weight vector $ {w}_0 $ . This process yields $ {N}_{X1} $ loss function values $ {J}_{j1}\left({Z}_{tr,j1}^t,{X}_{tr,j1}^t,{w}_0\right) $ . We average these as $ {J}_0=1/{N}_{X1}{\sum}_{j1=1}^{N_{X1}}{J}_{j1}\left({Z}_{tr,j1}^t,{X}_{tr,j1}^t,{w}_0\right) $ and take as target likelihood value $ {l}_0=\exp \left(-{J}_0\right) $ . By defining $ {l}_0 $ this way, we obtain a representative likelihood value for trained weight vectors. Consequently, the set of trained weight vectors obtained from the perturbed training datasets can be interpreted as an unbiased sample from the posterior distribution of the weights, $ p\left(w|{X}_{tr,j}^t,{Z}_{tr,j1}^t\right) $ .

With the above construction, where each weight vector $ {w}_{j1,k}\approx {\tilde{w}}_{j1,k} $ , we have ensured that the proposal densities $ q\left({w}_{j1,k}|{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ are identical across all samples. Additionally, the corresponding likelihoods $ l\left({Z}_{tr,j1}^t|{w}_{j1,k},{b}_k,{X}_{tr,j1}^t\right) $ are the same for each weight vector. As a result, all importance weights $ {\beta}_{j1,k} $ are equal, implying that each weight vector $ {w}_{j1,k} $ is equally probable in the posterior distribution $ p\left(w|{X}_{tr,j}^t,{Z}_{tr,j1}^t\right) $ .

With these equally probable samples from $ p\left({w}_{j1,k}|{X}_{tr,j1}^t,{Z}_{tr,j1}^t\right) $ , Eq. (16) can be rewritten as:

(33) $$ p\left(z|x,{\theta}_{tr},{\theta}_{te}\right)=\frac{1}{N_x}\sum \limits_{i=1}^{N_x}\frac{1}{N_{X1}}\sum \limits_{j1=1}^{N_{X1}}\frac{1}{N_{X2}}\sum \limits_{j2=1}^{N_{X2}}\frac{1}{N_w}\sum \limits_{k=1}^{N_w}p\left(z|{w}_{j1,k},{x}_i^t,{X}_{te,j2}^t,{Z}_{te,j2}^t\right). $$

3.4. Uncertainty in the trained neural network

As shown in Eq. (33), the final step is to estimate $ p\left(z|{w}_{j1,k},{x}_i^t,{X}_{te,j2}^t,{Z}_{te,j2}^t\right) $ , which captures the uncertainty due to the neural network’s imperfection. This imperfection is evaluated using the testing data $ \left({X}_{te,j2}^t,{Z}_{te,j2}^t\right) $ for each index $ j2 $ . Since the evaluation is based on the statistical relation described by the input–output pairs of testing data, we rewrite

(34) $$ p\left(z|{w}_{j1,k},{x}_i^t,{X}_{te,j2}^t,{Z}_{te,j2}^t\right)={p}_{\theta_{te,j2,i}}\left(z|{w}_{j1,k},{x}_i^t\right), $$

i.e., the pdf of $ z $ given the input $ {x}_i^t $ and weight vector $ {w}_{j1,k} $ , determined using the input–output pairs in the testing data $ j2 $ . The extra index $ i $ in the testing data denotes that we only use testing data for which the input data is close to $ {x}_i^t $ . Specifically, we calculate the pdf of the output testing data $ {Z}_{te,j2} $ for which the input testing data is within a region around the new input, so $ {X}_{tr,j2}^t\in {D}_{x_i^t} $ . This pdf can be approximated by known functions (e.g., a Gaussian or a Gaussian mixture), by a histogram, by the set of samples $ {Z}_{te,j2,m} $ , or by new samples generated from $ {Z}_{te,j2,m} $ (e.g., via a so-called diffusion process (Sohl-Dickstein et al., Reference Sohl-Dickstein, Weiss, Maheswaranathan and Ganguli2015)).

The size of $ {D}_{x_i^t} $ is governed by the uncertainty in the input testing data. In our experience, setting the radius of $ {D}_{x_i^t} $ between 0.5 and 2 times the uncertainty standard deviation in the input testing data yields robust and consistent output pdf estimates.

3.5. The resulting algorithm and practical considerations

Using Eq. (34) into Eq. (33), we have a complete uncertainty quantification methodology that can be written as:

(35) $$ p\left(z|x,{\theta}_{tr},{\theta}_{te}\right)=\frac{1}{N_x}\sum \limits_{i=1}^{N_x}\frac{1}{N_{X1}}\sum \limits_{j1=1}^{N_{X1}}\frac{1}{N_{X2}}\sum \limits_{j2=1}^{N_{X2}}\frac{1}{N_w}\sum \limits_{k=1}^{N_w}{p}_{\theta_{te,j2,i}}\left(z|{w}_{j1,k},{x}_i^t\right). $$

This equation shows that we average the pdfs $ {p}_{\theta_{te,j2,i}}\left(z|{w}_{j1,k},{x}_i^t\right) $ over the input samples $ i $ , the trained weight vector samples $ k $ , the training data samples $ j1 $ , and the testing data samples $ j2 $ , to find the full uncertainty pdf of the output of the deep learning procedure. The following provides the set of steps needed in this calculation:

1. 1) Initial training: Train a neural network using the original, unperturbed training data. Apply any standard technique to avoid overfitting, but do not use dropout. Store the resulting weight vector as $ {w}_0 $ .

2. 2) Generate perturbed training data: Create an ensemble of $ {N}_{X1} $ training datasets drawn from the uncertainty distributions of the input and output training datasets. This results in $ {N}_{X1} $ perturbed versions of the full training dataset.

3. 3) Determine the target likelihood value $ {l}_0 $ : Use the perturbed input training datasets from step 2 to calculate their loss function values using the corresponding perturbed training output dataset and the fixed weight vector $ {w}_0 $ from step 1. Average the $ {N}_{X1} $ loss function values and store the result in $ {J}_0 $ . Calculate the target likelihood value $ {l}_0=\exp \left(-{J}_0\right) $ .

4. 4) Determine an ensemble of weight vectors: For each perturbed training dataset, train $ {N}_w $ neural networks with different initial untrained weight vectors and batch orders. This leads to $ {N}_{X1}\times {N}_w $ different neural networks. These networks should all be trained such that their loss function value is close to the target loss function value $ {J}_0 $ . This process results in $ {N}_{X1}\times {N}_w $ trained weight vectors.

5. 5) Generate perturbed new input data: Create $ {N}_x $ perturbed new input values drawn from the uncertainty distribution of the new input value. This results in $ {N}_x $ perturbed versions of the new input value.

6. 6) Generate perturbed testing data: Create $ {N}_{X2} $ perturbed testing datasets drawn from the uncertainty distribution of the testing data. This results in $ {N}_{X2} $ perturbed versions of the testing data.

7. 7) Generate samples for the pdf: For each perturbed new input value from step 5, find all the perturbed testing input data from step 6 that are close to it. Pass the selected testing input data through the $ {N}_{X1}\times {N}_w $ neural networks from step 4 and subtract the result from their corresponding testing output data. This results in a large number of output samples for each perturbed new input data.

8. 8) Determine final pdf: Calculate histograms for each perturbed new input data based on the samples calculated in step 7. Average the histograms to find the final uncertainty pdf for output $ z $ .

3.6. Practical tips

1. 1. Ensemble sizes. We tested our algorithm on several test cases, including a real-world application described below. In our experience, the predictive pdf converged for $ {N}_x={N}_{tr}={N}_{te}={N}_w=20 $ in all applications. However, we expect these numbers to be application dependent.

2. 2. Accuracy of likelihood calculations. The likelihood is calculated from the loss function by multiplying it by minus one and then taking the exponential. When a large amount of training data is used, the loss function values will be very large, and taking the exponent of minus that large number can lead to underflow issues. Since we are only interested in relative likelihood values, we can subtract the same constant value from each loss function before taking the exponent. In our algorithm, we used the $ {J}_0 $ value.

3. 3. Simplifications for Gaussian output uncertainties. The above algorithm can be simplified significantly if $ l\left({Z}_{tr,j1}^t|{w}_{j1,k},{b}_k,{X}_{tr,j1}^t\right) $ is assumed to be Gaussian as in Eq. (21), and the uncertainty in the output training data $ {Z}_{tr} $ is also Gaussian distributed. In that case, we do not need to introduce perturbed output training data $ {Z}_{tr}^t $ in Step 2. Instead, we can work directly with $ p\left(w|{X}_{tr,j1}^t,{Z}_{tr}\right) $ . Then the likelihood of the weight vectors becomes (see Appendix A):

(36) $$ l\left({Z}_{tr,j1}^t|w,b,{X}_{tr,j1}^t\right)=A\exp \left[-\frac{1}{2}\sum \limits_{m=1}^{N_{tr}}\frac{{\left({Z}_{tr,m}^t-f\left({X}_{tr,j1,m},w,b\right)\right)}^2}{\sigma_f^2+{\sigma}_{z,m}^2}\right], $$

in which $ {\sigma}_{z,m}^2 $ is the uncertainty variance in the output training data $ {Z}_{tr,m} $ . All steps and equations from Eq. (21) to Eq. (33) remain unchanged, with the only modification being that the perturbed output training vector $ {Z}_{tr,j1}^t $ is replaced everywhere with the original, unperturbed output vector $ {Z}_{tr} $ . This is a situation often encountered in the geosciences.

5.1. Description of the training and testing datasets

In this study, we leverage the training and testing datasets used in Chiu et al. (Reference Chiu, Yang, Van Leeuwen and Feingold2021). The input variables are derived from aircraft cloud probe measurements collected during the Aerosol and Cloud Experiments in the Eastern North Atlantic (ACE-ENA) field campaign, operated by the Atmospheric Radiation Measurement (ARM) user facility. These include cloud water content ( $ {q}_c $ ), cloud droplet number concentration ( $ {N}_c $ ), drizzle water content ( $ {q}_r $ ), and drizzle drop number concentration ( $ {N}_r $ ), with respective uncertainties of 30%, 50%, 30%, and 20% (Glienke and Mei, Reference Glienke and Mei2019, Reference Glienke and Mei2020; Mei et al., Reference Mei, Wang and Comstock2020). The output variable, autoconversion rate, is calculated by applying the in-situ cloud droplet size distribution to the stochastic collection equation formulated as a two-moment bin model. The rate spans a wide dynamic range, from $ 7.7\times {10}^{-24} $ to $ 2.8\times 1{0}^{-6}\hskip0.33em \mathrm{kg}\hskip0.33em {\mathrm{m}}^3{\mathrm{s}}^{-1} $ , capturing diverse stages of cloud evolution and supporting robust training. The uncertainty in the autoconversion rate is approximately 53%, based on the analytical autoconversion equation (6) from Chiu et al. (Reference Chiu, Yang, Van Leeuwen and Feingold2021). The full dataset, consisting of approximately 10,612,000 samples, is partitioned into training, validation, and testing subsets, with a 60–20-20% data split. To reduce dynamic range and satisfy the assumption of additive Gaussian noise, a base-10 logarithmic transformation is applied to all input and output variables.

5.2. Training procedure

Building on the neural network architecture used by Chiu et al. (Reference Chiu, Yang, Van Leeuwen and Feingold2021), we implement a feedforward network comprising six fully connected hidden layers, each containing 16 nodes. Following the procedure outlined in Section 3.5, we first trained the network using the unperturbed training dataset, with randomly initialized weights and a shuffled batch order. To mitigate overfitting, we applied the early stopping criterion proposed by Prechelt (Reference Prechelt1998), halting training when the validation loss no longer decreased with additional epochs. Given the high variability of the validation loss across epochs, we smoothed the loss curve using a moving average over a 50-epoch window. We then identified the epoch at which the smoothed loss stabilized, and used the corresponding model weights as the final trained weights, $ {w}_0 $ .

Following the initial training with unperturbed data, we then trained separate neural networks for each of the 20 perturbed training datasets $ {X}_{tr,j1},{Z}_{tr,j1} $ , using 20 different sets of perturbed initial weights and batch orders. This yielded $ {N}_w=20 $ models per perturbed dataset. Training was stopped when the loss function for the training data reached $ {J}_0 $ . As a result, we obtained a total of 400 neural networks (20 networks for each of the 20 perturbed datasets) to quantify uncertainty in both the model weights and the training data.

In practice, exact convergence to $ {J}_0 $ is complicated by the fact that many machine learning frameworks report loss only at the end of each epoch. To address this, we identified the two epochs where the loss values bracketed $ {J}_0 $ , and then applied a bisection method to approximate a loss value within 0.01 of $ {J}_0 $ . This threshold ensures that the likelihood values are sufficiently close to $ {l}_0=\exp \left(-{J}_0\right) $ , leading to approximately uniform importance weights across the ensemble. As shown in Figure 3, most importance weights are close to the expected value of $ 1/400=0.0025 $ , which corresponds to a perfect equal-weight ensemble. A small number of networks did not converge precisely, which is expected due to the high dimensionality of the optimization space. For these networks, training was stopped at slightly less precise loss values, resulting in importance weights deviating by no more than 6%, which is acceptable for our purposes.

Figure 3. The normalized likelihood values (i.e., the importance weights) for the $ 400 $ neural networks.

5.3. Results

To evaluate our methodology on new input data, we selected 50 distinct input vectors spanning a wide range of values. Each input vector was perturbed 100 times using random noise sampled from the input uncertainty distributions described in Section 4.1. This resulted in 40,000 output samples per input vector ( $ {N}_x=100 $ , $ {N}_w=20 $ , $ {N}_{X1}=20 $ ). For each input, we computed the predictive uncertainty distribution using the testing data, as described in Eq. (29), and averaged the resulting histograms.

As mentioned in Section 3.4, the region ( $ {D}_{x_i^t} $ ) around each perturbed new input $ {x}_i^t $ , must be defined to compute the local empirical distribution from the test data. We used all output testing data that had corresponding input testing data close to $ {x}_i^t $ within $ 0.2 $ (in $ \log 10 $ space), for each of the four variables in the input vector. We experimented with thresholds ranging from 0.1 to 0.5 and found that this choice had no significant visual impact on the resulting distributions.

The results for all $ 50 $ input vectors are displayed in Appendix B. To highlight the key findings, Figure 4 shows the total uncertainty pdfs as a function of the $ \log 10 $ autoconversion rate for four representative input vectors. In each panel, the blue curve is the total uncertainty pdf from Eq. (30), and the red bar is the point prediction from the baseline neural network (i.e., using $ {w}_0 $ ) without uncertainty quantification. The black bars show the predictions from Bagging, obtained by creating an ensemble of 20 networks with different initial weights and batch orders.

Figure 4. Examples of total uncertainty pdf in the output autoconversion rate for four input vectors. The blue curves are the total uncertainty pdfs; the red bar is the autoconversion rate calculated by the baseline neural network ( $ {w}_0 $ ) without uncertainty quantification; black bars are the Bagging output samples. Note the wide variety of shapes of the blue uncertainty pdfs and the small spread in the black bars, demonstrating the inadequacy of the Bagging approach to represent uncertainty.

Figures 4(a) and (b) show that the predictive pdfs span roughly an order of magnitude in autoconversion rate, which is not surprising given the relatively large measurement uncertainty in the input variables (ranging between 20% and 50%). Although the pdf in Figure 4(a) can be approximated by a Gaussian distribution, a non-negligible skewness is observed in Figure 4(b). The pdfs in Figure 4(c) and (d) are broader, showing large uncertainty and bimodality.

To disentangle the contributions of different uncertainty sources, we first isolated the effect of neural network weight uncertainty. This was done by disabling all other uncertainty sources in Eq. (30) and computing outputs using only the ensemble of weight vectors. The resulting distributions were very narrow and closely resembled the Bagging ensemble outputs (black bars), indicating that weight uncertainty is not a major contributor.

Next, we assessed the combined impact of input uncertainty and model uncertainty from the testing data. This involved evaluating the corresponding terms in Eq. (30) while ignoring the other uncertainties. These two sources contributed similar amounts of variability and together accounted for most of the spread in the predictive distribution, as shown by the black curves in Figure 5, which peak at similar values to the full pdfs (blue curves). However, the combined input and model uncertainty distributions were typically unimodal and symmetric (in $ {\log}_{10} $ space), failing to reproduce the skewness or bimodality observed in the full predictive distributions. This suggests that the observed skewness and bimodality are primarily caused by uncertainty in the training and testing datasets.

Figure 5. Same as Figure 4, showing full uncertainty pdfs (blue curves) and the contribution from input uncertainty and model uncertainty (black curves).

To further investigate the skewness and bimodality in the predictive pdfs, we analyzed all 50 cases in Appendix B. We found that positive skewness occurred predominantly at high autoconversion values (above $ 10-12 $ ), while strong negative skewness was associated with intermediate values ( $ 14-17 $ ). Bimodal distributions were mainly observed at low autoconversion values (around $ -18 $ ). Attempts to relate these features to specific patterns in the 4D input space were inconclusive, and even 3D visualizations of input subsets failed to reveal consistent structure.

Returning to Figure 4, we observe the limitations of Bagging. The resulting ensemble predictions (black bars) exhibit little spread and fail to capture features such as bimodality, leading to significant underestimation of uncertainty. In fact, the shortcomings of Bagging extend beyond this. As noted in Section 3.3 (following Eq. 21), ensemble methods like Bagging, Deep Ensembles, and MC Dropout implicitly assume equal likelihood for all network realizations. However, in practice, the importance weights of these networks, derived from their likelihood values, vary exponentially.

To illustrate this further, we normalized all likelihoods by the largest value. This resulted in one network receiving an importance weight of 1, while the next highest weight was $ \exp \left(-824\right) $ , effectively zero. All remaining networks contributed negligibly. Thus, only the highest-likelihood network meaningfully affects the predictive distribution, rendering the ensemble effectively a single-model prediction. This outcome stands in stark contrast to our method, which produces near-uniform importance weights (see Figure 3). Accordingly, if standard ensemble-based methods were applied in a statistically consistent manner, they would collapse to a single model and fail to provide meaningful uncertainty quantification.