# Evidential Deep Learning: Enhancing Predictive Uncertainty Estimation for
Earth System Science Applications

###### Abstract

Robust quantification of predictive uncertainty is critical for understanding
the factors the drive different outcomes in weather and climate predictions.
Ensembles provide predictive uncertainty estimates and can be decomposed in
physically relevant ways, but both physics-based and machine learning
ensembles are computationally expensive and can underestimate the probability
of extreme outcomes. Parametric probabilistic deep learning models can produce
an uncertainty estimate with a single model by predicting the parameters of a
known probability distribution but do not account for epistemic, or
uncertainty from the model, when doing so. Evidential deep learning, a newer
technique that extends parametric deep learning to higher-order distributions,
can account for both aleatoric and epistemic uncertainty with a single model.
This study compares the uncertainty metrics derived from evidential neural
networks to those obtained from calibrated ensembles. Through example
applications of classification of winter precipitation type and regression of
surface layer fluxes, we show evidential deep learning models attaining
predictive accuracy rivaling standard methods, while robustly quantifying both
sources of uncertainty uncertainty. We evaluate the uncertainty in terms of
how well the predictions are calibrated and how well the uncertainty
discriminates between correct and incorrect predictions. Analyses of
uncertainty in the context of the prediction inputs reveal sensitivities to
underlying meteorological processes, facilitating interpretation of the
models. The conceptual simplicity, interpretability, and computational
efficiency of evidential neural networks make them highly extensible, offering
a promising approach for reliable and practical uncertainty quantification in
Earth system science modeling. In order to encourage broader adoption of
evidential deep learning in Earth System Science, we have developed a new
Python package, MILES-GUESS (https://github.com/ai2es/miles-guess), that
enables users to train and evaluate both evidential deep learning and ensemble
deep learning models.

##  1 Introduction

Uncertainty is an inherent aspect of any prediction (Abdar et al., 2021), yet
it is often overlooked or not communicated effectively alongside the
prediction itself (Gneiting et al., 2007). The ability to provide well-
calibrated and robust predictive uncertainty estimates can be highly valuable,
allowing users to understand the reliability of predictions and make more
informed decisions based on them (e.g., Nadav-Greenberg and Joslyn, 2009;
Kendall and Gal, 2017). For model developers, accurate predictive uncertainty
estimates can help identify challenging cases and determine when a model may
be operating outside its training domain (Kendall and Gal, 2017; Karpatne,
2018). Additionally, by connecting uncertainty estimates with other analysis
tools, researchers can gain insights into the input sensitivities that
influence uncertainty levels, allowing a better understanding of the factors
that drive these uncertainties (Herman and Schumacher, 2018; Liu et al.,
2020). In the realm of machine learning, various approaches have emerged for
quantifying total predictive uncertainty. These include well-established
methods such as bagging (Breiman, 1996), Gaussian processes (Rasmussen and
Williams, 2006), quantile regression (Koenker, 2005), as well as newer
conformal methods (Romano et al., 2019; Stankeviciute et al., 2021;
Angelopoulos and Bates, 2021). However, it is essential to recognize that
existing methods often come with inherent limitations. Conventional
techniques, such as bagging, may encounter difficulties in capturing the full
spectrum of uncertainty, particularly when dealing with intricate, multimodal
distributions. Moreover, many established methods often lack the ability to
decompose predictive uncertainty into its underlying components, hindering a
deeper understanding of uncertainties. Additionally, handling custom
probability distributions can be challenging for some of these methods,
limiting their adaptability to specific problem domains.

Traditionally, uncertainty quantification (UQ) within the Earth System Science
community has been pursued through the development and enhancement of physics-
based numerical models, which generate probabilistic forecasts using ensembles
of deterministic forecasts that vary in initial conditions, boundary
conditions, or model specifications (Leith, 1974). One notable advantage of
these methods is their strong foundation in the true physics of
atmospheric/oceanic motion. However, deterministic numerical model ensembles
come with considerable computational costs and often lack proper uncertainty
calibration (Vannitsem et al., 2018). To combat the computational expense and
lack of calibration, UQ has been attempted through statistical linking
functions. Two prominent techniques are ensemble model output statistics
(Gneiting et al., 2005), in which a parametric distribution is prescribed and
fit, and Bayesian model averaging (Raftery et al., 2005), where the UQ takes
the form of a weighted mixture distribution.

The use of modern Machine Learning (ML) for Earth system UQ has been the focus
of much recent research (McGovern et al., 2017; Haynes et al., 2023),
especially within the forecast post-processing community (Haupt et al., 2021;
Schulz and Lerch, 2022; Vannitsem et al., 2020). Popular techniques include
parametric fitting (Ghazvinian et al., 2021; Rasp and Lerch, 2018; Chapman et
al., 2022; Guillaumin and Zanna, 2021; Barnes and Barnes, 2021; Foster et al.,
2021; Delaunay and Christensen, 2022; Gordon and Barnes, 2022), quantile-based
probabilities transformed to full predictive distributions (Scheuerer et al.,
2020), or creating direct approximations of the quantile function via
regression based on Bernstein polynomials (Bremnes, 2020).

Many users of uncertainty rely on a single metric, such as probability or
ensemble spread, but decomposing uncertainty into different components
provides valuable insights into its sources and nature (Kendall and Gal,
2017), and can help validate how well a model captures the different
uncertainty sources. In statistics, the law of total variance decomposes
uncertainty into aleatoric, arising from inherent data randomness, and
epistemic, arising from model limitations and insufficient training data
(Kendall and Gal, 2017). High aleatoric uncertainty indicates the lack of a
clear relationship between the model inputs and the target and can only be
reduced with the addition of more informative input variables (Herman and
Schumacher, 2018). High epistemic uncertainty can be reduced by accumulating
more data in sparse areas of the input space or by reducing the complexity or
flexibility of the model. Such distinction aids in assessing model learning
capacity, generalization, and guiding hyperparameter optimization (Karpatne,
2018). Advancing this field requires techniques that efficiently decompose
uncertainty while achieving general predictive reliability.

While parametric probabilistic machine learning models, such as those that
predict the parameters of a categorical or Gaussian probability distribution,
can express aleatoric uncertainty (Nix and Weigend, 1994), they do not account
for epistemic uncertainty (Amini et al., 2020). The predicted variance only
accounts for data variance, not model variance. Ensembles of deterministic ML
models (Lakshminarayanan et al., 2016) and sampling methods, such as Monte
Carlo dropout (Srivastava et al., 2014; Gal and Ghahramani, 2016) approximate
epistemic uncertainty by deriving their spread from model perturbations, but
if their predictions are single labels or values, then they are not accounting
for spread in the data. On the other hand, Bayesian neural networks (Neal,
2012), which treat every weight as a random variable, can accurately estimate
both aleatoric and epistemic uncertainty but are computationally demanding and
challenging to implement for complex architectures. Ensembles of parametric
probabilistic models can also be used to approximate aleatoric and epistemic
uncertainty (Delaunay and Christensen, 2022) but only with a sufficiently
large ensemble size (Shaker and Hüllermeier, 2020).

The recent rise of evidential neural networks (Sensoy et al., 2018; Amini et
al., 2020; Ulmer et al., 2023) offers a promising solution that strikes a
balance between efficiency and accuracy while providing an effective approach
to estimate both sources of uncertainty. Evidential neural networks use a
single deterministic neural network while modifying the prediction task to
estimate the parameters of a higher-order evidential distribution, which draws
relevance from Bayesian data analysis principles (Gelman et al., 2013). This
distribution treats the parameters of the target distribution as random
variables and models them with an assumed prior distribution (Sensoy et al.,
2018; Amini et al., 2020). For multinomial (categorical) and Gaussian
distributions, analytical formulations of conjugate prior distributions, such
as the Dirichlet and Normal Inverse Gamma (Sensoy et al., 2018; Amini et al.,
2020), enable the construction of exact loss functions for neural network
training. These loss functions consist of a negative log-likelihood component
to maximize the fit to the data and a regularizer term to minimize evidence
allocated to incorrect predictions and inflate the conditional uncertainty
(Sensoy et al., 2018).

In this work, we introduce the concept of evidential deep learning to the
Earth system science community. Evidential deep learning represents a
relatively recent machine learning technique known for its ability to provide
predictive uncertainty estimates with practical advantages, as discussed by
Sensoy et al. (2020). Our primary focus is on the application of evidential
neural networks within the weather and climate domain, where the accurate
estimation of uncertainty holds significance in decision-making, as
highlighted in previous studies (Shepherd, 2009; Bauer et al., 2015). To
assess the utility of evidential neural networks, we employ them in both
classification and regression tasks. In the classification domain, we showcase
an evidential neural network trained to predict winter precipitation type
based on simulated atmospheric soundings. In regression tasks, our evaluation
centers on assessing the model’s performance in estimating surface energy
fluxes using observed meteorological variables (McCandless et al., 2022).

In addition to practical applications, we have defined three key objectives:
(1) To quantitatively assess and compare the predictive capabilities of
evidential versus deterministic neural networks using essential metrics such
as Brier Skill Score (BSS) and Root Mean Square Error (RMSE). (2) To evaluate
the calibration of predicted uncertainties derived from evidential models and
ensembles through various analysis techniques. (3) To explore calibration
tuning approaches tailored to evidential networks, including parameter
adjustments such as a loss regularization weight for regression and the fine-
tuning of the dropout rate for Monte Carlo (MC) ensembles. These objectives
aim to provide a balanced understanding of the potential of evidential
architectures in generating uncertainty estimates while maintaining accuracy,
and emphasize the significance of uncertainties in assessing the limitations
of ML models.

##  2 Methods

###  2.1 Example Problems of Interest

Figure 1: (a-b) Example of a precipitation type sounding with temperature,
dewpoint, and wind profiles along with (c-d) different predictive uncertainty
representations.

In the realm of supervised learning, a neural network is employed to predict
the most likely label or outcome based on a given set of input predictor
variables. Figure 1 provides an illustration of model inputs and outputs for a
classification problem, based on one of the example problems explored in this
paper. The inputs are shown in Figure 1(a) and (b) and include temperature
($T$; ∘C), dewpoint temperature ($T_{\text{dew}}$; ∘C), and zonal ($U$; m s-1)
and meridional ($V$; m s-1) wind, at equally spaced heights in the atmosphere.
The output is a winter precipitation-type label derived from a set of
precipitation-type probabilities: $\mathbf{p}$ =
$(p_{rain},p_{snow},p_{sleet},p_{frzrain})$ (Figure 1(c)).

The contrast between the two modeling approaches investigated here (ensemble
versus evidential) when predicting the ‘sleet’ label is illustrated in Figure
1(c). In the case of the evidential model, the predictions are represented by
a smooth curve, offering a continuous and comprehensive quantification of
possible outcomes. Conversely, an ensemble approach such as k-fold cross
validation or Monte Carlo sampling, provides discrete point estimates, as seen
in the figure, which are obtained by collecting multiple model predictions.
It’s important to note that these point estimates require a full ensemble to
approximate the continuous probabilistic insights provided directly by the
evidential model. This distinction highlights the effectiveness of the
evidential model in capturing the full spectrum of uncertainty for the full
set of labels, as opposed to the more discrete and ensemble-based
representation of uncertainty.

The data used for winter precipitation-type modeling consists of two main
sources. The target data source is the NOAA National Severe Storms Laboratory
Meteorological Phenomena Identification Near the Ground, or mPING (Elmore et
al., 2014). mPING is a crowd-sourced, GPS-location-linked, weather reporting
project that allows anyone to report precipitation type and other weather
hazards anonymously across the globe. For this project, we have selected only
snow, rain, freezing rain, and sleet/ice pellet reports and have narrowed the
data window to be between 1 January 2015 and 30 June 2022. Data after 1 July
2020 were withheld for testing and training and validation data was split
randomly after being grouped by day. The value of these reports come from
being able to gather information from any populated location rather than at
fixed weather stations, but the reports contain biases and uncertainties
related to users’ abilities to discern different precipitation types from each
other.

The input data are simulated vertical atmospheric profiles from the NOAA Rapid
Refresh (RAP) numerical weather prediction model analysis (Benjamin et al.,
2016). The RAP has a 13 km horizontal resolution on a Lambert Conformal grid
over North America. RAP temperature, relative humidity, wind, and pressure are
interpolated from pressure levels ranging every 25 hPa from 1000 to 100 hPa to
height levels every 250 m from 0 to 5 km above ground level. Dewpoint is
derived from relative humidity, temperature and pressure using the metpy
package (May et al., 2022). Surface temperature, dewpoint, and wind are
extracted from the RAP 2 m temperature and dewpoint and 10 m wind components.
mPING observations are grouped by hour within the bounds of each RAP grid
cell. The most common precipitation type is used as the target for the machine
learning models.

Within the matched mPING and RAP data, a small percentage of cases contained
reports that were physically inconsistent with the RAP profile. These cases
were found after initial neural networks trained on the data reported nonzero
probabilities of frozen precipitation even when near-surface temperatures were
well above 0 ∘C. Further investigation revealed many profiles labeled as sleet
and freezing rain but lacking a profile where the wet bulb temperature goes
above 0 ∘C before returning to below freezing. Some reports were labeled with
a frozen type but did not have any freezing layers near the surface or near
surface moisture. When matched with radar, some of these reports were not
colocated with any precipitation on radar. To address these inconsistencies,
we filtered examples where there was a frozen precipitation type (snow,
freezing rain, or sleet) reported but the surface wet-bulb temperature
exceeded 3 ∘C or if rain was reported and the surface wet-bulb temperatures
was less than -3 ∘C. We tested other wet-bulb temperature thresholds up to 10
∘C and found that 3 ∘C filtered most of the inconsistent reports without
removing too many plausible reports.

Figure 2(a) illustrates a multi-layer perceptron (MLP) neural network
architecture tasked with predicting precipitation-type class probabilities
that we will investigate here. In a similar vein, as illustrated in Figure
2(b), we employ another MLP architecture, but this time for a regression
problem. The model predicts energy fluxes from Earth’s surface given near-
surface atmospheric conditions (McCandless et al., 2022; Muñoz-Esparza et al.,
2022b). The ultimate goal of this problem is to train a machine learning
parameterization of the surface layer from observations that could be used to
replace existing parameterizations in weather and climate models. In this
example, the model takes as input four quantities: wind speed (m s-1) at 10 m
above ground level (AGL), the change in potential temperature (K s-1) between
the surface and 10 m AGL, the Bulk Richardson number (unitless) valid at 10 m
AGL, and the change in water vapor mixing ratio (g kg-1 m-1) between the
surface and 10 m AGL. The illustrated architecture predicts up to three
outputs: friction velocity (m s-1), kinematic sensible heat flux (kg m s-1),
and kinematic latent heat flux (kg m s-1).

The dataset used in the surface layer energy flux problem (SL) provides long-
term observations from the Royal Netherlands Meteorological Institute (KNMI)
Cabauw Experimental Site for Atmospheric Research located in Cabauw,
Netherlands (Bosveld et al., 2020). The site has been utilized for validating
surface-layer parameterizations and land-surface models since 1972. The area
surrounding the site consists of flat agricultural land with minimal elevation
changes and low brush coverage.

We used the Cabauw data gather from 2003-2017 for training and testing the ML
models. The dataset includes various meteorological measurements noted earlier
that are used as inputs to the model as well as model training targets
(McCandless et al., 2022; Muñoz-Esparza et al., 2022a). The outgoing longwave
radiation, converted to skin temperature, is derived from a device located
approximately 200 meters from the flux tower. The data were grouped by day and
then split randomly into training, validation, and testing sets. The same
random seed was used for splitting to ensure that all ML models used the same
testing set. McCandless et al. (2022) recently showed that the Cabauw dataset
could be used to train effectively several machine learning models including a
MLP like the one illustrated in Figure 2(a) with good performance.

Figure 2: (a) Deterministic and evidential MLP architectures for predicting
class probabilities in the precipitation-type categorical dataset. (b)
Architectures for predicting parameters in the surface-layer regression
dataset, including Gaussian ($\mu$, $\sigma^{2}$) and Normal-Inverse Gaussian
($\gamma$, $\nu$, $\alpha$, $\beta$) distributions. In both architecture
diagrams, the neural network is represented by $\text{f}_{\textbf{w}}$ where w
are trainable parameters.

###  2.2 Classification Challenges

Once trained, the ML models are used to make predictions for a series of
events. We define an event as the space and time extent over which a single
prediction is valid, such as from 1800 to 1900 UTC on 2 February 2023 in the
area within 13 km of the center of Boulder, Colorado. In classification
problems, ML models output a number associated with each of $K$ known outcomes
for an event, given a vector of inputs $\mathbf{x}$. If the numbers each range
between 0 and 1 and sum to 1, then the numbers behave as probabilities and
thus represent the parameters of a categorical distribution describing a
constant probability of outcome occurrence during the event.

$p(\hat{y}=y_{i}|\mathbf{x})=p_{i}(\mathbf{x})$ (1)

Although predictions cover the volume of an event, the observations that
verify the predictions generally occur only at a single point in space and
time within that event. In the winter precipitation type case, the observed
outcome is a series of crowd-sourced reports occurring within a given area and
timeframe that may vary from the most likely outcome due to a mix of subscale
meteorological and societal factors. To model the potential variability in
observed outcomes for an event, we can use the multinomial distribution, which
generalizes the categorical distribution to the outcomes of $n$ repeated
observations under the same event conditions. Mathematically, we express the
multinomial probability for a set of $n$ observations with a fixed number for
each outcome $\mathbf{y}=\\{y_{i},...,y_{K}\\}$ where $\sum_{i=1}^{K}y_{i}=n$,
given a fixed set of probabilities for each outcome
$\mathbf{p}=\\{p_{i},...,p_{K}\\}$ as

$p(\mathbf{y}|n,\mathbf{p},K)=\frac{{n!}}{{\prod_{i=1}^{K}y_{i}!}}\prod_{i=1}^{K}p_{i}^{y_{i}}.$
(2)

To train the ML model to estimate $\mathbf{p}$ for a given event, we construct
a training dataset $\mathcal{D}$ = $\\{x_{j},y_{j}\\}_{j=1}^{N}$ that we then
input into the ML model repeatedly in order to adjust the parameters (weights)
$\mathbf{w}$ of the ML model to maximize their likelihood given the training
data. This statistical optimization process, called maximum likelihood
estimation, involves the model ingesting the inputs and predicting outputs y
to minimize the negative log-likelihood loss function

$\mathcal{L}(\mathbf{w}|\mathcal{D})=-\sum_{i=1}^{K}y_{i}\ln(\hat{y}_{i}(x)).$
(3)

Maximum likelihood estimation assumes that there we are seeking one fixed set
of $\mathbf{w}$ that best explains our data due to its origins in frequentist
statistics (MacKay, 2003). Thus, these predicted probabilities only account
for aleatoric uncertainty from the variation in the outcome for each input.

###  2.3 Model Ensembles

Additional model weights could also plausibly explain the data and these
plausible model weights could be found by repeating the maximum likelihood
estimation process with different random initializations or by resampling the
training data. Both approaches can produce ensembles of probabilistic ML
models (Lakshminarayanan et al., 2016; Gal and Ghahramani, 2016; Dietterich,
2000). Another approach is the Bayesian neural network (Neal, 2012), which
treats the weights of the neural network as random variables that are each
represented with a parametric probability distribution. Both approaches allow
for the estimation of both aleatoric and epistemic uncertainty, extending the
model’s ability to capture and express uncertainty in predictions. Ensembles
are straightforward to train but often require a large number of members or
many repeated samples to produce robust uncertainty estimates. Bayesian neural
networks require more weights per model compared with standard MLPs and are
less likely to converge to results that both perform well deterministically
and produce robust uncertainty estimates.

To generate ensembles, we will consider three methods: cross-validation
splitting, deep ensembling (Lakshminarayanan et al., 2016), and Monte Carlo
dropout (Gal and Ghahramani, 2016). Cross-validation is a standard approach
for training and evaluating the performance of neural models. It involves
dividing the dataset into multiple subsets or folds, using a fold as a
validation set while the remaining folds form the training set. The ensemble
model is then both trained and evaluated multiple times, with each fold
serving as the validation set for a different model. Ensembles of k-fold=20
were used.

Deep ensembling (Lakshminarayanan et al., 2016) extends the concept of
ensemble methods to deep learning architectures. Instead of training a single
neural network, multiple instances of the same network are independently
trained. Each network explores different regions of the parameter space due to
random initialization or training dynamics. To create the ensemble, we employ
different weight parameterizations of the neural networks using fixed seeds,
ensuring that all models in the ensemble are trained on the same data split.
Similarly, ensembles of size 20 were used.

Finally, Monte Carlo dropout (Gal and Ghahramani, 2016) involves randomly
deactivating neurons during model training using a dropout rate. During
inference, the model is used to obtain $N_{MC}$ predictions while keeping
dropout enabled. We used $N_{MC}=100$ in all cases to create the ensembles of
predictions. Note that all of these ensemble generation methods are applicable
to both categorical and regression problems.

###  2.4 Evidential Classification

The recent evidential deep learning approaches find a compromise between
ensemble methods and Bayesian neural networks by using fixed weights for the
neural network but predict instead the parameters for a higher-order posterior
distribution that describes the space of possible lower-order distributions
that could only be partially sampled by ensembles (Sensoy et al., 2018; Amini
et al., 2020; Meinert and Lavin, 2021; Oh and Shin, 2022; Meinert et al.,
2023). In Bayesian inference, prior information is combined with the observed
data using Bayes’ theorem to update our beliefs and obtain the posterior
distribution. Recall that Bayes’ theorem states that the posterior probability
of a parameter given the data is proportional to the product of the prior
probability and the likelihood of the data given the parameter: Posterior
$\sim$ Prior $\cdot$ Likelihood. By including a prior distribution in our
inference process, we can define the space of possible models and thus include
the epistemic uncertainty. However, in order to make the inference problem
tractable, we must pick a prior distribution that can be analytically related
to our posterior distribution.

Given these constraints, the best choice of prior for our
multinomial/categorical target distribution is the Dirichlet distribution
(Murphy, 2007; Hoffman et al., 2013). It is described by hyperparameters
$\mathbf{\alpha}$= $(\alpha_{1},\alpha_{2},...,\alpha_{N})$, the probability
density function has the following form

$f(\mathbf{p}|\bm{\alpha})=\begin{cases}\frac{1}{B(\bm{\alpha})}\prod_{i=1}^{K}p_{i}^{\alpha_{i}-1}&\text{for
}\textbf{p}\in
S_{K},B(\alpha)=\frac{\prod_{i=1}^{K}\Gamma(\alpha_{i})}{\Gamma(\sum_{i=1}^{K}a_{i})}\\\
0&\text{otherwise},\end{cases}$ (4)

where $S_{K}$ is the unit simplex defined as $S_{K}$ =
$\\{\mathbf{p}|\sum_{i=1}^{K}p_{i}=1~{}\text{and}~{}0\leq p_{1},...,p_{K}\leq
1\\}$, and $B(\mathbf{\alpha})$ is the K-dimensional multinomial beta function
(Kotz et al., 2004). Note the similar form to Equation 2. The expected
probability for the $k$th outcome is the mean of the Dirichlet posterior
distribution

$E[p_{k}]=\hat{p}_{k}=\frac{\alpha_{k}}{S}.$ (5)

The Dirichlet is chosen because it is the conjugate prior to the multinomial,
meaning that the posterior when the Dirichlet prior is combined with a
multinomial likelihood is also a Dirichlet. This drastically simplifies the
updating process during Bayesian inference because each posterior $\alpha_{i}$
can be computed by summing the prior $\alpha_{i}$ with each observed or
predicted $y_{i}$. In this context, each $\alpha_{i}$ is a pseudocount of how
many outcomes of a certain type are expected to occur during the event.

Sensoy et al. (2018) also motivated the use of the Dirichlet posterior
distribution through the Dempster-Shafer Theory of Evidence (DST), which is an
extension of Bayesian statistics (Dempster, 1968) to decision-making with
uncertain evidence. Evidence in the DST classification context is represented
as a non-negative number $e_{i}$ with higher values indicating stronger
evidence for a particular outcome. The uninformed evidence prior is that each
outcome has evidence of 1, so summing the prior evidence with how much
evidence is observed for each outcome from data or a predictive model, we
receive a posterior $\alpha_{k}$ = $e_{k}+1$ and the total amount of evidence
$S=\sum_{i=1}^{K}(e_{i}+1)$. The $\alpha$ values can then be plugged into a
Dirichlet distribution to derive probabilities for each outcome.

DST is extended for decision-making in classification problems through
Subjective Logic (Jsang, 2018). Subjective Logic, as a formal framework for
modeling uncertainty and subjective beliefs, utilizes belief masses to
quantify belief strength and enables nuanced representation of subjective
confidence levels. Within this framework, a belief mass $b_{i}$ is the amount
of evidence assigned to a particular outcome, and uncertainty mass is the
amount of evidence not allocated to any outcome,

$\displaystyle b_{k}$ $\displaystyle=\frac{e_{k}}{S},$ (6) $\displaystyle u$
$\displaystyle=\frac{K}{S}.$ (7)

$K$ is usually defined as the number of outcomes, assuming no prior evidence
for any outcome. Belief masses look similar to probabilities but only have to
sum to 1 when $u$ is included:

$u+\sum_{i=1}^{K}b_{k}=1.$ (8)

The color bar in Figure 1(d) illustrates the belief masses plus $u$ for a
precipitation-type evidential model.

By incorporating prior knowledge through base rates, Subjective Logic
contextualizes beliefs within a broader information framework. In the
integration of DST and Subjective Logic, the Dirichlet distribution provides a
natural framework for representing belief masses, which are used to quantify
the strength of subjective beliefs. The Dirichlet distribution also enables
flexible modeling of uncertainty by capturing the covariance among different
belief masses.

Figure 2(a)(ii) shows an MLP architecture employing the ‘evidential’
categorical model. The parametric neural network parameterized by
$\mathbf{w}$, and given an input sample $x_{i}$, now predicts the evidence
vector, $e_{i}$, represented by $f(x_{i}|\mathbf{w})$, hence the name
assignment ‘evidential MLP’. Compared to a deterministic classifier, the only
architectural difference is the output activation function, which is taken to
be ReLU following Sensoy et al. (2018) rather than softmax, to filter non-
negative evidence. Accordingly, the Dirichlet distribution corresponding with
this evidence has parameters $\alpha_{i}$ = $f(x_{i}|\mathbf{w})+1$. The
predicted expected probabilities for the classes are then computed as
$\alpha_{i}/S_{i}$.

Sensoy et al. (2018) explored using one loss based on Type-II Maximum
Likelihood, and Malinin and Gales (2018) focused on computing the Bayesian
risk with respect to the class predictor. Here we only present one of the risk
losses as the authors noted the superior performance. Specifically, the
parametric neural network employs the following loss during training

$\displaystyle L_{i}(\textbf{w})$
$\displaystyle=\int||\mathbf{y}_{i}-\mathbf{p}_{i}||^{2}\frac{1}{B(\mathbf{\alpha})}\prod_{j=1}^{K}p_{ij}^{\alpha_{ij}-1}d\mathbf{p}_{i}$
(9)
$\displaystyle=\sum_{j=1}^{K}\left(y_{ij}-\hat{p}_{ij}\right)^{2}+\frac{\hat{p}_{ij}(1-\hat{p}_{ij})}{S_{i}+1}.$
(10)

Sensoy et al. (2018) prove three propositions to explain how Eq. 9 optimizes
the models: (1) The proposed loss function in the neural network prioritizes
data fit by minimizing the variance of the generated Dirichlet experiment for
each sample in the training set. (2) It decreases the prediction error when
new evidence is added to the Dirichlet parameters and increases when evidence
is removed, ensuring effective learning from new information. (3)
Additionally, the loss function exhibits learned loss attenuation, reducing
prediction error by removing evidence from the largest Dirichlet parameter for
a sample, except when the class label is the correct one.

By combining these propositions, it is concluded that the neural network with
the proposed loss function generates more evidence for the correct class
labels and avoids misclassification by removing excessive misleading evidence.
The loss function also aims to decrease the variance of predictions on the
training set by increasing evidence, but only if it leads to a better data
fit. The loss over a batch of training samples is computed by summing the loss
for each individual sample. The model learns to generate evidence for specific
class labels based on patterns in the data, and it adjusts the generated
evidence based on backpropagation when counter-examples are observed. However,
limited counter-examples may result in an increased overall loss when reducing
the evidence magnitude, leading to some misleading evidence for incorrect
labels.

To address this limitation, a Kullback-Leibler (KL) divergence term is
incorporated into the loss function

$\mathcal{L}=\sum_{i=1}^{N}\mathcal{L}_{i}(\textbf{w})+\nu_{t}\sum_{i=1}^{N}KL\left[D(\mathbf{p}_{i}|\widetilde{\alpha_{i}})||D(\mathbf{p}_{i}|\langle
1,...,1\rangle)\right],$ (11)

where $\nu_{t}$ is an annealing coefficient that depends on the epoch number
$t$, $D(\mathbf{p}_{i}|\langle 1,...,1\rangle)$ represents a uniform Dirichlet
distribution, and
$\widetilde{\mathbf{\alpha}_{i}}=\mathbf{y}_{i}+(1-\mathbf{y}_{i})\odot\mathbf{\alpha}_{i}$
represents the Dirichlet parameter after the misleading evidence has been
removed for sample $i$. This term regularizes the predictive distribution and
penalizes divergences from a state of uncertainty (“I don’t know”) that do not
contribute to data fit. The loss function with this regularization term is
used to guide the model’s behavior and achieve the desired total evidence
shrinkage to zero for samples that cannot be correctly classified. The KL
divergence term can be reduced to

$\displaystyle
KL\left[D(\mathbf{p}_{i})|\widetilde{\mathbf{\alpha}_{i}}||D(\mathbf{p}_{i})|\mathbf{1}\right]$
$\displaystyle=\log\left(\frac{\Gamma\left(\sum_{l=1}^{K}\widetilde{\alpha}_{il}\right)}{\Gamma(K)\prod_{l=1}^{K}\Gamma(\widetilde{\alpha}_{il})}\right)$
(12)
$\displaystyle+\sum_{l=1}^{K}\left(\widetilde{\alpha}_{il}-\mathbf{1}\right)\left[\psi(\widetilde{\alpha}_{il}-\mathbf{1})-\psi\left(\sum_{j=1}^{K}\widetilde{\alpha}_{ij}\right)\right],$

where $\psi$ represents the digamma function, $\Gamma$ is the Gamma function,
and $\mathbf{1}$ is a vector of ones. By incrementally increasing the impact
of the KL divergence in the loss function via the annealing coefficient, the
neural network is allowed the opportunity to explore different parameter
configurations. This strategy helps prevent early convergence to the uniform
distribution for misclassified samples, allowing the network to potentially
classify them correctly in subsequent training iterations.

###  2.5 Evidential Regression

Thus far, our discussion has exclusively centered around classification
problems. However, when it comes to predicting continuous targets, such as in
regression problems, the concept of the “I don’t know” outcome is not
explicitly defined because the range of possible outcomes is unbounded.
Subjective logic is primarily tailored to handle uncertainty in categorical or
discrete scenarios, where belief masses and degrees of uncertainty can be
assigned to various hypotheses or propositions. Amini et al. (2020) approached
the regression problem similarly to Sensoy et al. (2018) but adjusted the
formulation to handle continuous outcomes by changing both the target
distribution and the loss function.

In the context of regression tasks, we work with a dataset $\mathcal{D}$ =
$\\{x_{i},y_{i}\\}_{i=1}^{N}$ where the targets $y$ are assumed to be
independently and identically distributed (i.i.d.) according to a Gaussian
distribution characterized by mean $\mu$ and variance $\sigma^{2}$ (Gelman et
al., 1995). The deterministic model shown in Figure 2(b)(i) predicts just
$\hat{y}$ and is trained using an appropriate loss function such as the mean-
squared error. Figure 2(b)(ii) shows a parametric model for predicting
parameters $\mu$ and variance $\sigma^{2}$. Here, as with the multinomial
distribution for categorical problems, if maximum likelihood is used to derive
an optimization loss, the likelihood parameters $\mu$ and $\sigma^{2}$ are
treated as fixed and deterministic, hindering the estimation of model
uncertainty. Amini et al. (2020)’s evidential regression assumes the mean
$\mu$ is drawn from a Gaussian distribution, while the variance $\sigma^{2}$
follows an Inverse-Gamma distribution

$\displaystyle(y_{1},...,y_{N})\sim\mathcal{N}(\mu,\sigma^{2})$ (13)
$\displaystyle\mu\sim\mathcal{N}(\gamma,\sigma^{2}\nu^{-1})$ (14)
$\displaystyle\sigma^{2}\sim\Gamma^{-1}(\alpha,\beta),$ (15)

where $\mathbf{m}$ = $(\gamma,\nu,\alpha,\beta)$, $\gamma\in\mathbb{R}$,
$\nu>0$, $\alpha>1$, and $\beta>0$. This formulation results in a higher-order
distribution referred to as the “evidential distribution”, denoted as
$p(\theta|\mathbf{m})$, which can be represented by a Normal-Inverse-Gamma
distribution

$\displaystyle p(\mu,\sigma^{2}|\gamma,\nu,\alpha,\beta)$
$\displaystyle=\frac{{\sqrt{\nu}}}{{\sqrt{2\pi\sigma^{2}}}}\frac{{\beta^{\alpha}}}{{\Gamma(\alpha)}}\left(\frac{1}{{\sigma^{2}}}\right)^{\alpha+1}\exp\left(-\frac{{2\beta+\nu(\mu-\gamma)^{2}}}{{2\sigma^{2}}}\right)$
(16)
$\displaystyle=\text{St}\left(y_{i};\gamma,\frac{\beta(1+\nu)}{\nu\alpha},2\alpha\right),$
(17)

where $\text{St}\left(y;\mu_{St},\sigma^{2}_{St},\nu_{St}\right)$ is the
Student-t distribution evaluated at $y$ with scale $\sigma^{2}_{St}$ and
$\nu_{St}$ degrees of freedom. For multi-task models, there are four
parameters for each model target, as is shown in Figure 2(b)(iii). By learning
the parameters $\mathbf{m}$ through training, evidential deep learning models
define full distributions over the likelihood parameters $(\mu,\sigma^{2})$,
allowing for a comprehensive representation of uncertainty in the model’s
predictions. The loss used to train a parametric neural network for predicting
$\mathbf{m}$, is computed by taking the negative logarithm of Equation 16

$\mathcal{L}_{i}^{NLL}(\mathbf{w})=\frac{1}{2}\log\left(\frac{\pi}{\nu}\right)-\alpha\log(\Omega)+\left(\alpha+\frac{1}{2}\right)\log\left((y_{i}-\gamma)^{2}\nu+\Omega\right)+\log\left(\frac{\Gamma(\alpha)}{\Gamma(\alpha+\frac{1}{2})}\right),$
(18)

where $\Omega=2\beta(1+\nu)$. Following the approach employed by Sensoy et al.
(2018), Amini et al. (2020) included an additional regularizer term to
suppress evidence (or raise the uncertainty) in support of incorrect
predictions

$\mathcal{L}_{i}^{R}(\mathbf{w})=|y_{i}-\gamma|(2\nu+\alpha),$ (19)

where the first term represents the model error while the second term is
proportional to the total evidence accumulated by the learned posterior. The
total loss used during training is finally

$\mathcal{L}_{i}(\mathbf{w})=\mathcal{L}_{i}^{NLL}(\mathbf{w})+\lambda\mathcal{L}_{i}^{R}(\mathbf{w}),$
(20)

where the parameter $\lambda$ is selected to best calibrate the model. As we
show in the results section (and as others have shown), if $\lambda$ is too
small, the model tends to over-fit the data, while overly large values of
$\lambda$ lead to uncertainty over-inflation (Soleimany et al., 2021).

###  2.6 Law of Total Variance and Entropy

The law of total variance (LoTV; Casella and Berger, 2002) states that for two
random variables $X$ and $Y$ on the same probability space, the variance of
variable $Y$ may be decomposed as

$\text{Var}(Y)=\mathbb{E}[\text{Var}(Y|X)]+\text{Var}(\mathbb{E}[Y|X]),$ (21)

where $\text{Var}(Y)$ represents the total variance of the random variable
$Y$, $\mathbb{E}[\text{Var}(Y|X)]$ denotes the expected value of the
conditional variance of $Y$ given $X$, and $\text{Var}(\mathbb{E}[Y|X])$
represents the variance of the conditional means of $Y$ given $X$. The first
term is often referred to as the “aleatoric” uncertainty or “uncertainty in
the data,” while the latter represents the “epistemic” uncertainty or
“uncertainty in the model’s predictions.”

For the Dirichlet distribution, these quantities can be computed as

$\displaystyle\underbrace{\mathbb{E}[\text{Var}(Y|p_{j})]}_{\text{Aleatoric}}$
$\displaystyle=\mathbb{E}[p_{j}(1-p_{j})]$ (22)
$\displaystyle=\frac{\alpha_{j}}{S_{\alpha}}-\left(\frac{\alpha_{j}}{S_{\alpha}}\right)^{2}-\frac{\frac{\alpha_{j}}{S_{\alpha}}\left(1-\frac{\alpha_{j}}{S_{\alpha}}\right)}{S_{\alpha}+1}$
(23)
$\displaystyle\underbrace{\text{Var}(\mathbb{E}[Y|p_{j}])}_{\text{Epistemic}}$
$\displaystyle=\text{Var}(p_{j})$ (24)
$\displaystyle=\frac{\frac{\alpha_{j}}{S_{\alpha}}\left(1-\frac{\alpha_{j}}{S_{\alpha}}\right)}{S_{\alpha}+1}.$
(25)

The epistemic uncertainty computed with the LoTV and the quantity $u$ from DST
do not have the same form.

For the Normal-Inverse Gamma distribution, we have

$\displaystyle\underbrace{\mathbb{E}[\mu]=\gamma}_{\text{Prediction}}$ (26)
$\displaystyle\underbrace{\mathbb{E}[\sigma^{2}]=\frac{\beta}{\alpha-1}}_{\text{Aleatoric}}$
(27)
$\displaystyle\underbrace{\text{Var}[\mu]=\frac{\beta}{\nu(\alpha-1)}}_{\text{Epistemic}}$
(28)

Note the relationship between the two uncertainty quantities: the epistemic
uncertainty is the aleatoric uncertainty divided by parameter $\nu$, which is
interpreted as a “virtual observation count” that controls how strongly the
prior influences the posterior to the evidence acquired from the data. As
$\nu$ increases, the prior dominates, so the data must be very persuasive to
move the posterior away from the prior prediction $\gamma$. With lower $\nu$,
the data readily overrides the prior, allowing the posterior to diverge from
$\gamma$. So $\nu$ modulates the impact of the data versus the prior, rather
than literally counting observations.

The LoTV is not limited to stochastic models; it can also be applied to
ensembles of probabilistic models. For instance, in the case of a Gaussian
parametric model that predicts ($\mu$, $\sigma^{2}$), an ensemble can be
created using various approaches discussed below. Conversely, in a categorical
problem, the variance can be estimated directly from the predicted
probabilities and the true labels, as shown in Equation S5. Given an input $x$
to the ensemble of trained models, the output consists of a list of predicted
means and variances. By using Equation 21, we can compute the aleatoric and
epistemic components.

In ensembles involving categorical models, additional statistical quantities
can be computed alongside Equation 21. One such quantity is entropy, denoted
as $H$, which serves as a measure of uncertainty encompassing both the data
distribution and the model’s predictions. The formulation for entropy, as
given by Equation 30, is expressed as follows:

$\displaystyle p_{j}=\frac{1}{N}\sum_{i=1}^{N}p_{ij}$ (29) $\displaystyle
H=-\sum_{j=1}^{K}p_{j}\log(p_{j}),$ (30)

In this equation, $p_{j}$ represents the ensemble average over $N$ ensemble
members. The range of values for entropy spans from 0 to $\log(K)$, where 0
signifies perfect certainty (indicating that only one class has a probability
of 1), and $\log(K)$ signifies maximum uncertainty (implying that all classes
have equal probabilities).

###  2.7 Evaluation Metrics

To quantitatively assess the performance and uncertainty of our models trained
on categorical (precipitation type) and regression (SL) datasets, we adopt the
evaluation methods outlined by Haynes et al. (2023). This assessment framework
centers around four key objectives described below. First, the accuracy of
mean predictions generated by our models is prioritized and relies on metrics
such as the Brier Skill Score (BSS) and Root Mean Square Error (RMSE). These
metrics, visualized through the attributes diagram, help to gauge how well
each model produces an accurate and calibrated mean forecast.

For classification problems, the Brier score (BS) measures the mean squared
difference between the observed outcomes ($o_{ij}$) and the predicted
probabilities ($p_{ij}$) for each class $i$ and outcome $j$, given by

$\text{BS}=\frac{1}{NK}\sum_{i=1}^{N}\sum_{j=1}^{K}(o_{ij}-p_{ij})^{2},$ (31)

where $N$ is the total number of events and $K$ is the total number of
possible outcomes (class labels). A Brier score of zero indicates perfect
accuracy, while a score of 1 represents the worst possible performance. We
also utilize the BSS, which is computed as $BSS$ =
$1-\frac{BS_{\text{forecast}}}{BS_{\text{climatology}}}$, where
$BS_{\text{forecast}}$ is the Brier score of the predicted forecast, while
$BS_{\text{climatology}}$ is the Brier score of the climatological forecast,
which is the mean squared difference between the observed frequency of the
event and the predicted probability of the event based on climatology (the
long-term historical frequency of the event). The BSS ranges from negative
infinity to 1. A BSS of 1 indicates a perfect forecast (no error), and a
negative BSS suggests that the forecast is worse than the climatology. A BSS
of 0 means the forecast is as skillful as climatology.

In the regression problem, we employ the root mean square error (RMSE) as the
performance metric, defined as

$\text{RMSE}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}},$ (32)

where $y_{i}$ represents the ground truth labels and $\hat{y}_{i}$ denotes the
model predictions.

Secondly, we quantify the calibration of uncertainty estimates, ensuring that
predicted uncertainties closely align with observed outcomes. This alignment
is assessed using the Probability Integral Transform (PIT) plot, which
assesses if the model produces unbiased and well-calibrated distributions. The
PIT represents the quantile of the predicted distribution at which the
observed value occurs. It is calculated by approximating the cumulative
distribution function (CDF) of the predicted distribution and evaluating it
with the observed value. The PIT histogram displays the distribution of PIT
values over multiple data samples, aiming for a uniform histogram in a
perfectly calibrated model. The deviation score for the PIT (Nipen and Stull,
2011) can be computed as

$\text{PITD}=\sqrt{\frac{1}{M}\sum_{i=1}^{N_{M}}\left(\frac{N_{M}}{N}-\frac{1}{M}\right)^{2}}.$
(33)

In this formula, $M$ represents the number of bins used for the PIT histogram,
$N_{M}$ is the count of samples in each bin, and $N$ denotes the total number
of samples. The PITD value ranges between 0 and 1, with 0 indicating the best
possible calibration. However, the score is sensitive to the number of bins
used for the evaluation. To control for this sensitivity, we evaluate the
calibration performance in comparison to the worst-case scenario, we define
the PITD Skill Score as follows

$\text{PITD Skill Score}=1-\frac{\text{PITD}}{\text{PITD}_{\text{worst}}},$
(34)

where $\text{PITD}_{\text{worst}}$ represents the worst possible PITD score,
which assumes that all the forecasts end up in one of the outermost bins. The
PITD Skill Score quantifies how well the probabilistic forecasts align with
the observed probabilities, with higher scores (1 being the best) indicating
better calibration relative to the worst-case scenario. Both a uniform PITD
and a 1-1 relationship in the spread-skill diagram indicate that a model is
calibrated according to uncertainty, and not a uniform PITD alone (Hamill,
2001).

Thirdly, we investigate discrimination performance, particularly concerning
the confidence levels in a models’ predictions. The discard fraction plot aids
in understanding how the model’s ability to differentiate between various
outcomes improves as prediction confidence grows, involves sorting the data by
uncertainty and progressively removing data points with higher uncertainties
to assess any improvement in model performance (Dürr et al., 2020; Barnes and
Barnes, 2021). This test provides insights into the correlation between
prediction error and uncertainty, highlighting whether excluding instances
with high uncertainty enhances model performance. In operational scenarios,
the discard test is particularly valuable for setting model spread thresholds
based on acceptable error levels.

Lastly, in the regression problem we calculate the dependency of the RMSE on
the predicted spread as an alternative to traditional error metrics, offering
a means to normalize model performance against a baseline. This approach,
depicted in the skill-score diagram, helps to provide further insights into
the relationship between RMSE and forecast spread (Delle Monache et al.,
2013).

##  3 Results

###  3.1 precipitation-type

Figure 3: Attribution curves for the (a) MC ensemble using the deterministic
model and (b) the evidential model. The columns show the result for each
precipitation-type. In each subpanel, the diagonal, horizontal, and vertical
dashed lines indicate the 1-to-1 line, no-resolution line, and climatology
line. The blue-shaded area indicates skill relative to climatology. Red-shaded
rectangles illustrate the reliability of each model in predicting each class.

Figure 3 presents a comparison between a deterministic MLP with Monte Carlo
Dropout in (a) and an evidential MLP in (b) using attributes curves. Both
models demonstrate similar performance for all precipitation-types in terms of
task prediction. Specifically, when predicting rain, both models provide
predictions that are approximately calibrated with high probabilities (greater
than 0.9), which account for the majority of rain predictions. However, at low
probabilities, the models tend to make over-confident and incorrect
predictions, deviating from the ideal 1-1 line. This indicates a lack of
proper calibration in these regions, where the models exhibit excessive
certainty despite making erroneous predictions. Overall the model performs the
best on rain, then snow, freezing rain, and finally sleet. These results are
summarized by the confusion matrix plotted in Figure S1 for both deterministic
and evidential MLPs, which again show very similar performance between the two
modeling approaches.

The attributes diagram shows that for snow and freezing rain, both models
exhibit largely calibrated forecasts, with most points closely aligning with
the diagonal. However, concerning sleet predictions, an unexpected downward
trend at high probabilities is observed in the output of both models. The
predictions with the highest predicted probability of sleet often have a
different reported outcome. This issue may be attributed to inconsistent
labels in the mPING dataset that were not captured by the wet-bulb temperature
filter. Some instances of wet snow may be mistakenly labeled as sleet. While
not directly addressing data quality issues here, this finding provides an
opportunity to explore the effectiveness of computed and predicted
uncertainties when dealing with noisy data.

Figure 4: 2D histograms illustrate the relationship between the top-1
aleatoric and top-1 epistemic uncertainties for (a) the MC-ensemble
deterministic model and (b) the evidential model. In (c) and (d), the total
summed quantities for these uncertainties are compared.

Figures 4(a) and 4(b) compare aleatoric and epistemic uncertainties for the
true class (referred to as ‘top-1’) in an Monte Carlo ensemble of
deterministic MLPs and an evidential MLP model, respectively. Figures 4(c) and
4(d) display the summed variances across the four precipitation-types.
Throughout the analysis, the predicted aleatoric uncertainty consistently
surpasses the predicted epistemic uncertainty (see x/y tick magnitudes),
indicating higher variability in the RAP and mPING training datasets compared
to the model fit.

The histograms in Figures 4(b) and (d) are shaped using Equations 23 and 25,
respectively. The ensemble MLP’s histogram shape resembles that of the
evidential model but appears more spread out, suggesting possible stochastic
approximation in the ensemble model. On the other hand, the results in Figure
4(b) and (d) are computed using derived expressions. Comparing the ensemble
MLP to the evidential MLP, the top-1 versus the summed distributions show
greater similarity in the ensemble model, while the evidential model exhibits
a more complex relationship, in contrast to the generally monotonic behavior
observed in the summed distribution (see below).

Figure 5: 2D histograms compare uncertainty estimates from the LoTV and DST.
Figures (a) and (b) plot the top-1 predicted aleatoric and epistemic
uncertainties, respectively, against the evidential uncertainty $u$ from DST.
Figures (c) and (d) show the total aleatoric and total epistemic
uncertainties, summed over all classes, versus $u$. The dashed line indicates
a 1-to-1 relationship in each panel.

Figure 5 compares the aleatoric and epistemic contributions to the total
uncertainty according to the LoTV, and DST estimates of uncertainty. The LoTV
quantifies variance for each precipitation-type, while DST produces a single
uncertainty value $u$ across all classes per input $x$. Although both aim to
measure model prediction uncertainty, they differ in focus. The LoTV targets
model limitations through epistemic variance. In contrast, DST reflects
uncertainties in the evidence itself and resulting belief in hypotheses,
rather than just model shortcomings. Using both provides comprehensive
insights into uncertainty sources.

Figures 5(a) and (c) show the aleatoric uncertainty exceeds DST $u$ for most
events. This result implies that there is strong evidence for at least one of
the predicted precipitation types but that in many cases one cannot be
unambiguously determined due to the inherent noisiness of the relationship
between the RAP profile and the observed precipitation type. As a numerical
weather model, RAP has unavoidable errors from numerical discretization
approximations and sub-grid parameterizations. Its coarse spatiotemporal
resolution also likely misses some atmospheric and terrain variations, further
inflating aleatoric uncertainty. Additionally, high collinearity among RAP
inputs like temperature causes aleatoric amplification. Errors in one variable
propagate across correlated inputs, boosting uncertainty beyond what RAP
restrictions alone would produce. This collinearity helps explain aleatoric
uncertainty surpassing DST $u$.

When aleatoric uncertainty dominates, potential mPING label errors may not
substantially impact epistemic or DST uncertainties. In this regime, model
uncertainties stem more from RAP’s variability rather than its own
limitations. However, when DST $u$ is relatively high and exceeds aleatoric
uncertainty (right of the dashed line), mPING label errors could raise DST
uncertainty through inconsistent evidence. Here, truth label uncertainties
have greater influence on belief and plausibility functions, increasing DST
uncertainty compared to aleatoric.

As $u$ approaches unity, the model becomes more uncertain about any outcome
being more likely than any other. Evidence counts even out across classes,
observed as flatter summed aleatoric contributions in Figures 5(a) and (c).
With higher DST uncertainty, the model concentrates less high counts in one
class, slightly reducing top-1 aleatoric uncertainty.

In Figures 5(b) and 5(d), the model’s epistemic uncertainties are smaller than
corresponding DST $u$, as DST allows belief masses to exceed individual class
probabilities when considering evidence from multiple sources. In contrast,
epistemic uncertainty reflects inherent model limitations and is typically
narrower. The monotonic relationship between epistemic and DST $u$ suggests
that higher uncertainty about precipitation evidence is a result of the
model’s limitations. As the DST $u$ value increases, the uncertainty in the
model’s predictions grows, indicating an increasing lack of confidence in its
ability to make precise determinations. The linear relationship between small
epistemic and DST $u$ values indicates proportional scaling of lower
uncertainties. Larger DST $u$ values arise from its ability to combine
evidence from different sources, encompassing more uncertainties than
epistemic uncertainty, which solely focuses on model uncertainties. Thus,
accounting for multiple uncertainty sources in the model’s predictions is
crucial.

Figure 6: 2D histograms illustrate the relationship of Brier score versus
standard deviation computed using the total (summed) aleatoric, epistemic, and
total uncertainties, as well as DST/entropy estimates. In (a) the MC-ensemble
deterministic results are shown, while (b) shows the evidential model. The
dashed (red) line shows 1-to-1. The color bar shows the count and applies to
all subpanels.

Continuing the analysis, an essential step involves comparing the Brier score
with various uncertainty estimates to assess how well these uncertainties
align with true outcomes and the model’s predictive performance. Figure 6
displays 2D histograms that relate the computed Brier score to various
uncertainties. The histograms are presented in four columns: aleatoric,
epistemic, total, and entropy in (a), and DST $u$ in (b). Figure 6(a) shows
results from a deterministic MLP using Monte Carlo sampling, while 6(b) is
from an evidential MLP.

Both models generally exhibit aleatoric uncertainty values that surpass the
corresponding Brier scores, displaying an approximately sigmoidal correlation
between skill and uncertainty. This relationship is more pronounced in the
evidential model. Similarly, the epistemic uncertainty is usually outweighed
by the Brier score for both models, with a notable drop in the Brier score at
the highest epistemic uncertainties, particularly in the ensemble model. The
evidential model generally predicts smaller or similar Brier scores for the
same epistemic value.

The total uncertainty is comparable between the two approaches, with aleatoric
uncertainty exerting a more significant influence than epistemic uncertainty.
The ensemble model’s overall uncertainty shape resembles that of the
evidential model, though with more noise and a broader distribution along
total aleatoric uncertainty.

In the ensemble model, the computed entropy closely mirrors the total
uncertainty, primarily driven by aleatoric contributions. On the other hand,
the shape of the DST $u$ histogram aligns more with epistemic uncertainty than
with the total uncertainty. Both relationships exhibit a somewhat sigmoidal
pattern rather than a linear one. When assuming a linear relationship, entropy
becomes over-dispersed initially, while the evidential model remains under-
dispersed. However, both eventually demonstrate over-dispersed behavior as $u$
increases, indicating the uncertainty surpasses the skill score. Additionally,
both models show a slight over-dispersion for aleatoric uncertainty and under-
dispersion for epistemic uncertainty.

In the Supplementary Information, we analyze the relationship between dropout
rate and uncertainty calibration for the ensemble of deterministic MLPs.
Figures S2(a) and S2(b) demonstrate the extent to which the dropout rate
influences calibration. Overall, the dropout can be used to fine-tune the
calibration, but the role is limited. Additionally, for the evidential model,
the loss weight on the KL divergence term (Equation 11) can be used to some
extent to adjust calibration for the precipitation-type data set (results not
shown). Notably, these hyperparameters’ dependence is not very strong, and
significant adjustments are generally unnecessary for effective uncertainty
calibration. In fact, the uncertainties for the evidential model are
predominantly well-calibrated through model optimization on the F1 score
alone, without extensive tuning of dropout or loss weighting (see Supplemental
Information for training and optimization details).

Conversely, calibrating the deterministic model requires the creation of an
ensemble of sufficient size to estimate uncertainty distributions, which may
grow with the number of classes. As a result, the evidential model offers
distinct advantages and is considered the preferred model due to its ability
to achieve very similar performance to the ensembles of a deterministic model
and well-calibrated uncertainties without the need for complex ensemble
techniques or turning other hyper-parameters Schreck and Gagne (2021).

Figure 7: Discard-test diagrams show the fraction of data points removed from
the test set versus the BSS computed on the remaining subset. Dashed lines
illustrate the evidential model while solid lines show the MC-ensemble
deterministic model. Points above the horizontal dashed (black) line have
skill higher relative to climatology.

Figure 7 illustrates the discard-test for the ensemble and evidential models
across aleatoric, epistemic, total, and entropy/DST uncertainties. In this
approach, a threshold is set for the model’s predicted uncertainty. During
deployment, any predictions above this uncertainty threshold are discarded or
considered “out of confidence.” Rather than providing potentially misleading
or unreliable results in uncertain situations, the model simply abstains from
making predictions. This conservative strategy ensures the model does not
output predictions when it is uncertain or lacks sufficient evidence.

The BSS is used as the skill-score metric for the discard test. The figure
shows that for rain and snow, the BSS consistently improves as more uncertain
data is withheld from evaluation. Freezing rain demonstrates a similar trend
but with a much lower overall BSS. In fact for the evidential model, only
about 15% of test data with the lowest freezing rain uncertainties had
positive skill scores, reflecting the model’s poorer performance in this
precipitation-type compared to the relatively easier rain and snow categories.
With the ensemble MLP, the result was even worse. Sleet appears as an outlier,
with the slope of the curve flipping between positive and negative across
uncertainty levels, signaling a lack of calibration. As noted previously, some
of the sleet performance and calibration issues likely stem from RAP profiles
with physically inconsistent mPING labels.

Figure 8: 17 December 2016 precipitation-type predictions and their
uncertainties from the evidential model are visualized for the central US.
Panel (a) displays the most likely precipitation-type prediction, while panels
(b-d) showcase the total aleatoric, DST, and epistemic uncertainty,
respectively. The red, purple, and navy contours indicate the 0 C isotherm at
the surface, 2 km AGL, and the 0-5 km AGL maximum, respectively.

Figure 8 presents a case study of a severe winter weather event that impacted
a significant portion of the country on December 17th, 2016. The evaluation of
the evidential model was conducted across the entire RAP domain, with the
model predicting the probability of each precipitation type conditioned on
precipitation occurring. In Figure 8(a), the model’s predictions for freezing
rain, sleet, snow, and rain are depicted. Freezing rain and sleet are
predicted over a vast region, stretching from western Kansas eastward and then
southward into South Carolina. Snow is primarily predicted for points north,
while rain is predicted for points south, corresponding to the mixed
precipitation types.

Figure 8(b) and (c) illustrate the patterns of the predicted aleatoric and
epistemic uncertainties, respectively. The highest uncertainties are observed
in regions where one precipitation-type transitions into another and along the
freezing line. Conversely, predictions of rain and snow are generally
accompanied by smaller uncertainties.

Of particular interest is the close correlation between the epistemic
uncertainty and the DST $u$ quantity. As mentioned earlier, the DST
uncertainty usually has a larger magnitude compared to the epistemic
uncertainty. However, Figure 8(d) clearly shows $u$ spiking to 1 in a
relatively small region over northern Georgia stretching up to southern
Kentucky, indicating that the model is completely uncertain about those
predictions. This situation also corresponds to the highest value of epistemic
uncertainty on the grid, showing that both uncertainties align during moments
of high uncertainty when the model is unsure about the precipitation-type.
Both DST $u$ and total epistemic uncertainty are highest along the rain to
freezing rain transition line and lower when transitioning to sleet and snow.
This higher epistemic uncertainty may be due both to near-surface temperature
uncertainty as well as uncertainty from mPING reporters being able to
discriminate between freezing rain and rain consistently.

Figure 9: (a) Composite soundings are presented lightest-to-darkest using the
median of the 10th, 20th, 50th, and 90th percentile least uncertain
predictions generated by the evidential model. The total uncertainty is
utilized for the MC-ensemble MLP (i) and the evidential MLP (ii), while the
DST uncertainty is employed for the evidential MLP in (iii). (b)(i) showcases
the breakdown of the top 10% least uncertain composites (depicting the median)
based on class and uncertainty. Panels (b)(ii-iii) illustrate the densities of
the 90th percentile least uncertain predictions, using DST $u$ and LoTV total
uncertainty, respectively.

Lastly, Figure 9 shows composite soundings for various models and uncertainty
percentiles. The panels in Figure 9(a)(i-iii) demonstrate that lower
uncertainty predictions correspond with more prototypical soundings for the
corresponding prediction type. DST uncertainty shows the best performance
since the ground level temperature of the composite sleet soundings remains
below the freezing line. The other figures, for the MC ensemble and total
uncertainty, have an above-freezing composite at the ground level.

The discrepancies between the different types of uncertainty for the
evidential neural network are exhibited by Figure 9(b)(i). In this figure,
only sleet has varying composite soundings between the different types of
uncertainty (for the 90th percentile of least uncertainty predictions). This
occurs because even if a prediction is in the 90th percentile of epistemic or
aleatoric uncertainty, the other uncertainty type could still be high. As
noted earlier, the aleatoric uncertainty typically dominates the epistemic
uncertainty, and thus is why the median sleet sounding of the total
uncertainty in Figure 9(a)(iii) is close to the median aleatoric/epistemic
sounding in Figure 9(b)(i). For the density plots in Figures 9(b)(ii-iii),
which show sleet, while the range of the soundings is about equal, the
distribution using DST uncertainty is skewed more negative than using total
uncertainty. Overall the DST uncertainty seems to identify more physically
reasonable soundings for sleet.

###  3.2 Surface layer flux problem

Figure 10: The PITD score as a function of the evidential loss coefficient
($\lambda$, Equation 19) for (a) one multi-task model and (b) three single-
task models.

In this section, we present a similar evaluation of the evidential regression
model trained on the SL data. The key point of comparison, detailed below, is
that while the evidential regression model provides a solution for modeling
uncertainties like the categorical approach, its calibration is much more
sensitive to the weight $\lambda$ multiplying the KL term in Equation 20
during training. Furthermore, we observe the regression model infrequently
producing unrealistic uncertainties, even when calibrated, potentially posing
operational challenges. In addition to the evidential regression results, we
also present results from ensembles created using cross-validation, deep
ensembling, and Monte Carlo dropout for comparison.

We initially address the issue of model calibration by examining Figure 10(a),
which displays the PITD score as a function of the KL weight ($\lambda$) for a
three-task model that predicts fraction velocity, latent heat, and sensible
heat. The dataset’s characteristics prevent the model from achieving the
maximum PITD score at the same $\lambda$ value for all three tasks, and none
of the tasks approached a perfect skill score of one. To further explore this
relationship, we analyze Figure 10(b), which depicts the same dependency
computed using three single-task models. Effective utilization of the
evidential regression model on the SL dataset requires three single-task
models to determine the optimal $\lambda$ values (illustrated in the figure),
and furthermore, the values differ relative to the three-task model shown in
Figure 10(a). Recall that a flat PIT histogram (PITD equal to zero) does not
necessarily guarantee calibration (Chapman et al., 2022; Haynes et al., 2023;
Hamill, 2001).

Figure 11: Attribution curves for the (a) data ensemble using the parametric
Gaussian model and (b) the evidential model. The red line shows the
reliability curve. The columns show the result for each model task. Similar to
Figure 3, In each subpanel, the diagonal, horizontal, and vertical dashed
lines indicate the 1-to-1 line, no-resolution line, and climatology line. The
blue-shaded area indicates skill relative to climatology.

With the calibration weights identified, the three models were trained and
compared against an ensemble of Gaussian models created using cross-validation
splits. Figure 11 displays the attributes diagram for the evidential model and
the Gaussian model (the results for other ensembles are shown in Figure S5).
Overall, both models show similar performance, although there are some
discrepancies, such as the lower performance of the Gaussian model, which
still remains within the skill region depicted in the figure. Thus, similar to
the categorical model, the evidential variation on the standard MLP for the
regression problem appears to be an effective substitute for the principal
prediction tasks at hand.

Figure 12: The spread-skill relationship is depicted using a 2D histogram,
illustrating the relationship between the standard deviation ($\sigma$) and
RMSE for the three single-task evidential models. Each row in (a-c) represents
a specific model task, while columns (i-iii) display the results for
aleatoric, epistemic, and total uncertainties. Each subpanel features a 1-to-1
line (dashed) and all panels share the color bar indicating the count of each
2D bin.

How effective is the calibration? Figure 11 shows 2D histograms quantifying
the relationship between the computed RMSE for the three model tasks and
aleatoric, epistemic, and total uncertainty. Unlike the precipitation-type
problem, the dominant uncertainty is epistemic uncertainty which contributes
the most to the total uncertainty. The 1-1 relationship is observed for
epistemic and total uncertainties when both the RMSE and relative uncertainty
are generally small, and this relationship is observed for more than half of
the testing data in each case. However, as the quantities increase, the
computed RMSE flattens out for sensible and latent heat, while for friction
velocity, the relationship appears to be linear with an initial flattening of
the RMSE for relatively higher values of uncertainty, which then continues to
grow linearly. Note also that none of the models were calibrated according to
aleatoric uncertainty; in fact, all of them were under-dispersive. Overall,
the best PITD skill score and the best fit between the 1-1 line were achieved
for friction velocity.

Figure 13: (a) PIT histograms were generated for the three tasks using (i)
aleatoric, (ii) epistemic, and (iii) total uncertainty. (b) The discard-test
diagrams illustrate the relationship between the fraction of data points
removed from the test set and the remaining subset’s RMSE for the three
single-task models shown in (i-iii).

These results are further corroborated by the PIT histogram and the discard-
test, which are shown in Figures 13(a-b), respectively. The PIT histogram also
shows that friction velocity came the closest to calibration, having the
flattest histogram of the three models, followed closely by latent heat. Both
these quantities show a slight hump towards the middle of the histogram
indicating slight under-confidence by the model. Latent heat also has a small
hump on the right side of the histogram, further indicating some over-
confidence on a subset of the data. The sensible heat model on the other hand
clearly shows two humps on the left and right sides of the histogram, showing
that the model is overconfident. All of these observations are consistent with
Figure 12. All models were subject to extensive hyper-parameter optimization,
which suggests that there are limitations to the degree of calibration
possible with the evidential regression approach.

The discard test shows that for all three models, the predicted aleatoric,
epistemic, and total uncertainties are linked to the performance of the model
through the RMSE, where the more certain data points have lower RMSE values in
each case. This is still the case even though the models possessed different
degrees of calibration. Therefore, the uncertainty value can be thresholded
such that when the model is in operation, it can be used conservatively when
model predictions are too uncertain (e.g. the model will not return a
prediction if the predicted uncertainty is larger than the threshold).

While the above results demonstrate the strong performance of the evidential
model, we must emphasize a crucial consideration regarding the actual values
of the predicted uncertainties. It is essential to be cautious as they may not
consistently align with the range of the (in-distribution) task values (Ovadia
et al., 2019). For instance, Figure S8 illustrates histograms for the best-
calibrated friction velocity model. Notice the small hump at the largest
values for the predicted epistemic uncertainty, which appears inconsistent
with the physical nature of the quantity of interest (and should be relatively
straightforward to model). This discrepancy raises a concern about the
reliability of uncertainty estimates at extreme values.

Additionally, we conducted tests using a range of other model architectures
found through optimization, which demonstrated similar performance to the ones
presented here. However, these models produced different ranges for the
histogram values at the extreme ends, despite having comparable (but not
perfect) calibration. This finding highlights the need for careful
interpretation and validation of the uncertainty estimates, as different model
architectures can lead to varying distributions of predicted uncertainties,
even if the overall calibration is consistent. Therefore, while the evidential
model shows promise in terms of performance, it is crucial to critically
evaluate and validate the predicted uncertainties to ensure their alignment
with the true nature of the underlying data.

Figure 14: Figure (a-c) illustrates the median values of the friction
velocity, sensible heat, and latent heat in 10-minute increments. The curves
displayed in each panel represent the truth, predicted, aleatoric, and
epistemic quantities. All results are derived from the test data set.

In our investigation of the SL dataset, we sought to determine whether the
predicted uncertainties align with computed errors including the underlying
physical processes governing the three tasks. Figure 14 presents the computed
median predictions and uncertainties as a function of the time of day, with
the median helping to mitigate the impact of the outliers on the mean
calculation.

The results show a distinct trend in epistemic uncertainty throughout the day.
During nighttime, the epistemic uncertainty is relatively low, indicating a
higher level of confidence in estimating the quantities. However, as the sun
rises and solar radiation intensifies, the epistemic uncertainty rapidly
increases, peaking around midday. This rise in uncertainty is attributed to
the dynamic and complex atmospheric conditions associated with strong solar
heating and convectively-driven atmospheric turbulence during this time. As
the day progresses toward night and the sun begins to set, the epistemic
uncertainty gradually decreases following the stabilization of the surface
layer.

These patterns align with our expectations, as the latent heat, which is
influenced by phase changes and processes like evaporation and condensation,
exhibits higher uncertainty compared to sensible heat, which primarily relates
to temperature changes and has a simpler measurement process. However, even a
relatively small change to the value of $\lambda$ can flip this relationship
around (Figure S9). The findings emphasize the importance of selecting an
appropriate evidential loss weight to enhance the reliability and usefulness
of uncertainty estimates, particularly for latent heat and sensible heat
tasks. Proper calibration is crucial to ensuring the model’s uncertainty
predictions during the daytime remain reliable and physically sensible.

##  4 Discussion

Proper calibration of uncertainties is crucial for regression models employed
in applications like weather forecasting (Raftery et al., 2005). However,
continuous regression calibration presents distinct challenges compared to
classification. A key difficulty is defining calibration to capture whether
uncertainties reliably indicate expected errors, as highlighted by Kuleshov et
al. (2018). Existing methods exhibit flaws in distinguishing useful and
uninformative uncertainties (Kuleshov et al., 2018). Furthermore, balancing
calibration and the dispersion of predicted probabilities is non-trivial, as
techniques like isotonic regression can enforce calibration while increasing
dispersion (Levi et al., 2022). Other challenges include modeling complex
error distributions (Harakeh et al., 2020), incorporating distinct
uncertainties (Depeweg et al., 2018), and ensuring generalization (Amershi et
al., 2019). Addressing these issues is imperative for safety-critical systems
relying on regression uncertainties.

For categorical problems, DST quantification of uncertainty offers additional
benefits. It allows representation and reasoning about uncertain or
conflicting evidence, useful in real-world datasets with ambiguous or
uncertain ground truth labels. It also enables a principled combination of
evidence from multiple sources, beneficial when dealing with different data
types or integrating information from diverse sensors or models. The ability
of a neural network to parametrically represent partial beliefs is valuable in
situations with limited information or incomplete evidence, allowing a more
nuanced representation of uncertainty.

The study also demonstrates the potential of evidential deep learning for
predictive modeling and uncertainty quantification in meteorological
applications. The evidential approach provides calibrated estimates of
uncertainty alongside predictions without requiring model ensembles.

For the winter precipitation classification task, the evidential neural
network achieves comparable accuracy to standard classifiers while quantifying
aleatoric and epistemic uncertainties. The aleatoric uncertainty, driven by
inherent randomness in the RAP model which was used as model inputs, generally
exceeds the epistemic uncertainty stemming from model generalization errors.
This highlights the benefit of decomposing uncertainties to pinpoint dominant
sources. In regions with relatively high DST uncertainty, it closely aligns
with peaks in epistemic uncertainty, demonstrating how different uncertainty
metrics capture related aspects of model confidence.

The discard test validates the use of uncertainty to filter unreliable
predictions and improve operational reliability. Setting uncertainty
thresholds allows prudent model usage by having it abstain when unsure,
preventing potentially misleading high-confidence predictions in uncertain
scenarios.

For the surface layer energy flux regression problem, proper calibration is
essential for useful uncertainties. While the evidential model’s overall
performance is just as good as what a deterministic MLP can produce,
occasional unrealistic uncertainty values were predicted, emphasizing the need
for cautious interpretation due to calibration deficiencies. Architectural
choices significantly impact uncertainty characteristics as does the dropout
rate. The time-of-day analysis reveals sensible patterns in epistemic
uncertainty variation, offering reassurance about reliability. Uncertainty
spikes during the daytime when modeling intricacies of turbulent transport and
radiation are most challenging. This physics-based behavior appropriately
reflects fundamental processes governing the quantities of interest.

A key limitation is the extensive hyperparameter tuning required for
uncertainty calibration with the evidential regression model. In contrast, the
categorical evidential model provided well-calibrated uncertainties without
excessive tuning, highlighting the greater calibration challenges for
regression tasks, and motivating continued research. One strategy that we
found useful was to first find performant architectures for the relevant
prediction tasks, and then tune the loss weight while monitoring the PITD
score (this is how Figure 10 was generated) to calibrate as closely as
possible each model.

In general, calibrating ML models to account for epistemic uncertainty
requires some kind of calibration/validation dataset or prior assumption
derived from similar modeling similar data and cannot be determined alone from
the training data. While this is more obvious for post hoc calibration methods
such as isotonic regression or conformal prediction, it is also true for
evidential methods through the tuning of the KL divergence term in the loss.
As a consequence, if the model is applied to a similar dataset with a
distribution shift, such as training on analysis data and applying the same
evidential model to a forecast, then the evidential uncertainty estimates will
be inherently underdispersive. More research is needed to identify ways to
adjust uncertainty estimates based on the nature of the domain shift.

Overall, the study establishes evidential learning as an effective framework
combining probabilistic modeling and deep learning strengths. It efficiently
provides a means for calibrating predicted uncertainties, addressing standard
network limitations, and enhancing interpretability through uncertainty
analysis. With further advances, evidential neural networks hold significant
potential across meteorological domains.

##  5 Conclusions

In conclusion, this study demonstrates the potential of evidential deep
learning as an effective technique for predictive modeling and uncertainty
quantification in weather and climate applications, for both classification
and regression. The approach synergistically integrates the capabilities of
probabilistic modeling with the representational power of deep neural
networks. This enables the models to produce well-calibrated estimates of
uncertainty alongside predictions, overcoming the limitations of standard
deep-learning approaches. The ability to quantify different sources of
uncertainty provides valuable insights into model reliability and limitations
of the training data.

The decomposition of aleatoric and epistemic uncertainties facilitates
detailed analysis to identify challenging prediction cases and opportunities
for model improvements. With calibrated uncertainty estimates, evidential
neural networks have the potential to enhance understanding of forecast
reliability and inform critical decision-making across various meteorological
domains, from real-time severe event prediction to long-term climate
projections. Given the representational flexibility, computational efficiency,
and uncertainty quantification capabilities, evidential deep learning shows
promise for tackling a diverse array of prediction and uncertainty estimation
problems in the atmospheric and climate sciences.

##  6 Software development

The Machine Integration and Learning for Earth Systems (MILES) group
Generalized Uncertainty for Earth System Science (GUESS) package (MILES GUESS)
provides tools for estimating and analyzing different sources of uncertainty
in earth system science applications. Users working in weather and climate can
leverage MILES GUESS to train neural network models that quantify multiple
uncertainty types like aleatoric, epistemic, and DST uncertainty.

The package contains layers and losses for evidential deep learning, allowing
users to build neural networks that output distributions over targets rather
than just point estimates. This enables estimating both aleatoric uncertainty
from noise in the data and epistemic uncertainty from model limitations. The
code also supports Monte Carlo dropout ensembling for epistemic uncertainty
quantification.

Once models are trained with MILES GUESS, the package provides a range of
analysis and visualization tools tailored for earth system applications. These
include PIT calibration analysis, spread-error diagrams, coverage curves, and
more. The code is primarily written in Python and is designed to integrate
seamlessly with common earth science workflows based on TensorFlow/Keras and
PyTorch. MILES GUESS is accessible at https://github.com/ai2es/miles-guess
where example Jupyter Notebooks demonstrate how to use the package.

Supplementary Materials for “Evidential Deep Learning: Enhancing Predictive
Uncertainty Estimation for Earth System Science Applications”

##  Appendix SI Derivation of the Law of Total Variance

Let X and Y be random variables defined on the same probability space. The
variance of Y can be decomposed into expected values as follows. First, the
variance of random variable X is

$\displaystyle\text{Var}[Y]$
$\displaystyle=\mathbb{E}[Y^{2}]-\mathbb{E}[Y]^{2}$ (S1)

and the law of total expectation for the random variable $X$ is defined as
below:

$\displaystyle\mathbb{E}[Y]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}[Y|X]\right]$ (S2)

Next, we apply the law of total expectation and the definition of variance to
the first term on the right-hand side of the definition of variance:

$\displaystyle\mathbb{E}[Y^{2}]$
$\displaystyle=\mathbb{E}[\mathbb{E}[Y^{2}|X]]$
$\displaystyle=\mathbb{E}[Var[Y|X]+\mathbb{E}[Y|X]^{2}]$

Continuing, we apply the law of total expectation to the second term on the
right-hand side:

$\displaystyle\mathbb{E}[Y]^{2}$
$\displaystyle=\mathbb{E}[\mathbb{E}[Y|X]]^{2}$ (S3)

Combining the terms together, we can then rearrange the expectations since the
sum of expectations is the expectation of the sums.

$\displaystyle\mathbb{E}[Y]$
$\displaystyle=\mathbb{E}[Var[Y|X]]+\mathbb{E}[Y|X]^{2}]-\mathbb{E}[\mathbb{E}[Y|X]]^{2}$
$\displaystyle=\mathbb{E}[Var[Y|X]]+(\mathbb{E}[\mathbb{E}[Y|X]^{2}]-\mathbb{\mathbb{E}[Y|X]}^{2})$

Applying the definition of variance to the second part of the term, we get the
final formula for the law of total variance:

$\displaystyle\mathbb{E}[Y]$
$\displaystyle=\mathbb{E}[Var[Y|X]]+Var[\mathbb{E}[Y|X]]$ (S4)

##  Appendix SII Decomposition for ensembles

For an ML model predicting a categorical distribution, the expected value for
the probability of class $Y_{j}$ is $p_{j}$ and the variance is
$p_{j}(1-p_{j})$. For an ensemble of K models each predicting a probability
distribution across J classes, the total variance of the predictions can be
decomposed as:

$\displaystyle Var(Y_{j})$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}p_{j,k}(1-p_{j,k})+\frac{1}{K}\sum_{k=1}^{K}(p_{j,k}-\overline{p_{j}})^{2}$
(S5)

##  Appendix SIII Decomposition for evidential classifiers

For an evidential classifier model using the Dirichlet distribution as its
output distribution,

$\displaystyle\mathbb{E}[Y_{j}|\textbf{p}]$ $\displaystyle=p_{j}$
$\displaystyle Var[Y_{j}|\textbf{p}]$ $\displaystyle=p_{j}(1-p_{j})$
$\displaystyle E[p_{j}]$ $\displaystyle=\frac{\alpha_{j}}{S}$ $\displaystyle
Var[p_{j}]$
$\displaystyle=\frac{\frac{\alpha_{j}}{S_{\alpha}}\left(1-\frac{\alpha_{j}}{S_{\alpha}}\right)}{S_{\alpha}+1}$
$\displaystyle\mathbb{E}[Var[Y_{j}|\textbf{p}]]$
$\displaystyle=\mathbb{E}[p_{j}(1-p_{j})]$
$\displaystyle=\mathbb{E}[p_{j}]-\mathbb{E}[p_{j}^{2}]$
$\displaystyle=\mathbb{E}[p_{j}]-\mathbb{E}[p_{j}]^{2}-Var[p_{j}]$
$\displaystyle=\frac{\alpha_{j}}{S}-\left(\frac{\alpha_{j}}{S}\right)^{2}-\frac{\frac{\alpha_{j}}{S_{\alpha}}\left(1-\frac{\alpha_{j}}{S_{\alpha}}\right)}{S_{\alpha}+1}$
$\displaystyle Var[\mathbb{E}[Y_{j}|\textbf{p}]]$ $\displaystyle=Var[p_{j}]$
$\displaystyle=\frac{\frac{\alpha_{j}}{S_{\alpha}}\left(1-\frac{\alpha_{j}}{S_{\alpha}}\right)}{S_{\alpha}+1}$

##  Appendix SIV Model Training and Optimization

The training process for each MLP model (both evidential and deterministic)
begins with architecture and trainable weight initialization. The training
data is then randomly shuffled. A fixed number of data points called a batch,
is selected from the shuffled training data, and the model processes the batch
to generate predictions for each point. The mean-absolute error (MAE) was used
as the training loss for deterministic regression models, while the relevant
evidential loss (categorical or regression) was used for an evidential model.
Using the training loss, the model’s weights are updated through gradient
descent with back-propagation Rumelhart et al. (1986), using a pre-defined
learning rate to iteratively reduce the error. This training process, referred
to as one epoch, covers the entire training dataset. After each epoch, the
training data is again shuffled randomly, and the process is repeated for a
specified number of epochs.

To address the class imbalance in the p-type dataset, we utilized two
techniques: upsampling and class weights. Upsampling involves artificially
increasing the number of samples from under-represented classes, like sleet
and freezing rain, by duplicating existing samples. This helps mitigate
imbalance and improve model performance on minority classes. Class weights
assign higher importance to minority classes during training to counteract the
imbalance. Typically these approaches are used individually, but we took a
combined approach enabled by hyperparameter optimization (discussed in detail
below). We treated class balancing and weighting as tunable parameters,
letting the optimizer determine when to balance classes and what weight values
to use. This joint optimization of the two strategies produced improved models
compared to using either technique alone, for both deterministic and
evidential architectures.

The hyper-parameter optimization of all models was carried out using The Earth
Computing Hyperparameter Optimization package Schreck and Gagne (2021) which
is built upon Optuna Akiba et al. (2019) and designed to leverage the high-
performance computing clusters available at NCAR. The optimization process
begins by creating a “study” and conducting the first “trial,” wherein
hyperparameters are chosen from specified ranges, the model is trained, and
its performance is evaluated in box simulations. Each trial is independent of
others, and trials are executed until the optimization converges or a
predetermined number of trials are completed. Upon completion, the study
records the outcomes and associated metadata, enabling the estimation of
hyperparameter importance for model performance.

For sampling hyperparameters, we employ the Tree-structured Parzen Estimator
(TPE) Bergstra et al. (2011). The TPE method fits a Gaussian Mixture model
(GMM) $l(x)$ to the best Mean Absolute Error (MAE) values associated with
completed trials. A second GMM is also fitted to the remaining parameter
values, and the next parameter value is chosen based on maximizing $l(x)$ /
$g(x)$. Initially, we use random sampling instead of the TPE sampler to
provide observations that the GMM models can leverage for better-informed
parameter selections. This step aids in informing the TPE sampler effectively.
For further insights into hyperparameter optimization, refer to the Optuna
documentation.

All training and hyperparameter optimization configurations (including the
best configurations used for the models in this manuscript), may be found on
our GitHub site for the project: https://github.com/ai2es/miles-guess. For the
p-type models, we found that the macro-F1 metric enabled the optimizer to find
the “best” model for sleet and freezing rain performance while sacrificing
some performance in predicting rain and snow. Other metrics including expected
and mean calibration error estimates did not perform as well.

##  Appendix SV Additional Results

Figure S1: The confusion matrix illustrates the performance of the holdout
test split for (a) the deterministic model and (b) the evidential model.

Figure S1 presents a comparison of the performance results between the
deterministic and evidential models using a confusion matrix. The confusion
matrix is shown as unnormalized to show the model’s predictive capabilities
and highlights potential classification errors.

Figure S2: Dropout calibration analysis for the p-type deterministic model.
(a) The PITD is calculated for the predicted class probabilities against
varying dropout rates. (b) The macro Brier Skill Score is computed for each
bin in the discard test and then averaged. The figure depicts the correlation
between the skill score and the dropout rate. (c) Example discard-test
outcomes obtained at different dropout rates.

The dropout calibration analysis for the p-type deterministic model is shown
in Figure S2. In Figure S2(a) PITD is computed for the predicted class
probabilities against varying dropout rates, with the best dropout rate
identified as the point when the curves reach their minimums, which are not
the same value. In Figure S2(b) the macro Brier Skill Score, which peaks just
below 0.8, demonstrates a slightly lower optimal dropout rate when results are
aggregated. Finally, in Figure S2(c) The dropout rate shows limited
sensitivity to the model’s performance, and it becomes apparent that
calibration results in lower performance on rain and snow but higher
performance on sleet and freezing rain. These observations contribute to a
better understanding of the model’s behavior and its calibration with respect
to different precipitation types.

Figure S3: Comparison of the PITD against the dropout rate for the three tasks
computed using aleatoric, epistemic, and total uncertainties. The MAE and MSE
are also plotted in each sub-panel. (a) presents the results from the data
ensemble, while (b) displays those from the deep model ensemble.

Figure S3 illustrates the dropout calibration results for the three surface-
layer flux models. In general, the two models exhibit similar calibration
behavior, but there are notable distinctions for friction velocity. The deep
ensemble shows PITD minimums at lower dropout rates compared to the data
ensemble, suggesting better calibration in this specific case. Despite this
difference, the overall shape of the calibration curves remains similar
between the two models. For latent heat and sensible heat, the calibration
results are more closely comparable, indicating consistent calibration
performance for these flux models.

Figure S4: 2D histograms illustrating the aleatoric and epistemic
uncertainties for the models predicting friction velocity, sensible heat, and
latent heat.

Figure S4 shows 2D histograms comparing aleatoric and epistemic uncertainties
for the three surface-layer flux quantities. The relationship between these
uncertainties appears most complex for friction velocity, despite it being the
simplest of the three processes.

Figure S5: Attribution curves for (a) the deep ensemble utilizing the
parametric Gaussian model and (b) the MC dropout ensemble generated using the
best member from the deep ensemble. The reliability curve is represented by
the red line. Each column displays the results for a specific model task.
Similar to Figure 3 in the main text, the diagonal, horizontal, and vertical
dashed lines in each subpanel indicate the 1-to-1 line, no-resolution line,
and climatology line, respectively. The blue-shaded area signifies skill
relative to climatology.

Figure S5(a) displays attribute diagrams computed for the deep ensemble
consisting of 20 members. From this ensemble, the best model was selected, and
Monte Carlo dropout was performed using the optimal dropout rate determined
from Figure 10, creating an ensemble of size 100. The attribute diagram for
the selected best model is shown in Figure S5(b).

Figure S6: PIT histograms are depicted for the data and deep model ensembles
in (a)(i) and (b)(i), respectively. Panels (a)(ii) and (b)(ii) display a Monte
Carlo ensemble created using the ’best’ model from the corresponding ensemble
in (i), respectively.

Figure S6 presents the PIT histograms for each uncertainty type, computed
using: the data ensemble, a Monte Carlo ensemble created using the best member
from the data ensemble, a deep ensemble, and a Monte Carlo ensemble created
using the best member from the deep ensemble. In most cases, the aleatoric
uncertainty plays the dominant role in contributing to the total uncertainty,
except for the Monte Carlo ensemble generated with the best data ensemble
model where both quantities contribute to the total uncertainty.

Figure S7: Discard-test histograms are presented for the data and deep model
ensembles in (a)(i) and (b)(i), respectively. Panels (a)(ii) and (b)(ii)
display a Monte Carlo ensemble created using the ‘best’ model from the
corresponding ensemble in (i), respectively.

Complementary to Figure S6, Figure S7 presents the discard-test histograms for
each uncertainty type, computed using: the data ensemble, a Monte Carlo
ensemble created using the best member from the data ensemble, a deep
ensemble, and a Monte Carlo ensemble created using the best member from the
deep ensemble.

Figure S8: Distributions of aleatoric and epistemic quantities predicted by
the friction velocity evidential model are depicted in histograms.

Figure S8 plots histograms of the uncertainty predictions from the friction
velocity model. The small peak observed on the far right may indicate
potential outliers resulting from imperfect calibration of the model.

Figure S9: The median value of sensible heat in 10-minute increments is shown.
The curves displayed in each panel represent the aleatoric, and epistemic
quantities. All results are derived from the test data set.

Figure S9 shows the dependence of the median aleatoric and epistemic versus
time of the day. The title in each panel lists the value of the loss weight
$\lambda$ in Equation 19 in the main text.

###### Acknowledgements.

This material is based upon work supported by the National Science Foundation
under Grant No. ICER-2019758 and by the National Center for Atmospheric
Research, which is a major facility sponsored by the National Science
Foundation under Cooperative Agreement No. 1852977. We would like to
acknowledge high-performance computing support from Cheyenne and Casper
(Computational and Information Systems Laboratory, CISL, 2020) provided by
NCAR’s Computational and Information Systems Laboratory, sponsored by the
National Science Foundation. All data sets used in this study are available at
https://zenodo.org/record/8368187. The neural networks described here and the
simulation code used to train and test the models are archived at
https://github.com/ai2es/miles-guess.

## References

  * Abdar et al. (2021) Abdar, M., and Coauthors, 2021: A review of uncertainty quantification in deep learning: Techniques, applications and challenges. Information fusion, 76, 243–297.
  * Akiba et al. (2019) Akiba, T., S. Sano, T. Yanase, T. Ohta, and M. Koyama, 2019: Optuna: A next-generation hyperparameter optimization framework. Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining.
  * Amershi et al. (2019) Amershi, S., and Coauthors, 2019: Software engineering for machine learning: A case study. 2019 IEEE/ACM 41st International Conference on Software Engineering: Software Engineering in Practice (ICSE-SEIP), IEEE, 291–300.
  * Amini et al. (2020) Amini, A., W. Schwarting, A. Soleimany, and D. Rus, 2020: Deep evidential regression. Advances in Neural Information Processing Systems, 33, 20 004–20 016.
  * Angelopoulos and Bates (2021) Angelopoulos, A. N., and S. Bates, 2021: A gentle introduction to conformal prediction and distribution-free uncertainty quantification. arXiv preprint arXiv:2107.07511.
  * Barnes and Barnes (2021) Barnes, E. A., and R. J. Barnes, 2021: Controlled abstention neural networks for identifying skillful predictions for regression problems. Journal of Advances in Modeling Earth Systems, 13 (12), e2021MS002 575.
  * Bauer et al. (2015) Bauer, P., A. Thorpe, and G. Brunet, 2015: The quiet revolution of numerical weather prediction. Nature, 525 (7567), 47–55.
  * Benjamin et al. (2016) Benjamin, S. G., and Coauthors, 2016: A north american hourly assimilation and model forecast cycle: The rapid refresh. Mon. Weather Rev., 144 (4), 1669–1694, 10.1175/MWR-D-15-0242.1.
  * Bergstra et al. (2011) Bergstra, J., R. Bardenet, Y. Bengio, and B. Kégl, 2011: Algorithms for hyper-parameter optimization. Advances in neural information processing systems, 24.
  * Bosveld et al. (2020) Bosveld, F. C., P. Baas, A. C. M. Beljaars, A. A. M. Holtslag, J. V.-G. de Arellano, and B. J. H. van de Wiel, 2020: Fifty years of atmospheric Boundary-Layer research at cabauw serving weather, air quality and climate. Bound.-Layer Meteorol., 177 (2), 583–612, 10.1007/s10546-020-00541-w.
  * Breiman (1996) Breiman, L., 1996: Bagging predictors. Machine Learning, 24 (2), 123–140.
  * Bremnes (2020) Bremnes, J. B., 2020: Ensemble postprocessing using quantile function regression based on neural networks and bernstein polynomials. Monthly Weather Review, 148 (1), 403–414.
  * Casella and Berger (2002) Casella, G., and R. L. Berger, 2002: Statistical Inference. 2nd ed., Duxbury Press.
  * Chapman et al. (2022) Chapman, W. E., L. Delle Monache, S. Alessandrini, A. C. Subramanian, F. M. Ralph, S.-P. Xie, S. Lerch, and N. Hayatbini, 2022: Probabilistic predictions from deterministic atmospheric river forecasts with deep learning. Monthly Weather Review, 150 (1), 215–234.
  * Computational and Information Systems Laboratory, CISL (2020) Computational and Information Systems Laboratory, CISL, 2020: Cheyenne: HPE/SGI ICE XA System (NCAR Community Computing). Tech. rep., National Center for Atmospheric Research. URL https://doi.org/10.5065/D6RX99HX.
  * Delaunay and Christensen (2022) Delaunay, A., and H. M. Christensen, 2022: Interpretable deep learning for probabilistic MJO prediction. Geophys. Res. Lett., 49 (16).
  * Delle Monache et al. (2013) Delle Monache, L., F. A. Eckel, D. L. Rife, B. Nagarajan, and K. Searight, 2013: Probabilistic weather prediction with an analog ensemble. Monthly Weather Review, 141 (10), 3498–3516.
  * Dempster (1968) Dempster, A. P., 1968: A generalization of bayesian inference. Journal of the Royal Statistical Society: Series B (Methodological), 30 (2), 205–232.
  * Depeweg et al. (2018) Depeweg, S., J.-M. Hernandez-Lobato, F. Doshi-Velez, and S. Udluft, 2018: Decomposition of uncertainty in bayesian deep learning for efficient and risk-sensitive learning. International Conference on Machine Learning, PMLR, 1184–1193.
  * Dietterich (2000) Dietterich, T. G., 2000: Ensemble methods in machine learning. International workshop on multiple classifier systems, Springer, 1–15.
  * Dürr et al. (2020) Dürr, O., B. Sick, and E. Murina, 2020: Probabilistic deep learning: With python, keras and tensorflow probability. Manning Publications.
  * Elmore et al. (2014) Elmore, K. L., Z. L. Flamig, V. Lakshmanan, B. Kaney, V. Farmer, H. D. Reeves, and L. P. Rothfusz, 2014: mping: Crowd-sourcing weather reports for research. Bulletin of the American Meteorological Society, 95 (9), 1335–1342.
  * Foster et al. (2021) Foster, D., D. J. Gagne, and D. B. Whitt, 2021: Probabilistic machine learning estimation of ocean mixed layer depth from dense satellite and sparse in situ observations. Journal of Advances in Modeling Earth Systems, 13 (12), e2021MS002 474.
  * Gal and Ghahramani (2016) Gal, Y., and Z. Ghahramani, 2016: Dropout as a bayesian approximation: Representing model uncertainty in deep learning. International Conference on Machine Learning.
  * Gelman et al. (2013) Gelman, A., J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin, 2013: Bayesian data analysis, Vol. 2. Chapman and Hall/CRC.
  * Gelman et al. (1995) Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin, 1995: Bayesian data analysis. Chapman and Hall/CRC.
  * Ghazvinian et al. (2021) Ghazvinian, M., Y. Zhang, D.-J. Seo, M. He, and N. Fernando, 2021: A novel hybrid artificial neural network-parametric scheme for postprocessing medium-range precipitation forecasts. Advances in Water Resources, 151, 103 907.
  * Gneiting et al. (2007) Gneiting, T., F. Balabdaoui, and A. E. Raftery, 2007: Probabilistic forecasts, calibration and sharpness. Journal of the Royal Statistical Society Series B: Statistical Methodology, 69 (2), 243–268.
  * Gneiting et al. (2005) Gneiting, T., A. E. Raftery, A. H. Westveld, and T. Goldman, 2005: Calibrated probabilistic forecasting using ensemble model output statistics and minimum crps estimation. Monthly Weather Review, 133 (5), 1098–1118.
  * Gordon and Barnes (2022) Gordon, E. M., and E. A. Barnes, 2022: Incorporating uncertainty into a regression neural network enables identification of decadal state-dependent predictability in cesm2. Geophysical Research Letters, 49 (15), e2022GL098 635.
  * Guillaumin and Zanna (2021) Guillaumin, A. P., and L. Zanna, 2021: Stochastic-deep learning parameterization of ocean momentum forcing. Journal of Advances in Modeling Earth Systems, 13 (9), e2021MS002 534.
  * Hamill (2001) Hamill, T. M., 2001: Interpretation of rank histograms for verifying ensemble forecasts. Monthly Weather Review, 129 (3), 550–560.
  * Harakeh et al. (2020) Harakeh, A., M. Smart, and S. L. Waslander, 2020: Bayesod: A bayesian approach for uncertainty estimation in deep object detectors. 2020 IEEE International Conference on Robotics and Automation (ICRA), IEEE, 87–93.
  * Haupt et al. (2021) Haupt, S. E., W. Chapman, S. V. Adams, C. Kirkwood, J. S. Hosking, N. H. Robinson, S. Lerch, and A. C. Subramanian, 2021: Towards implementing artificial intelligence post-processing in weather and climate: proposed actions from the oxford 2019 workshop. Philosophical Transactions of the Royal Society A, 379 (2194), 20200 091.
  * Haynes et al. (2023) Haynes, K., R. Lagerquist, M. McGraw, K. Musgrave, and I. Ebert-Uphoff, 2023: Creating and evaluating uncertainty estimates with neural networks for environmental-science applications. Artificial Intelligence for the Earth Systems, 2 (2), 220 061.
  * Herman and Schumacher (2018) Herman, G. R., and R. S. Schumacher, 2018: Money doesn’t grow on trees, but forecasts do: Forecasting extreme precipitation with random forests. Monthly Weather Review, 146 (5), 1571–1600.
  * Hoffman et al. (2013) Hoffman, M. D., D. M. Blei, C. Wang, and J. Paisley, 2013: Stochastic variational inference. Journal of Machine Learning Research.
  * Jsang (2018) Jsang, A., 2018: Subjective Logic: A formalism for reasoning under uncertainty. Springer Publishing Company, Incorporated.
  * Karpatne (2018) Karpatne, A., 2018: Theory-guided data science: A new paradigm for scientific discovery from data. IEEE Transactions on Knowledge and Data Engineering, 30 (6), 1009–1021.
  * Kendall and Gal (2017) Kendall, A., and Y. Gal, 2017: What uncertainties do we need in bayesian deep learning for computer vision? Advances in neural information processing systems, 30.
  * Koenker (2005) Koenker, R., 2005: Quantile Regression. Cambridge University Press.
  * Kotz et al. (2004) Kotz, S., N. Balakrishnan, and N. L. Johnson, 2004: Continuous multivariate distributions, Volume 1: Models and applications, Vol. 1. John Wiley & Sons.
  * Kuleshov et al. (2018) Kuleshov, V., N. Fenner, and S. Ermon, 2018: Accurate uncertainties for deep learning using calibrated regression. International conference on machine learning, PMLR, 2796–2804.
  * Lakshminarayanan et al. (2016) Lakshminarayanan, B., A. Pritzel, and C. Blundell, 2016: Simple and scalable predictive uncertainty estimation using deep ensembles. Advances in neural information processing systems, 29.
  * Leith (1974) Leith, C. E., 1974: Theoretical skill of monte carlo forecasts. Monthly Weather Review, 102 (6), 409 – 418, https://doi.org/10.1175/1520-0493(1974)102¡0409:TSOMCF¿2.0.CO;2.
  * Levi et al. (2022) Levi, D., L. Gispan, N. Giladi, and E. Fetaya, 2022: Evaluating and calibrating uncertainty prediction in regression tasks. Sensors, 22 (15), 5540.
  * Liu et al. (2020) Liu, J., Z. Lin, S. Padhy, D. Tran, T. Bedrax Weiss, and B. Lakshminarayanan, 2020: Simple and principled uncertainty estimation with deterministic deep learning via distance awareness. Advances in Neural Information Processing Systems, 33, 7498–7512.
  * MacKay (2003) MacKay, D. J., 2003: Information theory, inference and learning algorithms. Cambridge university press.
  * Malinin and Gales (2018) Malinin, A., and M. Gales, 2018: Predictive uncertainty estimation via prior networks. Advances in neural information processing systems, 31.
  * May et al. (2022) May, R. M., and Coauthors, 2022: MetPy: A meteorological python library for data analysis and visualization. Bull. Am. Meteorol. Soc., 103 (10), E2273–E2284, 10.1175/BAMS-D-21-0125.1.
  * McCandless et al. (2022) McCandless, T., D. J. Gagne, B. Kosović, S. E. Haupt, B. Yang, C. Becker, and J. Schreck, 2022: Machine learning for improving surface-layer-flux estimates. Boundary-Layer Meteorology, 185 (2), 199–228.
  * McGovern et al. (2017) McGovern, A., K. L. Elmore, D. J. Gagne, S. E. Haupt, C. D. Karstens, R. Lagerquist, T. Smith, and J. K. Williams, 2017: Using artificial intelligence to improve real-time decision-making for high-impact weather. Bulletin of the American Meteorological Society, 98 (10), 2073–2090.
  * Meinert et al. (2023) Meinert, N., J. Gawlikowski, and A. Lavin, 2023: The unreasonable effectiveness of deep evidential regression. Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 37, 9134–9142.
  * Meinert and Lavin (2021) Meinert, N., and A. Lavin, 2021: Multivariate deep evidential regression. arXiv preprint arXiv:2104.06135.
  * Muñoz-Esparza et al. (2022a) Muñoz-Esparza, D., C. Becker, J. A. Sauer, D. J. Gagne, J. Schreck, and B. Kosović, 2022a: On the application of an observations-based machine learning parameterization of surface layer fluxes within an atmospheric large-eddy simulation model. Journal of Geophysical Research: Atmospheres, 127 (16), e2021JD036 214.
  * Muñoz-Esparza et al. (2022b) Muñoz-Esparza, D., C. Becker, J. A. Sauer, D. J. Gagne, II, J. Schreck, and B. Kosović, 2022b: On the application of an observations‐based machine learning parameterization of surface layer fluxes within an atmospheric large‐eddy simulation model. J. Geophys. Res., 127 (16), 10.1029/2021jd036214.
  * Murphy (2007) Murphy, K. P., 2007: Conjugate bayesian analysis of the gaussian distribution. def, 1 (2$\sigma$2), 16.
  * Nadav-Greenberg and Joslyn (2009) Nadav-Greenberg, L., and S. L. Joslyn, 2009: Uncertainty forecasts improve decision making among nonexperts. Journal of Cognitive Engineering and Decision Making, 3 (3), 209–227, 10.1518/155534309X474460, https://doi.org/10.1518/155534309X474460.
  * Neal (2012) Neal, R. M., 2012: Bayesian learning for neural networks, Vol. 118. Springer Science & Business Media.
  * Nipen and Stull (2011) Nipen, T., and R. Stull, 2011: Calibrating probabilistic forecasts from an nwp ensemble. Tellus A: Dynamic Meteorology and Oceanography, 63 (5), 858–875.
  * Nix and Weigend (1994) Nix, D. A., and A. S. Weigend, 1994: Estimating the mean and variance of the target probability distribution. Proceedings of 1994 ieee international conference on neural networks (ICNN’94), IEEE, Vol. 1, 55–60.
  * Oh and Shin (2022) Oh, D., and B. Shin, 2022: Improving evidential deep learning via multi-task learning. Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 36, 7895–7903.
  * Ovadia et al. (2019) Ovadia, Y., and Coauthors, 2019: Can you trust your model’s uncertainty? evaluating predictive uncertainty under dataset shift. Advances in neural information processing systems, 32.
  * Raftery et al. (2005) Raftery, A. E., T. Gneiting, F. Balabdaoui, and M. Polakowski, 2005: Using bayesian model averaging to calibrate forecast ensembles. Monthly weather review, 133 (5), 1155–1174.
  * Rasmussen and Williams (2006) Rasmussen, C. E., and C. K. I. Williams, 2006: Gaussian Processes for Machine Learning. The MIT Press.
  * Rasp and Lerch (2018) Rasp, S., and S. Lerch, 2018: Neural networks for postprocessing ensemble weather forecasts. Monthly Weather Review, 146 (11), 3885–3900.
  * Romano et al. (2019) Romano, Y., E. Patterson, and E. Candes, 2019: Conformalized quantile regression. Advances in neural information processing systems, 32.
  * Rumelhart et al. (1986) Rumelhart, D. E., G. E. Hinton, and R. J. Williams, 1986: Learning representations by back-propagating errors. Nature, 323 (6088), 533–536.
  * Scheuerer et al. (2020) Scheuerer, M., M. B. Switanek, R. P. Worsnop, and T. M. Hamill, 2020: Using artificial neural networks for generating probabilistic subseasonal precipitation forecasts over california. Monthly Weather Review, 148 (8), 3489–3506.
  * Schreck and Gagne (2021) Schreck, J. S., and D. J. Gagne, 2021: Earth computing hyperparameter optimization. GitHub, URL https://github.com/NCAR/echo-opt, https://zenodo.org/record/7713291.
  * Schulz and Lerch (2022) Schulz, B., and S. Lerch, 2022: Machine learning methods for postprocessing ensemble forecasts of wind gusts: A systematic comparison. Monthly Weather Review, 150 (1), 235–257.
  * Sensoy et al. (2020) Sensoy, M., V. Deli’c, and L. Kaplan, 2020: Uncertainty quantification in deep learning. arXiv preprint arXiv:2011.06225.
  * Sensoy et al. (2018) Sensoy, M., L. Kaplan, and M. Kandemir, 2018: Evidential deep learning to quantify classification uncertainty. Advances in Neural Information Processing Systems, Vol. 31.
  * Shaker and Hüllermeier (2020) Shaker, M. H., and E. Hüllermeier, 2020: Aleatoric and epistemic uncertainty with random forests. Advances in Intelligent Data Analysis XVIII: 18th International Symposium on Intelligent Data Analysis, IDA 2020, Konstanz, Germany, April 27–29, 2020, Proceedings 18, Springer, 444–456.
  * Shepherd (2009) Shepherd, J. G., 2009: Geoengineering the climate: science, governance and uncertainty. Royal Society.
  * Soleimany et al. (2021) Soleimany, A. P., A. Amini, S. Goldman, D. Rus, S. N. Bhatia, and C. W. Coley, 2021: Evidential deep learning for guided molecular property prediction and discovery. ACS Central Science, 7 (8), 1356–1367.
  * Srivastava et al. (2014) Srivastava, N., G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, 2014: Dropout: a simple way to prevent neural networks from overfitting. The journal of machine learning research, 15 (1), 1929–1958.
  * Stankeviciute et al. (2021) Stankeviciute, K., A. M Alaa, and M. van der Schaar, 2021: Conformal time-series forecasting. Advances in neural information processing systems, 34, 6216–6228.
  * Ulmer et al. (2023) Ulmer, D. T., C. Hardmeier, and J. Frellsen, 2023: Prior and posterior networks: A survey on evidential deep learning methods for uncertainty estimation. Transactions on Machine Learning Research.
  * Vannitsem et al. (2018) Vannitsem, S., D. S. Wilks, and J. Messner, 2018: Statistical postprocessing of ensemble forecasts. Elsevier.
  * Vannitsem et al. (2020) Vannitsem, S., and Coauthors, 2020: Statistical postprocessing for weather forecasts–review, challenges and avenues in a big data world. Bulletin of the American Meteorological Society, 1–44.