Fusion framework and multimodality for the Laplacian approximation of
  Bayesian neural networks
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.
Magnus Malmström, Isaac Skog, Daniel Axehill, Fredrik Gustafsson
# Fusion framework and multimodality for the Laplacian approximation of Bayesian neural networks

###### 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 nn (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 cfari10, 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, nn 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 nn in safety-critical applications, e.g., medical imaging and autonomous driving [5, 6, 7]. One of the main reasons nn 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, 31, 32, 33, 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

\[\mathcal{T}^{(c)}=\{x_{i},y_{i}\}_{i=1}^{N_{c}},\quad(x_{i},y_{i })\stackrel{{\text{\tiny def}}}{{\sim}}p(X,Y),\quad c=1,2,\ldots,C. \tag{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

\[p(Y=y|X=x)\approx f(x|\mathcal{T}^{(c)}),\quad c=1,2,\ldots,C. \tag{2}\]

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

\[\hat{p}_{lc}(y^{\star})=f(x_{l}^{\star}|\mathcal{T}^{(c)}),\quad l =1,\ldots L,\ c=1,\ldots,C. \tag{3}\]

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

\[p(y^{\star}|x_{1:L}^{\star})\triangleq p(y^{\star}|x_{1:L}^{\star},\mathcal{T }^{(1)},\ldots,\mathcal{T}^{(C)}). \tag{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.

### _Fusion using non-parametric models_

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

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

Thus, one reasonable way to fuse the pmf estimates is

\[\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)}. \tag{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

\[\hat{p}(y|x_{1:L}^{\star})\propto\prod_{l=1}^{L}\prod_{c=1}^{C} \hat{p}_{lc}^{w_{lc}}(y). \tag{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.

### _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_{l}^{\star}|\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}_{N}^{c}\) from the training data \(\mathcal{T}^{(c)}\). This dependence will be made explicit in the sequel, so we denote the conditional probabilities

\[\hat{p}_{lc}(y^{\star})\triangleq f(x_{l}^{\star}|\hat{\theta}_{N}^{c}), \quad\forall l,c. \tag{8}\]

The uncertainty in \(\hat{\theta}_{N}^{c}\) implies that the conditional probability \(f(x_{l}^{\star}|\hat{\theta}_{N}^{c})\) 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}_{N}^{c}\) with Monte Carlo (mc) samples \(\theta^{c(k)}\). That is, we use the approximation

\[\theta^{c(k)} \sim p(\theta|\mathcal{T}^{(c)}),\quad k=1,2\ldots,K, \tag{9a}\] \[\hat{p}_{lc}(y^{\star}) \approx\frac{1}{K}\sum_{k=1}^{K}f(x_{l}^{\star}|\theta^{c(k)}), \quad\forall l,c,\] (9b) \[\hat{p}(y^{\star}|x_{1:L}^{\star}) \propto\prod_{l=1}^{L}\prod_{c=1}^{C}\hat{p}_{lc}(y^{\star}). \tag{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_{1:L}^{\star})\) 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 [36, 38]. 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 uncertainty 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_{1:L}^{\star})\). 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.

### _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}. \tag{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. \tag{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.

### _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_{1:L}^{\star})\) 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:

\[z^{(k)} \sim p(z^{\star}|x_{1:L}^{\star}),\quad k=1,2,\ldots,K, \tag{12a}\] \[\hat{p}^{(k)} =\text{softmax}\big{(}z^{(k)}\big{)},\] (12b) \[\hat{p}(y^{\star}|x_{1:L}^{\star}) =\frac{1}{K}\sum_{k=1}^{K}\hat{p}^{(k)}. \tag{12c}\]

The uncertainty of the estimated class labels \(\hat{p}(y^{\star}|x_{1:L}^{\star})\) 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_{1:L}^{\star})\) 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_{1:L}^{\star})\), 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_{1:L}^{\star})\).

### _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}^{\star}\). 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]

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

Here

\[L_{N}(\theta)=\sum_{n=1}^{N}\ln f_{y_{n}}(x_{n}|\theta) \tag{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), \tag{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{)}^{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
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}. \tag{23c}\]

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

\[\hat{p}_{lc}(y^{\star})=\frac{1}{K}\sum_{k=1}^{K}\text{softmax}(z_{lc}) \tag{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

\[P_{N}^{g} =(H^{\top}R^{-1}H)^{-1}, \tag{25a}\] \[\tilde{g}_{N} =P_{N}^{g}H^{\top}R^{-1}\hat{\zeta} \tag{25b}\]

where

\[H=\begin{bmatrix}I_{M}\\ \vdots\\ I_{M}\end{bmatrix}\in\mathbb{R}^{CLM,M}. \tag{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

\[p(z^{\star}|x_{1:L}^{\star})\approx\mathcal{N}\big{(}z;\bar{g}_{N},\bar{P}_{N }^{g}), \tag{26}\]

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

\[z^{(k)}\sim\mathcal{N}\big{(}z;\hat{\bar{g}}_{N},P_{N}^{g}),\quad k\!=\!1,2, \ldots K. \tag{27}\]

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

1. Multiple nn independently train and classify the same input, i.e., \(L=1\).
2. 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].

## VI Ensemble linearized Laplacian approximation

Fusing the prediction using (25) to estimate \(p(z^{\star}|x_{1:L}^{\star})\) 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

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

That is, \(p(z^{\star}|x_{1:L}^{\star})\) 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_{1:L}^{\star})\) 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_{n}^{\circ},x_{n}^{\circ}\}_{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 ce 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.

### _Accuracy and calibration_

Calculate the \(J\)-bin histogram defined as

\[B_{j}=\left\{n:\frac{j-1}{J}\leq\max_{m}\hat{p}(m|x_{n}^{\circ})<\frac{j}{J}\right\} \tag{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_{n}^{\circ})\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}}1\big{(}\hat{y}_{n}^{\circ }=y_{n}^{\circ}\big{)}\] (30a) where \[\hat{y}_{n}^{\circ}=\operatorname*{arg\,max}_{m}\hat{p}(m|x_{n}^{\circ}) \tag{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

\[\text{conf}(B_{j})=\frac{1}{|B_{j}|}\sum_{n\in B_{j}}\max_{m}\hat{p}(m|x_{n}^{ \circ}), \tag{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 ce, defined as

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

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

### _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_{n}^{\circ})\ln\bigg{(}\hat{p}(m|x_{n}^{\circ}) \bigg{)}. \tag{33}\]

Here \(\hat{p}(m|x_{n}^{\circ})\) 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

\[\frac{1}{N_{\circ}}\sum_{n=1}^{N_{\circ}}\sum_{m=1}^{M}\big{(}\delta_{m,y_{n} ^{\circ}}-\hat{p}(m|x_{n}^{\circ})\big{)}^{2}, \tag{34}\]

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

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

Fig. 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.
### _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

\[P_{D}=P_{R} =Pr\{\text{detect }\mathcal{H}_{1}|\mathcal{H}_{1}\}, \tag{35a}\] \[P_{FA} =Pr\{\text{detect }\mathcal{H}_{1}|\mathcal{H}_{0}\},\] (35b) \[P_{P} =Pr\{\mathcal{H}_{1}|\text{detect }\mathcal{H}_{1}\}. \tag{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 aurp. The auroc and the aurp are between zero and one, and the value is one for a perfect detector.

### _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}_{lc}^{T}(y^{\star})\triangleq f(x_{l}^{\star}|\hat{\theta}_{N}^{c},T )=\text{softmax}\left(g(x^{\star}|\hat{\theta}_{N}^{c})/T\right). \tag{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.

## 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 cfari10[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

\begin{table}
\begin{tabular}{l|c c c} \hline \hline  & \multicolumn{3}{c}{mnist (in) fmnist (out)} \\ \cline{2-4} Method & \(\Delta_{E}\) (\(10^{3}\)) \(\downarrow\) & auroc \(\uparrow\) & aurp \(\uparrow\) \\ \hline 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 lla, Section V & -13.83 & **0.95** & **0.96** \\ ella \(w_{lc}=1/LC\), Section VI & -5.42 & 0.87 & 0.83 \\ \hline \hline \end{tabular}
\end{table} 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.

\begin{table}
\begin{tabular}{l|c c c c} \hline \hline \multirow{2}{*}{Method} & \multicolumn{3}{c}{mnist} \\ \cline{2-4}  & acc. \(\uparrow\) & lL (\(10^{3}\)) \(\uparrow\) & Brier score \(\downarrow\) & ecE \(\downarrow\) \\ \hline 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 lla, Section V & 94\(\%\) & 7.89 & 0.103 & 1.82 \\ ella \(w_{lc}=1/LC\), Section VI & -93\(\%\) & **8.37** & **0.071** & **1.06** \\ \hline \hline \end{tabular}
\end{table} 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.

Fig. 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.
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. Temperature scaling [43], see Section VII-D.
2. Deep ensemble [21]
3. One single lla, [29].
4. Fusion of the prediction of an ensemble of lla, see Section V.
5. 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.

### _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 ecc decreases. However, for the mnist data set, the ecc does not decrease for the fusion of llas compared to only using a single lla. That the fusion of llas does not decrease the ecc 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.

### _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

Fig. 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.
images, the fused prediction still indicates that it is a lynx in the image.

### _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.

## 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 llla. The second method extends the llla to represent multimodality in the predicted distribution. This extension is referred to as ella, where it combines a deep ensemble and llla.

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.
