# Fusion framework and multimodality for the Laplacian approximation of
Bayesian neural networks

Magnus Malmström, Isaac Skog , Daniel Axehill , Fredrik Gustafsson  This work
is supported by Sweden’s innovation agency, Vinnova, through project iQDeep
(project number 2018-02700). Magnus Malmström, Daniel Axehill, and Fredrik
Gustafsson, is with Linköping University (e-mail: {magnus.malmstrom,
daniel.axehill<EMAIL_ADDRESS>Isaac Skog, is with Uppsala
University (e-mail<EMAIL_ADDRESS>
###### Abstract

This paper considers the problem of sequential fusion of predictions from
neural networks (nn) and fusion of predictions from multiple nn. This fusion
strategy increases the robustness, i.e., reduces the impact of one incorrect
classification and detection of outliers the nn has not seen during training.
This paper uses Laplacian approximation of Bayesian nns (bnns) to quantify the
uncertainty necessary for fusion. Here, an extension is proposed such that the
prediction of the nn can be represented by multimodal distributions. Regarding
calibration of the estimated uncertainty in the prediction, the performance is
significantly improved by having the flexibility to represent a multimodal
distribution. Two class classical image classification tasks, i.e., mnist and
cfar10, and image sequences from camera traps of carnivores in Swedish forests
have been used to demonstrate the fusion strategies and proposed extension to
the Laplacian approximation.

##  I Introduction

This paper studies how to fuse the prediction from multiple neural networks
(nn) classifiers. Two problem scenarios are considered. Firstly, the problem
of fusing the predictions from multiple nn classifiers that attempt to
classify the same object. Secondly, to fuse the predictions from a single nn
classifier that are given a sequence of inputs known to belong to the same
class.

One often uses multiple algorithms and methods that work in parallel to have
redundancy in a decision process. The predictions are then combined to make
the decision more robust. However, without knowledge of the uncertainty in the
prediction, it is unclear how the different predictions should be weighted.
Hence, it is necessary to have good knowledge regarding the uncertainty in the
prediction from the model.

In recent years, nns have had great success generating images from text [1],
mastering board games such as GO [2], and in various control tasks [3, 4].
Despite their tremendous success, there is still limited use of nns in safety-
critical applications, e.g., medical imaging and autonomous driving [5, 6, 7].
One of the main reasons nns have not revolutionized autonomous vehicles yet is
the lack of knowledge of the uncertainty in their predictions. Take, for
example, the infamous accident in 2018 by one of Uber’s autonomous vehicles
[8]. Here, the lack of reliable classification uncertainty by the nn of the
surrounding object was a contributing factor to the fatal outcome of the
accident. To create a more robust decision, the fusion of the prediction from
sensors such as cameras and lidars could have been helpful in this scenario.

The problem of quantifying the uncertainty in the prediction of nn has lately
been gaining interest with many suggested methods [9, 10, 11, 12, 13, 14]. See
[15] for a survey over different methods. Broadly, one can separate these
methods to quantify uncertainty in the prediction into two categories. The
first category of methods is based on designing the structure of the nn such
that it learns its own uncertainty [16, 17, 18, 19]. The second category is
based on creating an ensemble of predictions, from which the uncertainty in
the prediction could be computed [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30].

This paper will consider methods from the second category. The naive approach
would be to independently train multiple nn on the same task [21]. This method
is frequently referred to as deep ensemble. The major disadvantage of this
approach is that it is very computationally costly since it requires training
multiple nn where only training one takes a lot of computational resources.
The so-called test time augmentation is another method to create the ensemble.
Here, a single nn creates an ensemble of predictions by predicting the
slightly modified version of the input [26]. Test time augmentation is typical
in medical image classification, where little data is available. Here, the
modification could, for example, be to rotate the image. However, this method
does not consider the correlation from using the same nn to predict the output
of all the augmented inputs.

Instead of assuming the parameters of the nn to be fixed, one can assume that
they have some distribution from which one can draw samples to create the
ensemble, i.e., a Bayesian nn (bnn) [31]. However, to train a bnn is not
straightforward and can be computationally expensive, and hence, one is often
referred to use some approximation method. For example, values of the
parameters could be sampled during the later part of the training [24, 25], or
the distribution can be assumed to be given by some already used
regularization techniques from which samples can be drawn during inference
[22, 23]. Another well-used approximation is the so-called Laplacian
approximation of a bnn, where the curvature of the likelihood function gives
the distribution of the parameters [30, 32, 33, 31, 34]. In particular, this
paper will focus on the linearized Laplacian approximation (lla) [29, 35],
where the so-called delta method is used to model the distribution for the
prediction rather than the distribution for the parameters [36, 37, 38].

A shortcoming with Laplacian approximation is that it is a local approach that
only quantifies the uncertainty in a neighborhood of some given parameters.
Hence, it lacks the expressiveness to represent a multimodal distribution. To
solve this shortcoming, we propose the ensemble lla (ella). The proposed
method combines deep-ensemble and Laplacian approximation, training multiple
nns making a Laplacian approximation for each of the nns.

The contribution of this paper is three-fold. Firstly, the paper presents a
method to fuse predictions from multiple classifications, which achieves a
more accurate classifier. Both for the case using multiple classifiers and
when the classification comes from a sequence of inputs known to belong to the
same class. Here, fusing the information is shown to lead to more robust
decisions. Secondly, the delta method for the Laplacian approximation of bnn
presented in [29] is extended such that the estimated probability can
represent multimodal distributions. Thirdly, the presented methods demonstrate
efficiency in detecting out-of-distribution examples.

##  II Fusion of Classifiers

We consider the problem of fusing several classifiers’ predicted probability
mass functions (pmf) from a Bayesian perspective. Let $y\in\\{1,\ldots,M\\}$
and $x\in\mathbb{R}^{n_{x}}$ denote the class label and input, respectively.
Further, assume that the joint probability distribution $p(X,Y)$ can describe
the underlying data-generating process. There are $C$ training data sets
generated as

$\displaystyle\begin{split}\mathcal{T}^{(c)}=\\{x_{i},y_{i}\\}_{i=1}^{N_{c}},\quad(x_{i},y_{i})\stackrel{{\scriptstyle\text{\tiny
iid}}}{{\sim}}p(X,Y),\quad c=1,2,\\!\dots,C.\end{split}$ (1)

available to train the classifiers. During the training-phase, the classifiers
try to learn (identify) a function $f(x|\mathcal{T}^{(c)})$ that approximates
the conditional distribution $p(Y=y|X=x)$ from the training data
$\mathcal{T}^{(c)}$. That is, after the classifiers have been trained, we have
the pmf approximations

$\displaystyle p(Y=y|X=x)\approx f(x|\mathcal{T}^{(c)}),\quad c=1,2,\dots,C.$
(2)

In the classification phase, $L$ inputs (assumed to be known to belong to the
same class) are generated as $x^{\star}_{l}\sim p(X,Y=y^{\star})$ for
$l=1,2,\dots,L$. Inserting the new inputs into the classifiers yields the pmf
estimates

$\displaystyle\hat{p}_{lc}(y^{\star})=f(x^{\star}_{l}|\mathcal{T}^{(c)}),\quad
l=1,\dots L,\ c=1,\dots,C.$ (3)

The fundamental question is how to fuse the estimates $\hat{p}_{lc}(y)$ given
that we know that all $x^{\star}_{l}$ are input samples corresponding to the
same, but unknown class $y^{\star}$? That is, to estimate

$\displaystyle p(y^{\star}|x_{1:L}^{\star})\triangleq
p(y^{\star}|x_{1:L}^{\star},\mathcal{T}^{(1)},\\!\\!\ldots,\\!\mathcal{T}^{(C)}).$
(4)

As a remark, most practical use cases have either $L=1$ (apply several
classifiers to one input) or $C=1$ (apply the same classifier to several
images, e.g., from a video stream), but both cases and their generalization
will be treated in parallel in the sequel.

###  II-A Fusion using non-parametric models

If all inputs and training data sets are independent, then the Bayes rule
gives that

$\displaystyle
p(y|x_{1:L}^{\star})\propto\prod_{l=1}^{L}\prod_{c=1}^{C}p_{lc}(y).$ (5)

Thus, one reasonable way to fuse the pmf estimates is

$\displaystyle\hat{p}(y|x_{1:L}^{\star})=\frac{\prod_{l=1}^{L}\prod_{c=1}^{C}\hat{p}_{lc}(y)}{\sum_{m=1}^{M}\prod_{l=1}^{L}\prod_{c=1}^{C}\hat{p}_{lc}(m)}.$
(6)

If the relative quality (uncertainty) of the estimates $\hat{p}_{lc}(y)$ are
known and represented by the scalar weights $w_{lc}$, $\sum_{lc}w_{lc}=1$,
they may also be fused as

$\displaystyle\hat{p}(y|x_{1:L}^{\star})\propto\prod_{l=1}^{L}\prod_{c=1}^{C}\hat{p}_{lc}^{w_{lc}}(y).$
(7)

Note that this fusion rule is equivalent to log-linear pooling [39]. Compared
to our Bayesian approach, this method is rather ad-hoc, and the weights
$w_{lc}$ reflect reliability in each classifier rather than a stochastic
measure of the expected performance. Furthermore, the weights do not depend on
the input, so the relative strength of each classifier in different regions is
not exploited.

###  II-B Fusion using parametric models

So far, we have assumed a black-box structure of the classifier, which only
depends on the training data set
$\hat{p}_{lc}(y)=f(x^{\star}_{l}|\mathcal{T}^{(c)})$. This assumption covers
simple classifiers such as the nearest neighbor. In the sequel, we will
consider parametric functions, such as nns, where the classifier explicitly
depends on a parameter $\theta$, which is estimated as $\hat{\theta}^{c}_{N}$
from the training data $\mathcal{T}^{(c)}$. This dependence will be made
explicit in the sequel, so we denote the conditional probabilities

$\displaystyle\hat{p}_{lc}(y^{\star})\triangleq
f(x^{\star}_{l}|\hat{\theta}^{c}_{N}),\quad\forall l,c.$ (8)

The uncertainty in $\hat{\theta}^{c}_{N}$ implies that the conditional
probability $f(x^{\star}_{l}|\hat{\theta}^{c}_{N})$ is itself a distribution.
This fact is the key point with our approach that allows the fusion of
classifiers where the reliability of each one may depend on both the input and
the quality of the training data, as reflected in the uncertainty of the
parameter estimate. A requirement to be able to do this is to specify the
posterior distribution of the parameters $p(\theta|\mathcal{T}^{(c)})$.

Now, a distribution of distributions is quite a complex object to handle,
where no analytical expressions can be expected. One remedy is to represent
the distribution of the estimated parameters $\hat{\theta}^{c}_{N}$ with Monte
Carlo (mc) samples $\theta^{c(k)}$. That is, we use the approximation

$\displaystyle\theta^{c(k)}$ $\displaystyle\sim
p(\theta|\mathcal{T}^{(c)}),\quad k=1,2\dots,K,$ (9a)
$\displaystyle\hat{p}_{lc}(y^{\star})$
$\displaystyle\approx\frac{1}{K}\sum_{k=1}^{K}f(x^{\star}_{l}|\theta^{c(k)}),\quad\forall
l,c,$ (9b) $\displaystyle\hat{p}(y^{\star}|x_{1:L}^{\star})$
$\displaystyle\propto\prod_{l=1}^{L}\prod_{c=1}^{C}\hat{p}_{lc}(y^{\star}).$
(9c)

Here, $K$ denotes the number of samples used in the mc sampling. One major
shortcoming with this formulation is that every classifier is represented with
only a point estimate, which does not utilize the fact that we are dealing
with a distribution of distributions. This is because the fusion is done after
the estimation using the mc samples. A remedy for this problem is to perform
the fusion before the sampling stage. However, quantifying
$p(\theta|\mathcal{T}^{(c)})$ is a difficult problem, hence also to fuse the
predictions from all the classifiers, both over the classifiers and also over
the sequence of inputs.

Another shortcoming with approximating $p(y^{\star}|x^{\star}_{1:L})$ by mc
sampling from (9) is that it comes with a high computational cost. This cost
is a result from the approximation requiring drawing samples from a high
dimensional Gaussian distribution and evaluating the whole nn multiple times.
The high computational complexity is particularly true for nn classifiers,
which have a huge dimension of the parameter space. In [29], a remedy to this
is presented, which reduces the dimension in the sampling to an
$M$-dimensional space. In this paper, it is done using the so-called delta
method [38, 36]. Then, the predicted values before the normalization are
sampled instead of the values of the parameters $\theta$ of the nn. Hence, no
forward passes of the nn are required, and the dimension of the distribution
from which the samples are drawn is significantly smaller. Apart from
decreasing the computational complexity, this might also lead to a more
straightforward fusion strategy compared to (9). This method to propagate the
uncerainty will further be described in the following sections.

##  III Neural Network Classifiers

This section will outline some basic concepts that are required for the two
presented methods in this paper to approximate the conditional pmf given a set
of classifiers, i.e., to estimate $\hat{p}(y^{\star}|x^{\star}_{1:L})$. In
this way, uncertainty in the parameters from the estimation step can be
preserved and transformed to, first, the output from the last layers of the
nns and, second, to the class probabilities. The goal is to estimate the
probability distribution of the class probabilities, including cross-
correlations.

###  III-A The Softmax Operator

It will be assumed that the nn classifier has a softmax function in the output
layer. This ensures that the model $f(x|\theta^{(c)})$ fulfills the properties
associated with a pmf, i.e., $f_{m}(x|\theta^{(c)})\geq 0$ $\forall m$ and
$\sum_{m}f_{m}(x|\theta^{(c)})=1$. Thus, it is assumed that the nn is
structured as

$f(x|\theta^{(c)})=\text{softmax}\left(g(x|\theta^{(c)})\right)$ (10a) where
$\text{softmax}(z)\triangleq\frac{1}{\sum_{m=1}^{M}e^{z_{m}}}\begin{bmatrix}e^{z_{1}}\\\
\vdots\\\ e^{z_{M}}\end{bmatrix}.$ (10b)

Here the family of functions $g(x|\theta)$ is unconstrained, while the softmax
function maps these functions onto the interval $[0,1]$.

It should here be noted that the softmax function is invariant to
translations, so $\text{softmax}(z)=\text{softmax}(z+\alpha)$ for all
$\alpha$. Here, a transformation is needed to make different classifiers
comparable. The chosen translation is arbitrary, but a natural choice is

$\bar{g}(x|\theta^{(c)})\triangleq
g(x|\theta^{(c)})-\max_{k}g_{k}(x|\theta^{(c)})\leq 0.$ (11)

This is a sound choice from a numerical point of view since the exponential
functions will operate on numbers smaller than zero, and there is no risk of
numerical overflow. However, our main reason for transformation is to enable
fusion.

###  III-B Monte Carlo Sampling of Class Probabilities

Define $z\in\mathbb{R}^{M}$ as the value of the classifiers before the
normalization done by the softmax function. Assume the distribution
$p(z^{\star}|x^{\star}_{1:L})$ has been quantified, i.e., a distribution for
the fusion of the classifiers before the normalization. By combining the mc
sampling technique in (9), we get an mc approximation of the class
probabilities in the following way:

$\displaystyle z^{(k)}$ $\displaystyle\sim p(z^{\star}|x^{\star}_{1:L}),\quad
k=1,2,\dots,K,$ (12a) $\displaystyle\hat{p}^{(k)}$
$\displaystyle=\mathrm{softmax}\big{(}z^{(k)}\big{)},$ (12b)
$\displaystyle\hat{p}(y^{\star}|x^{\star}_{1:L})$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\hat{p}^{(k)}.$ (12c)

The uncertainty of the estimated class labels
$\hat{p}(y^{\star}|x^{\star}_{1:L})$ is here explicitly represented by the
point cloud $\hat{p}^{(k)}$. From this, various probabilities can be computed,
e.g., the probability that $\hat{p}(y^{\star}|x^{\star}_{1:L})$ is true (count
the fraction of sample vectors $\hat{p}^{(k)}$ whose $c$’th element is the
largest one), or the risk that we decide $c_{1}$ while $c_{2}\neq c_{1}$ is
true.

Hence, we have laid the ground for the proposed Bayesian approach to estimate
the class probability distribution. The remaining task is to estimate
$p(z^{\star}|x^{\star}_{1:L})$, preferably without adding a substantial amount
of computations in the training and classification stages.

##  IV Uncertainty quantification

Next, a method to quantify $p_{lc}(z^{\star})$ for nns is described. After
that, the presentation of an extension where the correlation between
classifiers is included. This extension is required to estimate
$p(z^{\star}|x^{\star}_{1:L})$.

###  IV-A Laplacian approximation and the delta method

Representing $p(\theta|\mathcal{T}^{(c)})$ is challenging. Hence, one often
has to rely on approximations to represent $p(\theta|\mathcal{T}^{(c)})$. One
such approximation is the Laplacian approximation of bnns [31, 34]. It is a
local linear approach, that approximates the uncertainty in the parameters
using the curvature of the likelihood function around the estimated parameters
$\hat{\theta}_{N}^{c}$. Assume some prior of the parameters
$p(\theta)=\mathcal{N}(\theta;0,P_{0})$. Then the Laplacian approximation of
the posterior $p(\theta|\mathcal{T}^{(c)})$ yields [31]

$\displaystyle p(\theta|\mathcal{T}^{(c)})$
$\displaystyle=\mathcal{N}(\theta;\hat{\theta}_{N}^{c},P_{N}^{\theta,c}),$
(13a) $\displaystyle P_{N}^{\theta,c}$
$\displaystyle=\bigg{(}-\frac{\partial^{2}L_{N}(\theta)}{\partial\theta^{2}}\bigg{|}_{\theta=\hat{\theta}_{N}^{c}}+P_{0}^{-1}\bigg{)}^{-1}.$
(13b)

Here

$L_{N}(\theta)=\sum_{n=1}^{N}\ln f_{y_{n}}(x_{n}|\theta)$ (14)

denotes the cross-entropy likelihood function [40] where $y_{n}$ is used as an
index operator for the subscript $m$ of $f_{m}(x|\theta)$. The maximum a
posteriori estimate of the parameter $\hat{\theta}_{N}$ is given by

$\hat{\theta}_{N}^{c}=\operatorname*{arg\,max}_{\theta}p(\theta|\mathcal{T}^{(c)})=\operatorname*{arg\,max}_{\theta}L_{N}(\theta)+\ln
p(\theta),$ (15)

The choice of the prior distribution $p(\theta)=\mathcal{N}(\theta;0,P_{0})$,
corresponds to the so-called $L2$-regularization.

Under the assumption that the true model belongs to the considered model set,
Bernstein-von Mises theorem [41] gives us that the distribution of the maximum
posterior estimate asymptotically coincides with the Laplacian approximation
in (13). Furthermore, the second derivative of the likelihood function
approximates the inverse of the Fisher information matrix. In [29], it was
shown that for a classification problem, the Fisher information is given by

$\mathcal{I}^{\theta,c}\simeq\sum_{n=1}^{N}\sum_{m=1}^{M}\eta_{m,n}^{c}\frac{\partial
g_{m}\\!(x_{n}|\theta)}{\partial\theta}\Bigg{|}_{\theta=\hat{\theta}_{N}^{c}}\bigg{(}\\!\frac{\partial
g_{m}\\!(x_{n}|\theta)}{\partial\theta}\Bigg{|}_{\theta=\hat{\theta}_{N}^{c}}\\!\bigg{)}^{\\!\\!\top\\!}$
(16a) where $\eta_{m,n}^{c}\triangleq
f_{m}(x_{n}|\hat{\theta}_{N}^{c})(1-f_{m}(x_{n}|\hat{\theta}_{N}^{c})).$ (16b)

The delta method relies on linearizing a nonlinear model
$g(x^{\star}_{l}|\hat{\theta}_{N}^{c})$ to propagate the uncertainty in the
parameters to uncertainty in the prediction. Using the delta method

$p_{lc}(z^{\star})\approx\mathcal{N}\big{(}z;\hat{g}_{N}^{lc},P_{N}^{g,lc})$
(17a) where $\hat{g}_{N}^{lc}=\text{E}\\{z\\}\simeq
g(x^{\star}_{l}|\hat{\theta}_{N}^{c})$ (17b) and
$\begin{split}P_{N}^{g,lc}&=\mathrm{Cov}\\{z\\}\\\
&\simeq\bigg{(}\frac{\partial}{\partial\theta}g(x^{\star}_{l}|\theta)\big{|}_{\theta=\hat{\theta}_{N}^{c}}\bigg{)}^{\top}\\!P^{\theta,lc}_{N}\frac{\partial}{\partial\theta}g(x^{\star}_{l}|\theta)\big{|}_{\theta=\hat{\theta}_{N}^{c}}.\end{split}$
(17c)

Since the covariance is invariant to translations, the transformed functions
in (11) also have a Gaussian distribution

$p_{lc}(z^{\star})\approx\mathcal{N}\big{(}z;\hat{\bar{g}}_{N}^{lc},P_{N}^{g,lc}\big{)},$
(18a) where
$\displaystyle\hat{\bar{g}}_{N}^{lc}\simeq\bar{g}(x^{\star}_{l}|\hat{\theta}_{N}^{c}).$
(18b)

To summarize, the complete method to quantify uncertainty in the prediction
for one classifier is based on two linearizations. One to approximate the
uncertainty in the parameters, i.e., the Laplacian approximation in (13b), and
one to propagate the uncertainty to the prediction, i.e., the delta method in
(17). The complete method is here referred to as the lla.

###  IV-B Approximating the covariance

An nn often has millions of parameters, so it can be intractable to represent
the covariance $P^{\theta,c}_{N}$ for all the parameters. For example, the
covariance can be too large to be stored in the computer’s local memory.
Hence, one often has to rely on approximations of the covariance, e.g.,
neglecting some correlation between parameters [32], assuming some parameters
to be fixed [33], or neglecting all the cross-correlation between the
parameters (mean-field variational inference) [20]. In this work, we follow
the second approach where it will be assumed that the parameters of the
earlier part of the nn are fixed and only the parameters in the later layers
are assumed to be learned [33, 42, 29]. This approximation might be more or
less accurate depending on the number of included layers. A scaling of
$P^{\theta,c}_{N}$ with factor $T_{\theta}\geq 1$ can be introduced to
compensate somewhat for the approximation error. The scaling can be estimated
from validation data as is typically done for temperature scaling [43].

###  IV-C Multiple classifiers

There might be correlations between the classifications when using a single
classifier to classify a set of inputs $x^{\star}_{l}$, $l=1,\ldots,L$, known
to belong to the same class $y^{\star}$. This correlation can be taken into
account using

$\displaystyle\hat{p}_{c}(\zeta^{c}|x_{1:L}^{\star})\approx\mathcal{N}\big{(}\zeta^{c}(x^{\star});\mathbf{\hat{z}}^{c},R_{c}\big{)}$
(19)

where $\zeta^{c}\in\mathbb{R}^{LM}$ denotes the classifications from the
$c$’th classifier for the $L$ different inputs

$\displaystyle\hat{\zeta}^{c}=\begin{bmatrix}\hat{\bar{g}}_{N}^{1c}\\\
\vdots\\\ \hat{\bar{g}}_{N}^{Lc}\end{bmatrix}$ (20)

and the block $[R]_{i,j}\in\mathbb{R}^{M,M}$, $i,j=1,\ldots,L$, of the
covariance matrix is given by

$[R_{c}]_{i,j}=\frac{\partial}{\partial\theta}g(x^{\star}_{i}|\theta)^{\top}\big{|}_{\theta=\hat{\theta}_{N}^{c}}\\!P^{\theta}_{N}\frac{\partial}{\partial\theta}g(x^{\star}_{j}|\theta)\big{|}_{\theta=\hat{\theta}_{N}^{c}}.$
(21)

Extending the description by using multiple classifiers to classify the
sequence as

$\displaystyle
p(\zeta|x_{1:L}^{\star})\approx\mathcal{N}(\zeta;\mathbf{\hat{\zeta}},R)$
(22a) where
$\displaystyle\hat{\zeta}=\begin{bmatrix}\mathbf{\hat{\zeta}}^{1}\\\ \vdots\\\
\mathbf{\hat{\zeta}}^{C},\end{bmatrix}\in\mathbb{R}^{CLM},\quad
R=\begin{bmatrix}R_{1}&\ldots&R_{1,C}\\\ \vdots&\ddots&\vdots\\\
R_{C,1}&\ldots&R_{C}\end{bmatrix}.$ (22b)

Here the block $R_{i,j}\in\mathbb{R}^{LM,LM}$ encodes the cross-correlation
between the classifier $i$ and classifier $j$ for all inputs in the sequence.
This cross-correlation could, e.g., be caused by using the same data set to
train the different classifiers. It is, however, unclear how to compute this
cross-correlation between classifiers.

From aggregated state $\zeta$, the prediction for the different classifiers
can be recovered as

$\displaystyle\zeta^{(k)}\sim\mathcal{N}(\zeta;\hat{\zeta},R),\quad
k=1,2,\dots K$ (23a) $\displaystyle
z_{lc}^{(k)}=W_{lc}\zeta^{(k)},\quad\forall l,c$ (23b) where $W_{lc}$ is used
to find the classification of the $c$’th calasssifyer of the $l$’th input,
i.e.,
$W_{lc}=\begin{bmatrix}0&\ldots&0&I_{M}&0&\ldots&0\end{bmatrix}\in\mathbb{R}^{M,CLM}.$
(23c)

Then mc samples can be used to compute a point estimate of the pmf as

$\displaystyle\hat{p}_{lc}(y^{\star})=\frac{1}{K}\sum_{k=1}^{K}\text{softmax}(z_{lc})$
(24)

Afterward, the fusion formula (9c) can be applied to fuse prediction. However,
this is still a point estimate, and the fusion does not take that some
classifiers might have a larger covariance into consideration. Hence, two
strategies for aggregating the extended state of classifiers $\zeta$ that can
utilize the information that some classifiers might have a larger covariance
than others are presented in the following sections.

##  V Fusion of Neural Network Classifiers

It is possible to use the information from one prediction from an nn and fuse
it with a prediction from other sources of information (which might also be an
nn) having access to the uncertainty in the prediction for an nn. Given the
extended description in (22), for a set of classifications coming from $C$
classifiers and a sequence of $L$ inputs known to belong to the same class
$y^{\star}$, these predictions can be fused as follows

$\displaystyle P_{N}^{g}$ $\displaystyle=(H^{\top}R^{-1}H)^{-1},$ (25a)
$\displaystyle\hat{\bar{g}}_{N}$
$\displaystyle=P_{N}^{g}H^{\top}R^{-1}\hat{\zeta}$ (25b) where
$H=\begin{bmatrix}I_{M}\\\ \vdots\\\ I_{M}\end{bmatrix}\in\mathbb{R}^{CLM,M}.$
(25c)

Due to the Gaussian approximation in (22), the fused information will also be
Gaussian distributed. Hence, the distribution before the normalization can be
approximated as

$\displaystyle
p(z^{\star}|x^{\star}_{1:L})\approx\mathcal{N}\big{(}z;\bar{g}_{N},\bar{P}_{N}^{g}),$
(26)

where mc samples as in (12) can be applied to estimate the desired pmf. That
is, by changing (12a) to

$\displaystyle
z^{(k)}\sim\mathcal{N}\big{(}z;\hat{\bar{g}}_{N},P_{N}^{g}),\quad
k\\!=\\!1,2,\dots K.$ (27)

In this paper, two special cases of (22) are investigated:

  1. (i)

Multiple nn independently train and classify the same input, i.e., $L=1$.

  2. (ii)

An nn is used to classify a sequence of inputs known to belong to the same
class, i.e., $C=1$.

These are the edge cases, where a combination results in the general case. The
first case is motivated by safety-critical systems. There is often more than
one sensor measuring the environment. Multiple sensors measuring the same
thing are often used to add redundancy to the system.

For the second case, if one classification is uncertain, the additional
information that all the inputs in the sequence belong to the same class will
result in a more robust classification. In a standard object detection
algorithm, this setup can be used on cropped images of the tracked object.
After the object detection algorithm finds a bounding box locating an object,
it crops the image and sends it to a new classifier that classifies the
cropped image. Then, the sequential fusion can be used on the cropped images.
The second case of fusion strategies can be used in the case of test time
augmentation where a single nn is predicting the output from multiple
augmented versions of the same input [26].

Transform the output before the softmax step using (11).  Estimate
$p_{lc}(z^{\star})$ using lla in (18). Transform the output before the softmax
step using (11).  Estimate $p_{lc}(z^{\star})$ using lla in (18). Aggregate
the classifiers using (22). Estimate $p_{lc}(z^{\star})$ using lla in (18).
Transform the output before the softmax step using (11). Estimate
$p_{lc}(z^{\star})$ using lla in (18). Compress $\zeta$ using the fusion in
(25) or ella. Estimate the pmf by mc sampling using (12).
$x_{l},\hat{\theta}_{N}^{c}$$\bar{g}(x_{l}|\hat{\theta}^{c}_{N})$$p_{lc}(z^{\star})$$\zeta$$p(z^{\star}|x_{1:L}^{\star})$$\hat{p}(y^{\star}|x_{1:L}^{\star})$

Figure 1:  A block diagram of the proposed strategy to estimate
$\hat{p}(y^{\star}|x_{1:L}^{\star})$. Note that, two different methods are
presented in this paper to compress the aggregated classification to quantify
$p(z^{\star}|x_{1:L}^{\star})$. Namely, the fusion strategy proposed in
Section V or the ella proposed in Section VI.

##  VI Ensemble linearized Laplacian approximation

Fusing the prediction using (25) to estimate $p(z^{\star}|x^{\star}_{1:L})$
summarizes the information from multiple distributions with a single mode into
a new distribution with a single mode. The lla provides a local description of
the uncertainty in the prediction, from which efficient mc sampling can
estimate the pmf, marginalizing the uncertainty in the parameters. However,
the Laplacian approximation is a local approach requiring more flexibility to
represent multimodal distributions. Hence, if the loss of (14) is multimodal,
the approximation might be inaccurate. On the other hand, the ensemble method,
training multiple nns, e.g., deep ensemble [21], can represent multimodal
distributions. However, it is very costly to create a new sample, i.e., since
every new sample requires training of a new nn. Hence, we suggest combining
the two methods mentioned, i.e., training multiple nns and creating an lla for
each traned nn. This proposed extension of the lla is referred to as ella.
Instead of drawing samples from the distribution described by (26), one can
draw samples from the full description of the distribution in (23). Afterward,
(12) can be applied to estimate the desired pmf. That is, by changing (12a) to

$\displaystyle z^{(k)}=\sum_{l=1}^{L}\sum_{c=1}^{C}w_{lc}z_{lc}^{(k)}.$ (28)

That is, $p(z^{\star}|x^{\star}_{1:L})$ is approximated by $CL$ Gaussian
distributions where the parameter is used $w_{lc}$ to weigh the impact of the
different modes. For example, could $w_{lc}$ be chosen proportional to the
covariance $P_{N}^{g,lc}$ for the classifier. This results in a method that
can represent a multimodal distribution and is efficient to sample from. The
resulting method is similar to a Gaussian mixture model but with some key
differences [44]. Instead of independently sampling from the difference modes,
the sampling is done over the extended state $\zeta$, where the weighting is
done after the sampling. This enables us to include cross-correlation between
the modes.

Similar ideas to combine different methods to quantify the uncertainty have
been presented in, e.g., [45], where deep ensemble and mc dropout are
combined. Ideally, each realization should represent a different mode in the
distribution. To make this more likely, repulsive training can be used when
training the nns [46], where it is enforced during the training process that
the different nns converge to different parameters.

The complete strategy proposed in this paper to estimate
$\hat{p}(y^{\star}|x^{\star}_{1:L})$ is summarized in Fig. 1. Here, for the
compression of $\zeta$, either the fusion strategy in Section V or the ella
proposed in this section could be used.

##  VII Validation

Suppose access to some validation data set
$\mathcal{V}=\\{y^{\circ}_{n},x^{\circ}_{n}\\}_{n=1}^{N_{\circ}}$. The
question is how to validate that the estimated pmf resembles the true one. The
inherent difficulty is that the validation data, as the training data,
consists of inputs with corresponding class labels and not measures of a pmf.
This is one of the reasons why a unified qualitative evaluation metric is
lacking for the uncertainty in the prediction [15]. However, some of the most
common metrics used are classification accuracy, log-likelihood (ll),
prediction entropy, Brier score, expected calibration error (ece), area under
the receiver operating characteristic curve (auroc), and the area under the
precision-recall curve (aupr) when detecting out-of-distribution samples.

Negative ll, prediction entropy, and Brier score are all proper scoring rules,
i.e., they emphasize a careful and honest assessment of the uncertainty [47].
However, none of them measure the calibration, i.e., the reliability of the
estimated pmf. Hence, the most important metric is the ece when comparing the
method on in-distribution data. For out-of-distribution data, the auroc and
aupr are useful metrics, as well as the difference in predicted entropy for in
and out-of-distribution data.

###  VII-A Accuracy and calibration

Calculate the $J$-bin histogram defined as

$B_{j}=\bigg{\\{}n:\frac{j-1}{J}\leq\max_{m}\hat{p}(m|x^{\circ}_{n})<\frac{j}{J}\bigg{\\}}$
(29)

from the validation data. For a perfect classifier $B_{j}=\emptyset$ for
$j<J$. For a classifier that is just guessing, all sets are of equal size,
i.e., $|B_{j}|=|B_{i}|$ $\forall i,j$. Note that
$\max_{m}\hat{p}(m|x^{\circ}_{n})\geq 1/M$, so the first bins will be empty if
$J>M$.

The accuracy of the classifier is calculated by comparing the size of each set
with the actual classification performance within the set. That is,

$\text{acc}(B_{j})=\frac{1}{|B_{j}|}\sum_{n\in
B_{j}}\mathbb{1}\big{(}\hat{y}^{\circ}_{n}=y^{\circ}_{n}\big{)}$ (30a) where
$\hat{y}^{\circ}_{n}=\operatorname*{arg\,max}_{m}\hat{p}(m|x^{\circ}_{n})$
(30b)

Instead of certainty, from hereon, the standard and the equivalent notion of
confidence will be used [43, 48]. The mean confidence in a set is denoted
$\text{conf}(B_{j})$ and is defined as

$\displaystyle\text{conf}(B_{j})=\frac{1}{|B_{j}|}\sum_{n\in
B_{j}}\max_{m}\hat{p}(m|x^{\circ}_{n}),$ (31)

This measures how much the classifier trusts its estimated class labels. In
contrast to the accuracy, it does not depend on the annotated class labels
$y_{n}$. Comparing accuracy to confidence gives the ece, defined as

$\displaystyle\text{{ece}}=\sum_{j}^{J}\frac{1}{|B_{j}|}|\text{acc}(B_{j})-\text{conf}(B_{j})|.$
(32)

A small value indicates that the weight is a good measure of the actual
performance.

###  VII-B Log likelihood, predicted entropy, and Brier score

The total negative ll is given by (14), while the predicted entropy for one
input is given by

$E=-\sum_{m=1}^{M}\hat{p}(m|x^{\circ}_{n})\ln\bigg{(}\hat{p}(m|x^{\circ}_{n})\bigg{)}.$
(33)

Here $\hat{p}(m|x^{\circ}_{n})$ denotes some generic estimate of the pmf.
Hence, the difference between ll and entropy is that to calculate the
predicted entropy, the class label $y^{\circ}$ information is not required.
That means the predicted entropy can be computed for out-of-distribution
examples. For out-of-distribution detections, the difference in entropy
between the in-distribution data and out-of-distribution data is of interest,
i.e., $\Delta_{E}=E_{in}-E_{out}$.

The Brier score [47, 48] corresponds to the least squares fit

$\displaystyle\frac{1}{N_{\circ}}\sum_{n=1}^{N_{\circ}}\sum_{m=1}^{M}\bigl{(}\delta_{m,y_{n}^{\circ}}-\hat{p}(m|x^{\circ}_{n})\bigr{)}^{2},$
(34)

where $\delta_{i,j}$ denotes the Kronecker delta function.

Figure 2:  The images visualize the distribution for the fusion of prediction from multiple lla and the distribution generated by an ella. Here, the visualized using images from cfar10\.  TABLE I:  The computed performance measures for the two datasets. The arrows indicate whether a high or low value is preferable. For a description of measures, see Section VII.  | mnist
---|---
Method | acc. $\uparrow$ |  ll $(10^{3})$ $\uparrow$ | Brier score $\downarrow$ |  ece $\downarrow$
Temp. sc. [43] | 92$\%$ | 7.83 | 0.123 | 2.85
Deep ensemble [21] | 95$\%$ | 7.94 | 0.080 | 2.29
lla [29] | 92$\%$ | 7.81 | 0.125 | 1.33
Fusion of llas, Section V | 94$\%$ | 7.89 | 0.103 | 1.82
ella $w_{lc}=1/LC$, Section VI | 93$\%$ | 8.37 | 0.071 | 1.06
|  cfar10
Method | acc. $\uparrow$ |  ll $(10^{3})$ $\uparrow$ | Brier score $\downarrow$ |  ece $\downarrow$
Temp. sc. [43] | 82$\%$ | 7.45 | 0.246 | 1.42
Deep ensemble [21] | 86$\%$ | 7.87 | 0.197 | 1.19
lla [29] | 82$\%$ | 7.74 | 0.255 | 1.08
Fusion of llas, Section V | 86$\%$ | 8.04 | 0.197 | 0.67
ella $w_{lc}=1/LC$, Section VI | 87$\%$ | 8.20 | 0.191 | 0.63

###  VII-C Receiver operating characteristic, precision, and recall

The different methods of quantified uncertainty in the prediction can be
evaluated on how well out-of-distribution samples can be detected. How well
can the prediction of an input that the model has been trained on, in-
distribution, be distinguished from an entirely different input, i.e., out-of-
distribution examples. Three widely used measures to measure how well out-of-
distribution samples are detected are probability of detection, $P_{D}$ (also
known as sensitivity or recall $P_{R}$) probability of false alarm $P_{FA}$
and precision $P_{P}$ [49, 50]. Assume that the $\mathcal{H}_{0}$ hypothesis
is that the samples are out-of-distribution samples and that under the
$\mathcal{H}_{1}$ hypothesis, the samples are in-distribution samples. Then,
defining the three quantities as

$\displaystyle P_{D}=P_{R}$ $\displaystyle=Pr\\{\text{detect
}\mathcal{H}_{1}|\mathcal{H}_{1}\\},$ (35a) $\displaystyle P_{FA}$
$\displaystyle=Pr\\{\text{detect }\mathcal{H}_{1}|\mathcal{H}_{0}\\},$ (35b)
$\displaystyle P_{P}$ $\displaystyle=Pr\\{\mathcal{H}_{1}|\text{detect
}\mathcal{H}_{1}\\}.$ (35c)

The three measures in (35) will change depending on a set threshold. The
receiver operating characteristic (roc) curve is created by plotting the
$P_{D}$ against the $P_{FA}$ for different threshold values. Similarly, the
precision-recall curve is created by plotting the precision against the recall
for different threshold values. To evaluate how well different methods are at
detecting out-of-distribution samples from in-distribution samples, one often
refers to auroc and aupr. The auroc and the aupr are between zero and one, and
the value is one for a perfect detector.

###  VII-D Temperature scaling

The output of the softmax function can be used as a point estimate of pmf.
However, this does not consider the uncertainty in the parameters, and this
approach is known not accurately to reflect the uncertainty in the prediction,
i.e., the uncertainty needs to be calibrated [43, 48]. It is often
overconfident and underestimates the uncertainty. Hence assigning small
probability mass to unlikely classes.

A common approach used to calibrate the estimated pmf is the so-called
temperature scaling [43]. In temperature scaling, $g(x|\theta)$ is scaled by a
scalar quantity $T$ before the normalization by the softmax operator. With a
slight abuse of notation, introduce

$\hat{p}^{T}_{lc}(y^{\star})\triangleq
f(x^{\star}_{l}|\hat{\theta}_{N}^{c},T)=\text{softmax}\left(g(x^{\star}|\hat{\theta}_{N}^{c})/T\right).$
(36)

The scaling parameter $T$ is found after training on the model using
validation data. The subscript $T$ is used for the estimate of the pmf when
temperature scaling is used.

(a) Temp. sc.

(b) Deep ensemble.

(c) lla.

(d) Fusion of llas.

(e) ella.

Figure 3:  Distribution for the predicted entropy (33) for in-distribution (mnist) and out-of-distribution (fmnist) samples. Here, five different methods to create the estimated conditional pmf are investigated. TABLE II:  Performance measure for out-of-distribution detection. The arrows indicate whether a high or low value is preferable. See Section VII for descriptions of the measures.  |  mnist (in) fmnist (out)
---|---
Method |  $\Delta_{E}$ $(10^{3})$ $\downarrow$ |  auroc $\uparrow$ |  aupr $\uparrow$
Temp. sc. [43] | -6.36 | 0.82 | 0.70
Deep ensemble [21] | -13.12 | 0.94 | 0.91
lla [29] | -11.23 | 0.91 | 0.91
Fusion of llas, Section V | -13.83 | 0.95 | 0.96
ella $w_{lc}=1/LC$, Section VI | -5.42 | 0.87 | 0.83
Figure 4:  The image illustrates the fusion of a sequence of images from the
cfar10 dataset. Here, using the same nn to predict the class for all the
images in the sequence and using (25) to fuse the prediction. Hence, the
cross-correlation between the predictions is taken into consideration.

##  VIII Experiment study

The experiments in this section can be divided into two parts. The first part
investigates the fusion of multiple classifiers that classify the output for
the same input. Here, an ensemble of nns are trained on the two classical data
sets mnist [51] and cfar10 [52]. The nn trained on mnist will be evaluated on
Fashion mnist (fmnist) [53] to investigate how to detect out-of-distribution
examples. This part investigates the fusion strategy described in Section V
and ella described in Section VI. In the second part, we are simulating the
setup where a sequence of images known to belong to the same class by giving a
classifier a sequence of images from the cfar10 data set and data from a
camera trap that captures images of Swedish carnivores.

For the three different classification tasks, i.e., mnist, cfar10, and camera
trap images of Swedish carnivores, different structures of the nns were used.
For the mnist task, a five-layer fully connected nn was used; for the cfar10
task, using a LeNet5 [54] inspired structure; and for the Swedish carnivores,
using an AlexNet [55] inspired structure. Here, using the adam optimizer [56]
with the standard setting to estimate parameters $\theta$ for all the tasks.
The training was done for three epochs for the mnist task, ten epochs for the
cfar10, and six epochs for the Swedish carnivores.

This paper is investigating five different methods to quantify the
uncertainty, namely

  1. (i)

Temperature scaling [43], see Section VII-D.

  2. (ii)

Deep ensemble [21]

  3. (iii)

One single lla, [29].

  4. (iv)

Fusion of the prediction of an ensemble of lla, see Section V.

  5. (v)

An ella where every mode is equally weighted, i.e., (28) with $w_{lc}=1/LC$.

For the deep ensemble, the size of the ensemble was ten for the mnist task and
five for the cfar10 task. This was also the number of nns trained when fusing
the prediction from multiple llas and the number of modes in the ella.

###  VIII-A Fusion of multiple measurements

Having access to predictions from multiple independently trained nns, each
associated with an lla, it is possible to fuse their prediction to better
understand the predicted output of a given input. Table I shows the result of
the experiment done on mnist and cfar10\. For the cfar10 data set, all the
measures are improving for the fusion of llas compared to only using the
temperature scaling, deep ensemble, and a single lla. In particular, the ece
decreases. However, for the mnist data set, the ece does not decrease for the
fusion of llas compared to only using a single lla. That the fusion of llas
does not decrease the ece can be a consequence of that some information about
the distribution goes missing in the fusion stage. An example of the fusion of
the prediction can be seen in the second column in Fig. 2. Here, the mean and
standard deviation for the different modes are shown in red, and the fused
prediction is in blue. The fourth column shows the estimated pmf for the fused
prediction.

In Fig. 3, the distribution for the predicted entropy for images from both
images the nn is trained on (in-distribution) and images the nn is not trained
on (out-of-distribution). The in distribution is mnist, i.e., images of
numbers, and the out of distribution is fmnist, i.e., images of clothing
items. Compared to the other methods, fusing the prediction from multiple llas
seen in Fig. 3(d) provides a more evident visual difference between the in and
out distribution images. This claim is also supported in Table II where the
auroc and aupr are better for the fused predictions. The difference in the
predicted entropy for in and out of distribution samples for the fused
predictions is also significant; see Table II.

###  VIII-B Sequential fusion

In Fig. 4, a sequence of six images of aircraft from cfar10 is shown. Here,
classifying all the images in the sequence by the same nn. Hence, the
predictions are correlated according to (25). Assuming that all the images in
the sequence belong to the same class, this information can be used to improve
the prediction accuracy from images later in the sequence. For example, even
though the prediction of the fourth image is uncertain, thanks to the
information from the previous images, the fused prediction indicates that it
is likely that the object in image is an aircraft.

Fig. 5 illustrates an experiment where the prediction of multiple images from
a camera trap. Here, a lynx is passing the camera, and even though the last
image in the sequence is classified as a wolf, by using the information from
the previous images, the fused prediction still indicates that it is a lynx in
the image.

###  VIII-C Multimodal model

It is possible to have a multimodal representation of the prediction
distribution with access to an ella. To visualize the multimodality, one can
draw independent samples from the different modes created by the different nn
in the ensemble. Fig. 2 is visualizing this for two images from the cfar10
dataset. The third column shows a representation of the distribution before
the normalization. Here, independent mc-samples are taken from the different
modes to visualize the distribution.

In Table I, it can be seen that inclusion of the multimodality improves
calibration of the predicted uncertainty, i.e., leads to lower a ece. This
while still maintaining high accuracy and high ll. The use of the ella to
detect out-of-distribution samples can be seen in Fig. 3(e) and Table II.
Here, it is shown that multimodality improves the detection compared to
temperature scaling, but the improvements could be better when fusing the
prediction from multiple llas.

Figure 5:  The figure illustrates the fusion of the predicted class from an
image sequence of a lynx. Here, using the same nn to predict the object’s
class in all the images in the sequence, where (25) is used to fuse the
prediction over the sequence.

##  IX Summary and conclusion

This paper suggests two methods to estimate the conditional pmf given a
sequence of classifications when one can estimate the uncertainty in the
prediction from an nn. The first method is the fusion of the predictions,
where both the information from the previous and the cross-correlation between
the predictions are considered. The suggested method to include uncertainty in
the prediction from an nn used in the fusion strategy is called lla. The
second method extends the lla to represent multimodality in the predicted
distribution. This extension is referred to as ella, where it combines a deep
ensemble and lla.

The fusion strategies and the proposed extension are evaluated on
classification tasks with images from the classical data sets mnist and cfar10
and images from camera traps collected in the Swedish forests. It also
investigates how to use the methods to detect out-of-distribution samples.
Here, mnist is used as in-distribution and fmnist as out-of-distribution
samples. Sequential fusion of the predictions is shown to increase the
robustness of the prediction. Both when all predictions in the sequence are
uncertain and only one prediction is confident but incorrect. The results in
terms of ece show that the multimodal representation using the extended ella
method improves the performance compared to previous work. In most cases, the
fusion of multiple predictions improved the performance. Furthermore, when
detecting out-of-distribution samples, the fusion of predictions from multiple
nn leads to improvement in performance. This improvement is in terms of
increased auroc and aupr and the difference in predicted entropy. Here, the
performance using ella could be improved. This indicates that some modes might
be more critical in detecting out-of-distribution samples. Hence, methods to
give the modes different weights are a possible direction for future research.

## References

  * [1] A. Ramesh _et al._ , “Zero-shot text-to-image generation,” in _Proc. of the 37th Int. Conf. on Mach. Learn. (ICML)._ Online: PMLR, 23–29 Jun 2021, pp. 8821–8831.
  * [2] D. Silver _et al._ , “Mastering the game of go with deep neural networks and tree search,” in _Nature_ , vol. 529. Nature Publishing Group, 2016, pp. 484–489.
  * [3] Q. Li _et al._ , “Deep neural networks for improved, impromptu trajectory tracking of quadrotors,” in _Proc. of IEEE Int. Conf. on Robot. and Autom. (ICRA)_. Singapore, Singapore: IEEE, 2017, pp. 5183–5189, 19 May–3 June.
  * [4] X. Zhang, M. Bujarbaruah, and F. Borrelli, “Safe and near-optimal policy learning for model predictive control using primal-dual neural networks,” in _Proc. of IEEE American Contr. Conf. (ACC)_. Philadelphia, PA, USA: IEEE, 2019, pp. 354–359, 10-12 July.
  * [5] A. Paleyes, R.-G. Urma, and N. D. Lawrence, “Challenges in deploying machine learning: a survey of case studies,” in _Adv. in Neural Inf. Process. Syst. (NeurIPS) 34 Workshop: ML Retrospectives, Surveys & Meta-Analyses (ML-RSA)_, vol. 33, Virtual, 7–12 Dec 2020.
  * [6] S. A. Bagloee, M. Tavana, M. Asadi, and T. Oliver, “Autonomous vehicles: challenges, opportunities, and future implications for transportation policies,” in _J. of Modern Trans._ , vol. 24, no. 4. Springer, Dec. 2016, pp. 284–303.
  * [7] S. Grigorescu, B. Trasnea, T. Cocias, and G. Macesanu, “A survey of deep learning techniques for autonomous driving,” in _J. of Field Robotics_ , vol. 37, no. 3. Wiley Online Library, Jan. 2020, pp. 362–386.
  * [8] NTSB, “Preliminary Report Highway HWY18MH010,” National Transportation Safty Board (NTSB), Technical Specification (TS), 05 2018. [Online]. Available: https://data.ntsb.gov/Docket/?NTSBNumber=HWY18MH010
  * [9] A. D’Amour _et al._ , “Underspecification presents challenges for credibility in modern machine learning,” in _arXiv preprint arXiv:2011.03395_ , 2020.
  * [10] Z. Ghahramani, “Probabilistic machine learning and artificial intelligence.” in _Nature_ , vol. 521, no. 7553, Feb. 2015, pp. 452 – 459.
  * [11] K. Patel and S. Waslander, “Accurate prediction and uncertainty estimation using decoupled prediction interval networks,” in _arXiv preprint arXiv:2202.09664_ , 2022.
  * [12] Y. Ovadia _et al._ , “Can you trust your model’s uncertainty? evaluating predictive uncertainty under dataset shift,” in _Adv. in Neural Inf. Process. Syst. (NeurIPS) 33_ , Vancouver, Canada, 2019, 8–14 Dec.
  * [13] S. Lin, R. Clark, N. Trigoni, and S. Roberts, “Uncertainty estimation with a VAE-classifier hybrid model,” in _Proc. of IEEE Int. Conf. on Acoust., Speech and Signal Processing (ICASSP)_. Singapore, Singapore,: IEEE, 2022, pp. 3548–3552, 22–17 May.
  * [14] P. Izmailov, P. Nicholson, S. Lotfi, and A. G. Wilson, “Dangers of Bayesian model averaging under covariate shift,” in _Adv. in Neural Inf. Process. Syst. (NeurIPS) 35_ , vol. 34, 6–14 Dec 2021, virtual.
  * [15] J. Gawlikowski _et al._ , “A survey of uncertainty in deep neural networks,” in _arXiv preprint arXiv:2107.03342_ , 2021.
  * [16] B. Charpentier, D. Zügner, and S. Günnemann, “Posterior network: Uncertainty estimation without ood samples via density-based pseudo-counts,” vol. 33, Virtual, 07–12 Dec 2020.
  * [17] A. Kendall and R. Cipolla, “Modelling uncertainty in deep learning for camera relocalization,” in _Proc. of IEEE Int. Conf. on Robot. and Autom. (ICRA)_ , Stockholm, Sweden, 2016, pp. 4762–4769, 16–21, May.
  * [18] F. K. Gustafsson, M. Danelljan, G. Bhat, and T. B. Schön, “Energy-Based Models for Deep Probabilistic Regression,” in _Proc. of 16th European Conf. on Comput. Vision (ECCV)_ , Glasgow, UK/Online, 2020, pp. 325–343, 23-28 Aug.
  * [19] A. Eldesokey, M. Felsberg, K. Holmquist, and M. Persson, “Uncertainty-aware CNNs for depth completion: Uncertainty from beginning to end,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 12 014–12 023.
  * [20] C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra, “Weight Uncertainty in Neural Networks,” in _Proc. of the 32nd Int. Conf. on Mach. Learn. (ICML)._ , Lille, France, 2015, pp. 1613–1622, 6–11 Jul.
  * [21] B. Lakshminarayanan, A. Pritzel, and C. Blundell, “Simple and scalable predictive uncertainty estimation using deep ensembles,” in _Adv. in Neural Inf. Process. Syst. (NeurIPS) 31_. Curran Associates, Inc., 4–9 Dec 2017, Long Beach, CA, USA, 4–9 Dec.
  * [22] Y. Gal and Z. Ghahramani, “Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning,” in _Proc. of the 33td Int. Conf. on Mach. Learn. (ICML)._ , New York, NY, USA, 2016, pp. 1050–1059, 20–22 Jun.
  * [23] M. Teye, H. Azizpour, and K. Smith, “Bayesian Uncertainty Estimation for Batch Normalized Deep Networks,” in _Proc. of the 35th Int. Conf. on Mach. Learn. (ICML)._ , 6–11 Jul 2018, pp. 4907–4916, Stockholm, Sweden ,6–11 Jul.
  * [24] W. J. Maddox _et al._ , “A simple baseline for Bayesian uncertainty in deep learning,” in _Adv. in Neural Inf. Process. Syst. (NeurIPS) 33_ , Vancouver, Canada, 2019, 8–14 Dec.
  * [25] K. Osawa _et al._ , “Practical deep learning with Bayesian principles,” in _Adv. in Neural Inf. Process. Syst. (NeurIPS) 33_ , Vancouver, Canada, 2019, 8–14 Dec.
  * [26] M. S. Ayhan and P. Berens, “Test-time data augmentation for estimation of heteroscedastic aleatoric uncertainty in deep neural networks,” in _1st Conf. on Medical Imaging with Deep Learn. (MIDL)_ , Amsterdam, The Netherlands, 2018, 4–6 Jul.
  * [27] E. Ilg _et al._ , “Uncertainty estimates and multi-hypotheses networks for optical flow,” in _Proc. of 15th European Conf. on Comput. Vision (ECCV)_ , Munich, Germany, 2018, pp. 652–667, 8-14 Sep.
  * [28] G. Carannante, N. C. Bouaynaya, and L. Mihaylova, “An enhanced particle filter for uncertainty quantification in neural networks,” in _Proc. of IEEE 24th Int. Conf. on Inf. Fusion (FUSION)_. Sun City, South Africa/ Virtual: IEEE, 2021, pp. 1–7, Nov 1-4.
  * [29] M. Malmström, I. Skog, D. Axehill, and F. Gustafsson, “Uncertainty quantification in neural network classifiers–a local linear approach,” _arXiv preprint arXiv:2303.07114_ , 2023.
  * [30] A. Immer, M. Korzepa, and M. Bauer, “Improving predictions of Bayesian neural nets via local linearization,” in _Proc. of 24nd Int. Conf. on Artificial Intell. and Statistics. (AISTATS)_. San Diego, CA, USA: PMLR, 2021, pp. 703–711, 13-15 Apr.
  * [31] C. M. Bishop, _Pattern Recognition and Machine Learning_. Springer, 2006, new York, NY, USA.
  * [32] J. Martens and R. Grosse, “Optimizing neural networks with Kronecker-factored approximate curvature,” in _Proc. of the 32nd Int. Conf. on Mach. Learn. (ICML)._ , Lille, France, 2015, 6–11 Jul.
  * [33] A. Kristiadi, M. Hein, and P. Hennig, “Being Bayesian, even just a bit, fixes overconfidence in ReLU networks,” in _Proc. of the 37th Int. Conf. on Mach. Learn. (ICML)._ , Online, 2020, 13-18 July.
  * [34] D. J. MacKay, “A practical Bayesian framework for backpropagation networks,” in _Neural computation_ , vol. 4, no. 3. MIT Press, 1992, pp. 448–472.
  * [35] Z. Deng, F. Zhou, and J. Zhu, “Accelerated linearized Laplace approximation for Bayesian deep learning,” _arXiv preprint arXiv:2210.12642_ , 2022.
  * [36] H. Liero and S. Zwanzig, _Introduction to the theory of statistical inference_. Chapman and Hall CRC Texts in Statistical Science, Boca Raton, FL, USA, 2011.
  * [37] J. T. G. Hwang and A. A. Ding, “Prediction Intervals for Artificial Neural Networks,” in _J. Am. Stat. Assoc. (JSTOR)_ , vol. 92. Taylor & Francis, 1997, pp. 748–757.
  * [38] M. Malmström, “Uncertainties in neural networks a system identification approach,” Licentiate Thesis, Dept. Elect. Eng., Linköping University, Linköping, Sweden, Apr 2021.
  * [39] G. Koliander, Y. El-Laham, P. M. Djurić, and F. Hlawatsch, “Fusion of probability density functions,” _Proceedings of the IEEE_ , vol. 110, no. 4, pp. 404–453, 2022.
  * [40] A. Lindholm, N. Wahlström, F. Lindsten, and T. B. Schön, _Machine Learning: A First Course for Engineers and Scientists_. Cambridge University Press, 2022.
  * [41] I. Johnstone, “High dimensional bernstein-von mises: simple examples,” _Institute of Mathematical Statistics collections_ , vol. 6, pp. 87–98, 01 2010.
  * [42] M. Malmström, I. Skog, D. Axehill, and F. Gustafsson, “Detection of outliers in classification by using quantified uncertainty in neural networks,” in _Proc. of IEEE 25th Int. Conf. on Inf. Fusion (FUSION)_. Linköping, Sweden: IEEE, 2022, Jul 4-7.
  * [43] C. Guo, G. Pleiss, Y. Sun, and K. Q. Weinberger, “On calibration of modern neural networks,” in _Proc. of the 34th Int. Conf. on Mach. Learn. (ICML)._ Sydney, Australia: PMLR, 2017, pp. 1321–1330, 06–11 Aug.
  * [44] C. M. Bishop, _Pattern Recognition and Machine Learning_. Springer, 2006, new York, NY, USA, pages 111-112.
  * [45] R. Pop and P. Fulop, “Deep ensemble Bayesian active learning: Addressing the mode collapse issue in monte carlo dropout via ensembles,” _arXiv preprint arXiv:1811.03897_ , 2018.
  * [46] F. D’Angelo and V. Fortuin, “Repulsive deep ensembles are Bayesian,” in _Adv. in Neural Inf. Process. Syst. (NeurIPS) 35_ , vol. 34, 6–14 Dec 2021, pp. 3451–3465, virtual.
  * [47] T. Gneiting and A. E. Raftery, “Strictly proper scoring rules, prediction, and estimation,” in _J. of the American Statistic. Association_ , vol. 102, no. 477. Taylor & Francis, 2007, pp. 359–378.
  * [48] J. Vaicenavicius _et al._ , “Evaluating model calibration in classification,” in _Proc. of 22nd Int. Conf. on Artificial Intell. and Statistics. (AISTATS)_. Naha, Okinawa, Japan: PMLR, 2019, pp. 3459–3467, 16-18 Apr.
  * [49] S. M. Kay, _Fundamentals of statistical signal processing: Detection theory_. Prentice-Hall signal processing series, 1998, ch. 7, Upper Saddle River, NJ, USA.
  * [50] I. Goodfellow, Y. Bengio, and A. Courville, _Deep learning_. MIT Press, 2016, London, England.
  * [51] Y. LeCun, C. Cortes, and C. J. Burges, “The MNIST database of handwritten digits,” 1998. [Online]. Available: http://yann.lecun.com/exdb/mnist/
  * [52] A. Krizhevsky, “Learning multiple layers of features from tiny images,” Master’s thesis, University of Toronto, Canada, Apr. 2009.
  * [53] H. Xiao, K. Rasul, and R. Vollgraf, “Fashion-MNIST: a novel image dataset for benchmarking machine learning algorithms,” _arXiv preprint arXiv:1708.07747_ , 2017.
  * [54] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” _Proc. of the IEEE_ , vol. 86, no. 11, pp. 2278–2324, 1998.
  * [55] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Adv. in Neural Inf. Process. Syst. (NeurIPS) 25_ , Lake Tahoe, NVn USA, 2012, pp. 1097–1105, 3-8, Dec.
  * [56] D. P. Kingma and J. Ba, “ADAM: A Method for Stochastic Optimization.” in _Proc. of 3rd Int. Conf. for Learn. Representations (ICLR)_ , 7–9 May 2015, San Diego, CA, USA.