Direct comparison between Bayesian and frequentist uncertainty quantification for nuclear reactions

TL;DR

This paper systematically compares Bayesian and frequentist methods for uncertainty quantification in nuclear reactions, revealing that Bayesian approaches provide more accurate and flexible uncertainty estimates than traditional frequentist methods.

Contribution

It provides the first direct comparison between Bayesian and frequentist uncertainty quantification methods in nuclear reaction analysis, highlighting their differences and advantages.

Findings

Bayesian approach offers more flexible parameter exploration.

Frequentist uncertainties are narrower than Bayesian ones.

Bayesian uncertainties align better with empirical data.

Abstract

Until recently, uncertainty quantification in low energy nuclear theory was typically performed using frequentist approaches. However in the last few years, the field has shifted toward Bayesian statistics for evaluating confidence intervals. Although there are statistical arguments to prefer the Bayesian approach, no direct comparison is available. In this work, we compare, directly and systematically, the frequentist and Bayesian approaches to quantifying uncertainties in direct nuclear reactions. Starting from identical initial assumptions, we determine confidence intervals associated with the elastic and the transfer process for both methods, which are evaluated against data via a comparison of the empirical coverage probabilities. Expectedly, the frequentist approach is not as flexible as the Bayesian approach in exploring parameter space and often ends up in a different minimum. We also show that the two methods produce significantly different correlations. In the end, the frequentist approach produces significantly narrower uncertainties on the considered observables than the Bayesian. Our study demonstrates that the uncertainties on the reaction observables considered here within the Bayesian approach represent reality more accurately than the much narrower uncertainties obtained using the standard frequentist approach.

Tables

Table 1. Table 1: List of reactions studied: beam energy (col. 2), data reference (col. 3), percentage of the data that falls into the 95% confidence interval obtained in the frequentist approach (col. 4) and the Bayesian approach (col. 5).

Reaction	Energy (MeV)	Ref.	%Data(F)	%Data(B)
48Ca(p,p)	14.03	Lombardi et al. (1972)	70	100
48Ca(n,n)	12	Mueller et al. (2011)	61	100
48Ca(p,p)	25	McCamis et al. (1986)	51	94
90Zr(p,p)	12.7	Dickens et al. (1968)	66	100
90Zr(n,n)	10	Wang and Rapaport (1990)	92	100
90Zr(p,p)	22.5	Ball et al. (1964)	67	92
208Pb(p,p)	16	Makofske et al. (1972)	62	92
208Pb(n,n)	16.9	Floyd (1981)	62	76
208Pb(p,p)	35	van Oers et al. (1974)	44	88

Full text

uncertainty quantification, nucleon elastic scattering, transfer nuclear reactions, optical potential fitting

Introduction:

Over the past few years, Bayesian methods have rapidly drawn much attention in the field of low energy nuclear physics. They have opened new doors for theoretical predictions: providing a means to rigorously quantify uncertainties and the potential to help plan for future experiments at rare isotope facilities worldwide. Bayes statistics has been used in studies of the nuclear force Furnstahl et al. (2015); Melendez et al. (2017), of nuclear stability Neufcourt et al. (2018), and in nuclear astrophysics simulations (e.g. Orford et al. (2018); Lim and Holt (2018)). Our group has also been exploring Bayesian methods in the context of nuclear reactions Lovell and Nunes (2018), particularly in connection to assessing the uncertainties in the predicted cross section coming from the optical potentials. While the explosion of applications of Bayesian statistics to low energy nuclear theory is very exciting, it also calls for special scrutiny.

Up to a decade ago, there were very few studies done on uncertainty quantification in low energy nuclear theory, and those efforts relied primarily on $\chi^{2}$-minimization techniques stepping through parameter space in the maximal direction of the local gradient. Relevant to our subfield, the major global optical potentials being used currently in reaction calculations Becchetti and Greenlees (1969); Varner et al. (1991); Koning and Delaroche (2003) have been optimized in this way, which from now on we shall refer to as the frequentist approach. One important conclusion from this approach was the presence of strong correlations between some parameters of the optical potential, as discussed in Lovell et al. (2017); King et al. (2018).

In contrast, the Bayesian framework often relies on Markov Chain Monte Carlo (MCMC) methods to sample parameter space and obtain posterior predictions from the product of the likelihood function and prior distribution, without any normality assumption on the distribution of parameters. Thus, one should not assume that the minimum obtained from the frequentist approach matches the results obtained in the Bayes-MCMC method, and, moreover, that the uncertainties estimated by both methods are consistent. It is also possible that the insight obtained in the Bayesian analysis does not match the knowledge that we have established over the years using the frequentist methods. 111As a disclaimer, what we here call the frequentist or Bayesian approaches corresponds to a subset of a larger array of methods available under these umbrellas.

For all these reasons, it is important and timely to directly compare the two approaches. Despite the many recent applications mentioned above, to our knowledge, there exists no systematic and controlled comparison between these two methodologies in our field. This is the goal of the present study.

Methods and numerical details:

In this study, we focus on deuteron induced reactions on heavy ions, within a three-body model Thompson and Nunes (2009). The inputs to the cross sections calculated are nucleon-nucleus optical potentials $U_{opt}$. To capture the complex many-body effects of nucleon-nucleus scattering, the potential contains both a real part and an imaginary part to account for flux that leaves the elastic channel. Typically, optical potentials contain i) a real Woods-Saxon volume term with parameters $V$, $r$, $a$ for the depth, radius, and diffuseness ii) an imaginary Woods-Saxon volume term with parameters $W$, $r_{w}$, $a_{w}$, iii) a surface (derivative of a Woods-Saxon) imaginary term with parameters $W_{s}$, $r_{s}$, and $a_{s}$, iv) a spin-orbit term, and v) a Coulomb term for charged projectiles Thompson and Nunes (2009). In this work, we fitted the real volume, imaginary volume, and imaginary surface parameters and kept the spin-orbit and Coulomb terms fixed. We then used elastic scattering angular distribution data to constrain these free parameters.

We first consider neutron and proton scattering on 48Ca, 90Zr and 208Pb. Table 1 (columns 1-3) contains a summary of all the data used to constrain the necessary optical potentials. We chose only data sets with a wide angular distribution and small error bars, and assign a 10% uncertainty on all data to account for experimental error and model inadequacy; this also prevents overfitting. Because we eventually propagate the uncertainties to the transfer (d,p) channel, the energies for the neutron and proton elastic scattering are carefully chosen to be close to either half the deuteron beam energy or the energy of the proton in the exit channel. Starting from the original Becchetti and Greenlees (BG) global parameterization for the optical potential Becchetti and Greenlees (1969), we use both, the frequentist and Bayesian method, to obtain posterior distributions for the parameters and 95% confidence intervals for elastic scattering data. In the case of 90Zr, the resulting geometry parameters (radius and diffuseness) for proton and neutron optical potentials resulting from the frequentist approach were significantly different from each other, when starting from the BG parameterization. Since this is physically implausible, we use the parameters from the fitting of 90Zr(n,n) to initialize the 90Zr(p,p) reaction, only.

The frequentist approach follows the methods described in detail in Lovell et al. (2017); King et al. (2018), considering only the uncorrelated $\chi^{2}$ function. To initialize the fitting procedure, the original BG parameters were used and allowed to vary within a wide window (physical limits for each parameter). Then $1600$ parameter sets were pulled from the multivariate normal distribution induced by the $\chi^{2}$ around the optimal parameters (Lovell et al., 2017, Eq. (6)); 95% confidence intervals were constructed by removing the highest and lowest 2.5% of the predicted values for the cross sections at each calculated angle.

We use the Bayesian MCMC implementation Lovell and Nunes (2018), with the same numerical details. We choose wide Gaussian priors, centered at the original BG value, with a standard deviation equal to the mean value of the distribution. This choice has proven to ensure that the parameter space is appropriately sampled Lovell and Nunes (2018). Again a total of 1600 parameter sets were collected. In contrast to the frequentist method, 95% confidence intervals are constructed by taking the densest 95% of the 1600 cross section values calculated at each angle. The wrapper codes developed for both the frequentist and the Bayesian analyses make use of the reaction codes fresco and sfresco Thompson (1988).

Results:

We illustrate our results in detail for reactions on the 48Ca target. We start by analyzing nucleon elastic scattering on 48Ca in Fig.1. In the first row, we present 95% confidence intervals for neutron and proton elastic-scattering angular distributions as a function of scattering angle. For both methods, the uncertainty intervals appear to provide a good description of the data. However, the Bayesian approach results in wider confidence intervals than the frequentist. This is clear in Figs.1(panel d, e and f), where the percentage width of each uncertainty band is computed as a function of scattering angle. Here, $\epsilon$ was obtained by taking the width of the confidence interval at a given angle in the distribution and dividing by the best fit (average) cross section value in the frequentist (Bayesian) approach.

The reliability of the confidence intervals obtained from the two statistical methods can be assessed by the empirical coverage of the confidence intervals. We show in Figs.1 (g)(h) and (i) the empirical coverage probability curves of each model, namely the percentage of the experimental data which actually falls within the predicted confidence intervals with nominal value spanning $[0,100\%]$. Because this corresponds to training data, one might expect that reality matches predictions; this ideal situation corresponds to the black diagonal line in Figs.1(h, i, j). Points below this line correspond to underestimation of the confidence intervals, whereas points above the diagonal correspond to confidence intervals that are too large (a lesser harm). For each of the three reactions shown in Fig.1, the Bayesian model provides an accurate quantification of its uncertainty for the larger confidence intervals. In contrast, the frequentist approach undershoots when large confidence intervals are considered. This indicates that the 95% confidence intervals obtained with the frequentist approach are unrealistically narrow.

We now inspect the posterior distributions for the parameters and the correlations between parameters. In Fig.2, we show, along the diagonal, the posterior distributions predicted by both approaches (blue for the frequentist and orange for the Bayesian): 48Ca(n,n) at 12 MeV in Fig.2(a); 48Ca(p,p) at 14 MeV in Fig.2(b); and 48Ca(p,p) at 25 MeV in Fig.2(c). Note that these correlation plots require larger statistics: for Fig.2 we needed to collect $100,000$ parameter sets. The peak of the posterior distributions obtained for the parameters associated with the real part of the optical potential are nearly identical in both approaches, however the Bayesian procedure produce densities with large and sometimes asymmetric tails which thus differ significantly from Gaussian distributions. The differences in the peak of the posteriors for the parameters associated with the imaginary components are more noticeable. Note that some of these imaginary parameters could not be included in the fit for the frequentist approach, as the minimization procedure drove the parameters into unphysical values. For those cases, the parameters in the frequentist model were fixed at some intermediate physically plausible values, and only the Bayesian posteriors are shown.

Fig.2 also contains the scatter plots associated with each pair of parameters. Elongated elliptical forms demonstrate high correlations between parameter pairs, while circular scatter plots indicate low correlations. It is clear from Fig.2 that the frequentist approach finds strong correlations between parameters. This was also the conclusion in Lovell et al. (2017) and has been well established in the field of nuclear reactions Thompson and Nunes (2009). However, the Bayesian approach provides a very different picture: we find systematically that only the depth and radius of the real part are strongly correlated, in all cases here considered. All other scatter plots have approximately circular shape indicating that the parameters are weakly correlated.

The reduction of the number of fitted parameters in the frequentist approach could be responsible for the additional correlations seen in the scatter plots in the frequentist approach. To test this hypothesis, we have rerun the Bayesian procedure, fixing the same parameters as in the frequentist approach. The results lead to the same overall conclusion. Given that the posterior distributions for the Bayesian deviate strongly from Gaussian, one can thus expect that the bivariate distributions of the Bayesian model (corresponding to the scatter plots in Fig.2) would also deviate. Our hypothesis is thus that this is causing a bias in the frequentist approach, which leads to erroneous conclusions, both in the magnitude of the uncertainties and in the correlations between parameters.

Finally, we now look at the results when propagated through a three-body reactions model to determine the one-nucleon transfer (d,p) cross sections as done in Lovell and Nunes (2018); King et al. (2018). The adiabatic wave approximation (ADWA) Johnson and Tandy (1974) was the reaction model used to describe the transfer reaction. Apart from the optical potentials, all other inputs used in the calculations were the same as those used in Lovell and Nunes (2018). Using the parameter posteriors for each of the relevant nucleon-target optical potentials, we obtain 95% confidence intervals for the transfer reaction (our wrapper codes call the reaction code NLAT Titus et al. (2016) for this purpose). We show, in Fig.3(a), the 95% confidence intervals for the angular distributions of 48Ca(d,p)49Ca(g.s.) at 19.3 MeV. Again the blue-slanted-hashed area corresponds to the results using the frequentist approach and the orange-vertical-hashed area corresponds to those using the Bayesian approach. The transfer cross sections have been normalized to the data taken from Ref.Metz et al. (1975) at the peak of the angular distribution. Fig.3(b) shows that the percent uncertainty of the confidence intervals in the Bayesian approach is larger than the frequentist across all angles. Finally, we also show the percentage of this test data that falls within the given confidence interval in Fig.3(c). As previously with the training data, the Bayesian approach is more reliable, particularly when considering higher-levels of confidence.

In order to be able to draw general conclusions we expand our study to include reactions on 90Zr and 208Pb, as listed in Table 1. We include in columns 4 and 5, the percentage of data that actually falls in the 95% confidence interval predicted by the frequentist approach and the Bayesian approach, respectively. For all but one of the cases considered, the frequentist approach does not provide a reliable estimate for the confidence interval. Indeed, the confidence intervals predicted by the frequentist approach are too narrow and the fraction of the data that actually falls in them is much smaller than the nominal 95%. On the contrary, the Bayesian approach tends to be more conservative for the nucleon elastic scattering on 48Ca and 90Zr at the lower energies and makes predictions that are very close to the correct values for most of the other cases.

Conclusions:

In this work we compare directly two methods to evaluate uncertainties in reaction theory, namely the standard frequentist approach and the Bayesian framework. We perform a systematic study of nucleon scattering on three different target nuclei and use elastic scattering data to constrain the optical potential parameters. Our study is also controlled in that we make the same assumptions in both the frequentist and the Bayesian approaches. We conclude that the 95% confidence intervals generated by the frequentist approach are unrealistically narrow, as opposed to the Bayesian approach that provide a correct assessment of the uncertainty. In addition, the frequentist approach generates strong correlations between parameters that are not seen in the Bayesian approach. We attribute these differences to a strong deviation from the $\chi^{2}$-based Gaussian distribution. When propagating these uncertainty to (d,p) reactions, we arrive at the same conclusions deduced from the elastic scattering analysis.

There is, a priori, no reason for the stark differences found between these two uncertainty quantification approaches to be specific to reaction theory, therefore it is desirable that other subfields (e.g. atomic and molecular physics), perform similar studies.