Probabilistic neural networks for improved analyses with
  phenomenological models
Physics models typically contain adjustable parameters to reproduce measured
data. While some parameters correspond directly to measured features in the
data, others are unobservable. These unobservables can, in some cases, cause
ambiguities in the extraction of observables from measured data, or lead to
questions on the physical interpretation of fits that require these extra
parameters. We propose a method based on deep learning to extract values of
observables directly from the data without the need for unobservables. The key
to our approach is to label the training data for the deep learning model with
only the observables. After training, the deep learning model can determine the
values of observables from measured data with no ambiguities arising from
unobservables. We demonstrate this method on the phenomenological R-matrix that
is widely utilized in nuclear physics to extract resonance parameters from
cross section data. Our deep learning model based on Transformers successfully
predicts nuclear properties from measurements with no need for the channel
radius and background pole parameters required in traditional R-matrix
analyses. Details and limitations of this method, which may be useful for
studies of a wide range of phenomena, are discussed.
C. H. Kim, K. Y. Chae, M. S. Smith, D. W. Bardayan, C. R. Brune, R. J. deBoer, D. Lu, D. Odell
# Eliminating unobservable parameters in physics models using deep learning

###### Abstract

Physics models typically contain adjustable parameters to reproduce measured data. While some parameters correspond directly to measured features in the data, others are unobservable. These unobservables can, in some cases, cause ambiguities in the extraction of observables from measured data, or lead to questions on the physical interpretation of fits that require these extra parameters. We propose a method based on deep learning to extract values of observables directly from the data without the need for unobservables. The key to our approach is to label the training data with only the observables. After training, the deep learning model can determine the values of observables from measured data with no ambiguities arising from unobservables. We demonstrate this method on the phenomenological _R_-matrix that is widely utilized in nuclear physics to extract resonance parameters from cross section data. Our deep learning model based on Transformers successfully predicts nuclear properties from measurement with no need for the channel radius and background pole parameters required in traditional _R_-matrix analyses. Details and limitations of this method, which may be useful for studies of a wide range of phenomena, are discussed.

Models are routinely used to understand physical systems by fitting their predictions to measurement data and extracting "best fit" values for the model parameters. In some cases, these parameters correspond directly to system observables, i.e., features that can be measured. In other cases, the parameters do not correspond to measurable features, but rather may codify approximations made in the model or may be part of the framework of the model itself. The presence of these unobservable parameters can often lead to ambiguities in the extraction of observables from measured data or lead to questions on the physical interpretation of fits that require these extra parameters [1; 2; 3]. When possible, it is advantageous to extract values of observable parameters directly from features in the data without the need for unobservables.

One of the strengths of deep learning is the ability to identify features in data sets. For this reason and others, deep learning has been successfully used in a variety of science fields [4; 5; 6; 7; 8; 9; 10; 11; 12]. In this study, we propose to use deep learning to extract values of observable parameters directly from measured data sets, bypassing the need for unobservable parameters that can introduce ambiguities and increase uncertainties in the results. Fig. 1 shows a sketch of the method. First, training data is generated by running a physics model \(\mathcal{P}\) with various combinations of the observable \(\theta_{obs}\) and unobservable \(\theta_{un}\) parameter inputs. The outputs can be processed with a noise function \(\mathcal{N}\) to simulate the noise present in experimental measurements. Then, a deep learning model (neural network) \(\mathcal{D}\) is trained on this data to find particular observable features in the output of the physics model. The key to our method is to label the training data with _only_ the observable parameter values, so that the deep learning model finds the features without any knowledge of the unobservable parameters. By training over data generated with a wide range of unobservable parameter values, we ensure that the deep learning model works on data with any values of unobservable parameters. After training, the deep learning model is used to directly predict the values of observable parameters from measured data \(\mathcal{M}\). While a traditional analysis of the measured data with the physics model requires unobservable parameters to be specified, the deep learning model does not require the knowledge of the unobservables at the inference stage (Fig. 1**c**).

We demonstrate this method on the phenomenological _R_-matrix that is being critically used in nuclear physics [1; 3; 13; 14; 15; 16; 17; 18; 19]. The phenomenological _R_-matrix serves as a good example because it contains numerous and somewhat arbitrary unobservable parameters. It is based on separating the particle interaction space into internal and external regions. The complex many-body nuclear interactions are present in the internal region, whereas only the Coulomb interaction is considered in the external region. The boundary surface is defined by a sphere of radius known as the channel radius. The internal interactions are described by the _R_-matrix that is characterized by parameters of nuclear resonant states. From the _R_-matrix, expected observational data such as cross sections can be calculated. The parameters of nuclear states are determined by comparing the expected and experimental data. The extracted parameters provide the related nuclear properties, and the cross sections can be reproduced using the parameters for the purpose of extrapolation.

The main drawback of the phenomenological _R_-matrix arises from the dependency of fitted physical parame
ters on rather arbitrary choices of unobservable parameters [1; 3]. The _R_-matrix theory requires an infinite number of nuclear states in the calculation. However, the truncation of the number of states is necessary for the phenomenological approach, as generally only a limited number of nuclear states are known. The truncated states are compensated by the so-called background poles that are artificial states included to reproduce the effects of the truncated states. The strength of the background poles depends on the channel radius, which creates additional complications [1; 3; 13]. In practice, there are no strict physics rules to select valid choices of the channel radius and background poles. Different choices produce different fitting results because they create different theoretical predictions. This feature has caused difficulties in the phenomenological analysis and uncertainties of the results [1; 3].

We utilize the proposed method to extract observables without the requirement of unobservables used in the phenomenological \(R\)-matrix. One of the important applications of the phenomenological \(R\)-matrix is determining resonance parameters (e.g., spins \(J\), parities \(\pi\), energies \(E\), widths \(\Gamma\)) in nuclear reaction data. We built a deep learning model to find patterns of resonance parameters in reaction cross section spectra. Specifically, our goal was to extract the parameters of three resonance states from the well-known \({}^{12}\)C+\(p\) elastic scattering measurement of ref. [21] without using any information on the channel radius and background poles.

The training data for this task was prepared by the following specific process. We sampled possible sets of resonance parameters, along with random values of the channel radius (\(a_{c}\)) and random configurations of background poles (BGPs). We assumed a general situation where previous studies have partly determined the reso

Figure 1: Process of the current method to eliminate the unobservable parameters. The current method starts from (**a**) training data generation for (**b**) deep learning modeling to (**c**) inference on target measurement. The model requires only measurement data to extract physical objects at the inference stage.

\begin{table}
\begin{tabular}{c c c} \hline Parameters & Data Distribution & Model Performance \\ \hline \(J_{1}^{\pi}\) & \(\frac{1}{2}^{+}\) (Fixed) & _J_: -, \(\pi\): - \\ \(J_{2}^{\pi}\) & \(l\leq 2\) & _J_: 98.6 \%, \(\pi\): 100.0 \% \\ \(J_{3}^{\pi}\) & \(l=2\) & _J_: 99.0 \%, \(\pi\): - \\ \(E_{1}\) (MeV) & Normal(0.425, 0.003) & 0.00014 (0.03 \%) \\ \(E_{2}\) (MeV) & Normal(1.565, 0.010)\({}^{3}\) & 0.00118 (0.08 \%) \\ \(E_{3}\) (MeV) & Normal(1.610, 0.006)\({}^{3}\) & 0.00057 (0.04 \%) \\ \(\Gamma_{1}\) (keV) & Normal(30, 3) & 0.20 (0.67 \%) \\ \(\Gamma_{2}\) (keV) & Normal(50, 15) & 1.98 (4.44 \%) \\ \(\Gamma_{3}\) (keV) & Normal(60, 8) & 1.23 (2.08 \%) \\ \hline ANC (fm\({}^{-1/2}\))1  & Normal(1.81, 0.07) & - \\ \(a_{c}\) (fm) & 3 - 8 & - \\
5 BGPs [\(J^{\pi}\)] & \(l\leq 2\) & - \\
5 BGPs [\(E\)] & Uniform(3, 10) & - \\
5 BGPs [\(\Gamma\)] & Uniform(0, \(\Gamma_{W}\))2  & Uniform(0, \(\Gamma_{W}\))2  \\ \hline \end{tabular}
\end{table}
Table 1: Training data distribution and test results of model. Model performance is presented with accuracy and mean errors. Median values are taken to calculate the errors. The values in the parentheses are percentage errors. In total, 5 BGPs are included.
nance parameters. Table 1 shows the training data distributions of the parameters. Here, \(l\) is the relative orbital angular momentum of the particle pair. Possible spins and parities with the given \(l\) were uniformly sampled: e.g., \(J=0.5-2.5\) and \(\pi=+,-\) for \(l\leq 2\). We set the range of the channel radius to 3-8 fm that contains typical values being used in traditional \(R\)-matrix calculations. A level for each possible spin-parity was included as a background pole. The location and width of these levels were randomly determined within the given ranges.

We calculated cross section spectra corresponding to each set of parameters using \(R\)-matrix theory [16]. As the calculated spectra do not carry experimental noise, we added random noise to the spectra to reflect the measurement errors given by ref. [21]. We gave the log of the differential cross sections as the inputs to the model, while the labels are the corresponding (encoded) values of parameters \(J\), \(\pi\), \(\log(E)\), and \(\log(\Gamma)\) of three resonance states.

Transformers are currently one of the most effective deep learning architectures in a variety of applications [22; 23; 24; 25; 26; 27; 28]. We modify a Transformer to be suitable for this task. Details of the implemented architecture are described in the Methods section. Here, we summarize the primary specifications: First, as the input data are differential cross section spectra at three different angles, an extra dimension for the angles is needed in addition to the 2-dimensional cross section spectrum. Relations between the spectra of different angles are important in the analysis. However, Transformers consist of matrix multiplications that are operated only at the last two dimensions and are therefore limited within a spectrum at each angle. To accommodate this, we replaced the decoder input with the encoder input that is shifted in the angle dimension. By doing so, the matrix multiplications in the decoder are also operated between spectra at different angles. Second, because the different nature of our four parameters \(J\), \(\pi\), \(E\), and \(\Gamma\), we split the model output into four branches where two of them predict discontinuous values (classification) for \(J\) and \(\pi\), and the others predict continuous values (regression) for \(E\) and \(\Gamma\).

Figure 2: Errors of the model predictions for cross sections calculated with various unobservables. Res 1, 2, and 3 represent the first, second, and third resonances. The figures of the first background pole are given as representations.
Using deep ensembles, probabilistic deep learning is also used to quantify uncertainties of model predictions [29; 30; 31]. The idea of deep ensembles is straightforward and highly effective, featuring the use of multiple point-estimate models which are trained, using a negative log-likelihood loss function, from random initial model parameters and shuffled data samples [29; 32; 33]. The model prediction \(\mathbf{y}\) on given input data \(\mathbf{x}\) can be obtained from the predictions of these multiple models:

\[p(\mathbf{y}\,|\,\mathbf{x})=\frac{1}{N_{\text{model}}}\sum_{n=1}^{N_{\text{model}}}p( \mathbf{y}\,|\,\mathbf{x},\mathbf{\theta}_{n}), \tag{1}\]

where \(\mathbf{\theta}_{n}\) represents parameters of the \(n\)th model. The predictive probability for the classification of \(J\) or \(\pi\) is, therefore, the average of the probabilities predicted by the models. For the regression of \(\log(E)\) or \(\log(\Gamma)\), we assume that the likelihoods \(p(\mathbf{y}\,|\,\mathbf{x},\mathbf{\theta})\) follow Gaussian distributions. The predictive distribution \(p(\mathbf{y}|\mathbf{x})\) is further approximated as a Gaussian for simplicity, where its mean and variance are those of the mixture of the Gaussian likelihoods [31]. The final distribution is equivalent to the lognormal distribution of \(E\) or \(\Gamma\), which has several well-known advantages for representing positive-definite quantities [34]. More details of deep ensembles and the training process are presented in the Methods section.

The last column in Table 1 shows the average accuracies and errors of the trained model on the test dataset. The results show that the method and current model are quite effective. We also tested the model performance using the values of unobservable parameters used to calculate the cross sections during the training data generation. Fig. 2 shows errors of the predicted parameters for cross sections calculated with different values of eliminated parameters. The model predictions display little dependence - less than \(\sim\)2% and \(\sim\)1 keV for classification and regression, respectively - on these parameters. This demonstrates that the deep learning model performs well on data with any values of unobservable parameters. By using more training data or fine-tuning model hyperparameters, this performance can be further improved if needed.

Fig. 3 shows a calibration of the prediction uncertainty. For discrete parameters \(J\) and \(\pi\), the confidence level is the predicted probability from the softmax function, and the accuracy is the fraction of correct predictions. For continuous parameters \(E\) and \(\Gamma\), the confidence levels are calculated from their predicted ranges surrounding their median values, and the accuracies are the fraction of true data points included in the corresponding ranges. The dashed black line shows the ideal case where the confidence level is equal to the accuracy. It is clear that our probabilistic model achieves a performance close to

Figure 3: Reliability diagram of the current model. Each curve shows the estimated accuracy (fraction of true data points included in the confidence range) as a function of the corresponding confidence level. It is the standard method to evaluate the calibration of uncertainty.

Figure 4: Predictions of the \(E\) and \(\Gamma\) values from the measurement data. These lognormal distributions closely approximate Gaussian distributions. Previous results in the literature – M76 [21], A10 [16], D10 [1] – are also shown, where ref. [1] gave results using three different values of the channel radius.
this ideal case.

Fig. 4 shows the model predictions on the measurement data of ref. [21], along with the previous results of _R_-matrix analyses on the same data to extract \(E\) and \(\Gamma\)[16; 21; 1]. Our predictions of the spin-parities of second and third states are, at \(\sim\)100% probability, \(J_{2}=3/2\), \(\pi_{2}=-\), and \(J_{3}=5/2\), in agreement with the literature [35]. Our predictions of \(E\) and \(\Gamma\) are for the most part consistent with those of the previous studies. These consistent results on the well-known reaction data demonstrate the effectiveness of our model. While the previous works used different channel radius values and background pole configurations, the determinations of those parameters are not made in the current work. Our results are also tabulated in Table 2.

The current method is applicable to a wide variety of physics models. In this study, we focused on eliminating unobservable parameters, but this approach works for any subset of the parameters of a model. Interpolations or extrapolations of measured data are also possible by properly configuring the outputs of the deep learning model. Our approach of eliminating parameters could be used both when building generalized models that can address a variety of problems, as well as when building pre-trained models that need only to be slightly adjusted for a specific task [23; 36].

For cases where the values of the unobservable parameters in the physics model are unknown, the performance of the deep learning model - that is, its ability to accurately accurately predict observable parameters - can exceed that of the physics model. If the unobservables in the physics model could (by some means) be correctly specified, then the performance of the deep learning and physics models are comparable. The deep learning model can accurately identify distinctive features from the data, features that arise from variations of the observable parameters of the model. In cases where variations of unobservable physics model parameters can impact data features, the deep learning model predictions may have high uncertainties with diminished performance. By adding more training data samples or fine-tuning the deep learning model, these uncertainties can be reduced [37; 38; 39].

There are some additional benefits of our deep learning approach beyond the elimination of unobservable parameters. First, once a deep learning model is trained, the prediction calculation is normally much faster than typical Bayesian inference approaches that require a large number of Monte Carlo samplings [40; 41; 9]. Second, as shown in Fig. 3, the confidence level of the uncertainties in the deep learning model predictions can be explicitly tested through an analysis of a subset of data (a test dataset) [29; 32; 42; 43].

When extracting the values of observable model parameters from measured data, ambiguities can arise in some instances due to the presence of unobservable parameters. We demonstrate that a deep learning model can be trained without specifying such unobservable parameters by extracting values of observables directly from features in the dataset. By training the deep learning model without information on unobservables, these parameters are effectively eliminated. We demonstrate the effectiveness of this method on the phenomenological _R_-matrix model using a deep learning model based on Transformers. We trained this model to directly extract, from \({}^{12}\)C+\(p\) scattering measurement data, values of the observable parameters - resonance energies, widths, spins, and parities - with no need to specify values of the unobservable channel radius and background pole parameters used in traditional _R_-matrix analyses. We performed uncertainty quantification in deep learning using a probabilistic framework that provides calibrated uncertainties as tested on a subset of data samples. Our approach has broad applicability in physics to remove unobservables, or other subsets of model parameters, in analyzing datasets. Our approach has broad applicability to remove unobservables, or other subsets of model parameters, when analyzing data. We hope that this deep learning methodology can increase the discovery potential in a variety of fields.

## Acknowledgements

This work was supported by the National Research Foundation of Korea (NRF) grants funded by the Korea government (MSIT) (Grants No. 2016R1A5A1013277 and No. 2020R1A2C1005981). This work was also supported in part by the Institute for Basic Science (Grant No. IBS-R031-D1), and by the U.S. Dept. of Energy (DOE), Office of Science, Office of Nuclear Physics under contract DE-AC05-00OR22725. Computational works for this research were performed on the data analysis hub, Olaf in the IBS Research Solution Center. We thank D. Phillips for various discussions about this study.

## Methods

### Data Preparation

To generate the training data, we set ranges of the resonance parameters where the true parameter values are possibly located. These ranges can be obtained from previous studies. Alternatively, one can assume physically probable ranges. For example, one knows that a reaction with a high \(l\) value is largely suppressed by the centrifugal barrier, which allows an assumption of a \(l\) value less than a particular number. These assumptions prevail in the conventional analysis, where one does not try fitting with every possible \(l\) value but only some that are physically probable.

Each set of parameters was used to calculate the corresponding cross section using the _R_-matrix code, AZURE2 [16]. The target measurement data has differential cross section values on certain energy bin points at three different angles [21]. The calculated cross section values
on the same bin points at the three angles were taken as the input. As the locations of bin points and values of angles were fixed, they were null values in training and not included in the data. The input data samples were, therefore, 3-dimensional tensors with the shape of \((n_{\rm angle},n_{\rm bin},1)\). \(n_{\rm angle}\) is the number of measured angles, \(n_{\rm bin}\) is the number of bin points, and 1 represents the value of the cross section. A cross section value corresponds to a word in the natural language processing using the Transformer (see the next section) [22].

We added random noises in the calculated cross section values assuming Gaussian noises to reflect the measurement noises. The standard deviations of the Gaussians were obtained from the measurement error bars in ref. [21]. The addition of measurement noises is flexible so that any type of noises can be added. Non-Gaussian noises or noise on measured energy points can be handled by simply adding such noises to the training data.

### Model Architecture

The base architecture of the current model is from the Transformer [22]. Details of its description can be found in ref. [22]. Here, we briefly summarize its main features. The Transformer was originally constructed to handle sequence data for natural language processing. For translation purposes, the input and output are, for example, sentences from two different languages. The architecture is based on an encoder-decoder structure. Both encoder and decoder consist of multiple layers of multi-head attentions, fully connected feed-forward networks, residual connections, and layer normalizations [44; 45]. A multi-head attention splits the incoming data and applies multiple attention functions known as scaled dot-product attentions to let the model attend to the data points at multiple aspects [22]. An encoder layer contains a multi-head self-attention and feed-forward network. On the other hand, a decoder layer contains an additional multi-head attention known as an encoder-decoder attention performing with the encoder output.

The attentions find critical relations between two sequence data, known as a query and key, and use them to update another sequence data, known as a value. Queries, keys, and values in the self-attentions originate from the input. On the other hand, in the encoder-decoder attentions, queries originate from the decoder input, and keys and values originate from the encoder output. The obtained relations are often called attention (score) matrices, and they present locations of data points critical for the model predictions. Fig. 5 shows examples of the attention matrices from the current model in the prediction of the target measurement data.

The description of modifications to the original architecture is following. The inputs to the encoder and decoder are originally embedded and then passed to the positional encoding function to convert the sentences to vectors and encode the positional information of each word [22]. In this work, we removed such functions as the inputs are not exactly sentences but spectra with continuous values. We simply replaced the embedding and positional encoding with a linear layer that projects

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline Parameters & Current Model & Ref. [21]1  & Ref. [16] & Ref. [1]1  & Ref. [1]1  & Ref. [1]1  & Ref. [1]1  & Ref. [1]1  & Ref. [1]1  \\ \hline \(E_{1}\) (MeV) & 0.4247(2) & 0.424 & 0.426(3) & 0.427 & 0.427 & 0.427 \\ \(E_{2}\) (MeV) & 1.5592(12) & 1.558 & 1.556(1) & 1.560 & 1.559 & 1.558 \\ \(E_{3}\) (MeV) & 1.6023(7) & 1.604 & 1.602(2) & 1.603 & 1.604 & 1.606 \\ \(\Gamma_{1}\) (keV) & 34.2(3) & 33 & 34.1(8) & 33.8 & 32.9 & 30.9 \\ \(\Gamma_{2}\) (keV) & 55.2(9) & 55 & 57.9(17) & 51.4 & 51.4 & 51.3 \\ \(\Gamma_{3}\) (keV) & 50.2(8) & 50 & 48.3(19) & 48.1 & 48.1 & 47.8 \\ \(a_{c}\) (fm) & - & 4.0 & 3.4 & 4.0 & 5.0 & 6.0 \\ Background Poles & - & \(3/2^{+}\) state2  & 3 & 3 & 3 & 3 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Predictions of \(E\) and \(\Gamma\) for the measurement data from ref. [21] along with the previous results [1; 16; 21]. For the current model, the median values of the distributions are presented with the standard deviations as representations.
a cross section point to a vector. After the first linear layer, the shape \((n_{\text{angle}},n_{\text{bin}},1)\) was converted into \((n_{\text{angle}},n_{\text{bin}},d_{\text{model}})\), where \(d_{\text{model}}\) is a hyperparameter of the Transformer [22]. Additionally, we used the Gaussian Error Linear Units function as the activation function in the feed-forward networks of encoder and decoder layers [46].

The main operations in the Transformer are matrix multiplications that perform on the last two dimensions of data, whereas the current data has three dimensions. This feature limits the operations within a spectrum at each angle. The decoder in the Transformer performs the multi-head attention between its input and the encoder output. We gave the decoder input the same as the encoder input but shifted at the extra angle dimension, \([\mathbf{\theta_{1},\theta_{2},\theta_{3}}]\rightarrow[\mathbf{\theta_{3},\theta_{1}, \theta_{2}}]\). By doing so, the matrix multiplications in the attention were done between the data of different angles. The decoder input does not require the look-ahead mask originally demanded by the Transformer that uses target sequences as decoder input.

The hyperparameters in the model were:

* \(d_{\text{model}}=256\)
* \(d_{\text{df}}=256\)
* \(d_{\text{q}}=d_{\text{k}}=d_{\text{v}}=32\)
* \(h=8\)
* The number of encoder layers = 2
* The number of decoder layers = 2

We followed the notation of ref. [22].

The output of the modified Transformer was reshaped and passed to linear layers at the end of the model to return the four types of parameters.

### Probabilistic Model

We utilized a probabilistic framework to estimate the uncertainties of the model predictions. The method of deep ensembles is one of the most practical and effective uncertainty quantification approaches [29, 30, 31]. Deep ensembles consist of \(N_{\text{model}}\) point-estimate models trained with random initializations of model parameters and shuffling of data samples. In this study, we used 5 point-estimate models, \(N_{\text{model}}=5\). While it was first proposed as a non-Bayesian neural network, it has been shown that this method can be viewed as another Bayesian approach [29, 31, 47, 48]. Bayesian deep learning critically uses marginalization over possible sets of model parameters to obtain the predictive distribution [29, 48]:

\[p(\mathbf{y}\,|\,\mathbf{x},\mathbf{D})=\int_{\mathbf{\theta}}p(\mathbf{y}\,|\,\mathbf{x},\mathbf{\theta} ^{\prime})p(\mathbf{\theta}^{\prime}\,|\,\mathbf{D})d\mathbf{\theta}^{\prime}, \tag{2}\]

where \(\mathbf{x}\) and \(\mathbf{y}\) represents input and prediction, \(\mathbf{\theta}\) represents represents the model parameters, and \(\mathbf{D}\) represents training data. The prediction of deep ensembles in Equation 1 is equivalent to Equation 2 if the posterior is a mixture of delta functions. For this reason, deep ensembles can be seen as Bayesian marginalization with variational inference via a mixture of delta functions for the approximate distribution, \(p(\mathbf{\theta}\,|\,\mathbf{D})\sim q_{\mathbf{\phi}}(\mathbf{\theta})=(1/N_{\text{model}}) \sum_{p=1}^{N_{\text{model}}}\delta(\mathbf{\theta}=\mathbf{\theta}_{p})\). One of the advantages of deep ensembles is that they can capture different modes

Figure 5: Examples of the attention matrices when the target measurement data is passed as input. \(\mathbf{a}\), an attention matrix from a self-attention in the decoder. \(\mathbf{b}\), an attention matrix from an encoder-decoder attention. The top and left plots of each matrix show the corresponding input data. See the text for the explanation of attentions.
of posterior as multiple models are optimized to different optima [29, 47, 48]. More details on deep ensembles can be found in refs. [29, 31, 48].

We used negative log-likelihood (NLL) for the loss function in training. For regression, NLL for a data point \(i\) is

\[\text{NLL}= -\log p(y_{i}\left|f_{\text{NN}}(\mathbf{x}_{i};\mathbf{\theta})\right) \tag{3}\] \[= \frac{1}{2\sigma(\mathbf{x}_{i};\mathbf{\theta})^{2}}\|y_{i}-\mu(\mathbf{x}_ {i};\mathbf{\theta})\|^{2}\] (4) \[+\frac{1}{2}\text{log}(2\pi\sigma(\mathbf{x}_{i};\mathbf{\theta})^{2}),\]

where \(f_{\text{NN}}\) represents the neural network. The mean \(\mu\) and standard deviation \(\sigma\) are predicted by the neural network with the model parameter \(\mathbf{\theta}\) and input \(\mathbf{x}_{i}\). For classification, we also implemented the method proposed by ref. [38]. They presented a framework to combine aleatoric and epistemic uncertainties and showed the robustness of the method [37, 38, 39]. Logit vectors were sampled from a Gaussian distribution whose mean and standard deviation were predicted by the model. After passing them to a softmax function, the average of the softmax outputs was used to calculate NLL. For these reasons, the model in the current study had two outputs, \(\mu\) and \(\sigma\), for each parameter prediction.

We performed the training for 4,200 epochs with \(\sim\)21,000 training samples and \(\sim\)5,000 validation samples. The batch size was 768 on 3 GPU workers (NVIDIA Tesla V100 32GB SXM2). The number of samples is free to increase at the expense of training time. We used Adam optimizer with a weight decay of \(10^{-9}\)[49]. The learning rate was 0.0002 with a scheduler that decreased the learning rate by a factor of 0.99 at every 6 epochs. We constantly monitored the mean errors of predictions and NLL. After the training, we test the model with a relatively large number of samples, \(\sim\)30,000, to plot a reliability diagram with enough number of sample in each bin.

We plotted the reliability diagram for the performance test of the probabilistic model. For the classification tasks, the predicted probabilities were not evenly distributed to all probability ranges because of the high accuracies. For this reason, we set the bin locations to have an equal number of samples in each bin, not an equal bin interval. For regression tasks, we spanned the ranges from the median that corresponds to the mean of the original Gaussian \(N(\mu,\sigma^{2})\); for example, from \(e^{\mu-\sigma}\) to \(e^{\mu+\sigma}\) for the 68% confidence interval.
