Bayesian inference is an appealing approach for dealing with theoretical uncertainties and has been applied in different nuclear physics studies [1–16]. In the practice of Bayesian analyses, a sampling procedure is usually inevitable for approximating the posterior probability distribution of model parameters and for performing predictive computations. Various Markov chain Monte Carlo (MCMC) methods [17–21] are often used for this purpose, even for complicated models with high-dimensional parameter spaces. However, MCMC sampling typically requires many likelihood evaluations, which is often a costly operation in nuclear theory, and there is a need to explore other sampling techniques. In this paper, we review an established method called sampling/importance resampling (S/IR) [22–24] and demonstrate its use in realistic nuclear physics applications where we also perform comparisons with MCMC sampling.

In recent years, there has been an increasing demand for precision nuclear theory. This implies a challenge to not just achieve accurate theoretical predictions but also to quantify accompanying uncertainties. The use of ab initio many-body methods and nuclear interaction models based on chiral effective field theory (χEFT) has shown a potential to describe finite nuclei and nuclear matter based on extant experimental data (e.g. nucleon-nucleon scattering, few-body sector) with controlled approximations [25–29]. The interaction model is parametrized in terms of low-energy constants (LECs), the number of which is growing order-by-order according to the rules of a corresponding power counting [30–32]. Very importantly, the systematic expansion allows to quantify the truncation error and to incorporate this knowledge in the analysis [4–6, 10–14]. Indeed, Bayesian inference is an excellent framework to incorporate different sources of uncertainty and to propagate error bars to the model predictions. Starting from Bayes’ theorem p r θ | D ∝ L θ p r θ , (1) where p r ( θ | D ) is the posterior probability density function (PDF) for the vector θ of LECs (conditional on the data D ), L ( θ ) ≡ p r ( D | θ ) is the likelihood and pr( θ ) is the prior. Then for any model prediction one needs to evaluate the expectation value of a function of interest y( θ ) (target observables) according to the posterior. This involves integrals such as ∫ d θ y θ p r θ | D , (2) which can not be analytically solved for realistic cases. Fortunately, integrals such as Eq. 2 can be approximately evaluated using a finite set of samples { θ i } i = 1 N from p r ( θ | D ) . MCMC sampling methods are the main computational tool for providing such samples, even for high-dimensional parameter volumes [16]. However the use of MCMC in nuclear theory typically requires massive computations to record sufficiently many samples from the Markov chain. There are certainly situations where MCMC sampling is not ideal, or even becomes infeasible:

1) When the posterior is conditioned on some calibration data for which our model evaluations are very costly. Then we might only afford a limited number of full likelihood evaluations and our MCMC sampling becomes less likely to converge.

These are situations where one might want to use the S/IR method [23, 24], which can exploit the previous results of model evaluations to allow posterior probability updates at a much smaller computational cost compared to the full MCMC method. In the following sections we first review the S/IR method and then present both toy and realistic applications in which its performance is compared with full MCMC sampling. Finally, we illustrate limitations of the method by considering cases where S/IR fails and we highlight the importance of the so-called effective number of samples. More difficult scenarios, in which the method fails without a clear warning, are left for the concluding remarks.

The basic idea of S/IR is to utilize the inherent duality between samples and the density (probability distribution) from which they were generated [23]. This duality offers an opportunity to indirectly recreate a density (that might be hard to compute) from samples that are easy to obtain. Here we give a brief review of the method and illustrate with a toy problem.

Let us consider a target density h( θ ). In our applications this target will be the posterior PDF p r ( θ | D ) from Eq. 1. Instead of attempting to directly collect samples from h( θ ), as would be the goal in MCMC approaches, the S/IR method uses a detour. We first obtain samples from a simple (even analytic) density g( θ ). We then resample from this finite set using a resampling algorithm to approximately recreate samples from the target density h( θ ). There are (at least) two different resampling methods. In this paper we only focus on one of them called weighted bootstrap (more details of resampling methods can be found in Refs. [22, 23]).

Assuming we are interested in the target density h( θ ) = f( θ )/∫ f( θ ) d θ , the procedure of resampling via weighted bootstrap can be summarized as follows:

1) Generate the set { θ i } i = 1 n of samples from a sampling density g( θ ).

Intuitively, the distribution of θ * should be good approximation of h( θ ) when n is large enough. Here we justify this claim via the cumulative distribution function of θ* (for the one-dimensional case) p r θ * ≤ a = ∑ i = 1 n q i ⋅ H a − θ i = 1 n ∑ i = 1 n ω i ⋅ H a − θ i 1 n ∑ i = 1 n ω i n → ∞ → E g f θ g θ ⋅ H a − θ E g f θ g θ = ∫ − ∞ a f θ d θ ∫ − ∞ ∞ f θ d θ = ∫ − ∞ a h θ d θ , (3) with E g [ X ( θ ) ] = ∫ − ∞ ∞ X ( θ ) g ( θ ) d θ the expectation value of X(θ) with respect to g(θ), and H Heaviside step function such that H a − θ = 1 if θ ≤ a , 0 if θ > a . (4)

The above resampling method can be applied to generate samples from the posterior PDF h ( θ ) = p r ( θ | D ) in a Bayesian analysis. It remains to choose a sampling distribution, g( θ ), which in principle could be any continuous density distribution. However, recall that h( θ ) can be expressed in terms of an unnormalized distribution f( θ ), and using Bayes’ theorem 1) we can set f ( θ ) = L ( θ ) p r ( θ ) . Thus, choosing the prior pr( θ ) as the sampling distribution g( θ ) we find that the importance weights are expressed in terms of the likelihood, q i = L ( θ i ) / ∑ j = 1 n L ( θ j ) . Assuming that it is simple to collect samples from the prior, the costly operation will be the evaluation of L ( θ i ) . Here we make the side remark that an effective and computationally cost-saving approximation can be made if we manage to perform a pre-screening to identify (and ignore) samples that will give a very small importance weight. We also note that the above choice of g( θ ) = pr( θ ) is purely for simplicity and one can perform importance resampling with any g( θ ).

In Figure 1 we follow the above procedure and give a simple example of S/IR to illustrate how to get samples from a posterior distribution. We consider a two-dimensional parametric model with θ = (θ ₁, θ ₂). Given data D obtained under the model we have: p r θ 1 , θ 2 | D = L θ 1 , θ 2 p r θ 1 , θ 2 ∬ L θ 1 , θ 2 p r θ 1 , θ 2 d θ 1 d θ 2 . (5) For simplicity and illustration, the joint prior distribution for θ ₁, θ ₂ is set to be uniform over the unit square as shown in Figure 1A. In this example we also assume that the data D follows a multivariate Student-t distribution such that the likelihood function is L θ 1 , θ 2 = Γ ν + p / 2 Γ ν / 2 ν p / 2 π p / 2 | Σ | 1 / 2 1 + 1 ν θ − μ T Σ − 1 θ − μ − μ + p / 2 , (6) where the dimension p = 2, the degrees of freedom ν = 2, the mean vector μ = (0.2, 0.5) and the scale matrix Σ = [[0.02, 0.005], [0.005, 0.02]].

The importance weights q _i are then computed for n = 2000 samples drawn from the prior (these prior samples are shown in Figure 1A). The resulting histogram of importance weights is shown in Figure 1B. Here the weights have been rescaled as q ̃ i = q i / max ( { q } ) such that the sample with the largest probability mass corresponds to 1 in the histogram. We also define the effective number of samples, n _eff, as the sum of rescaled importance weights, n eff = ∑ i = 1 n q ̃ i . Finally, in Figure 1C we show N = 20, 000 new samples { θ i * } i = 1 N that are drawn from the prior samples { θ i } i = 1 n according to the probability mass q _i for each θ _i . The blue and green contour lines represent (68% and 90%) credible regions for the resampled distribution and for the Student-t distribution, respectively. This result demonstrates that the samples generated by the S/IR method give a very good approximation of the target posterior distribution.

Now that we have reviewed the basic idea of the S/IR method, we move on to present realistic applications of the resampling technique in nuclear structure calculations. Here we study Bayesian inference involving the ΔNNLO chiral interaction [33] with explicit inclusion of delta isobar degree of freedom at next-to-next-to-leading order. In Weinberg’s power counting the ΔNNLO interaction model is parametrized by 17 LECs, with four pion-nucleon LECs (c _1,2,3,4) that are inferred from pion-nucleon scattering data and 13 additional LECs that should be inferred from extant experimental data of low-energy nucleon-nucleon scattering and bound-state nuclear observables.

For this application we treat only a subset of the parameters as active and keep the other LECs fixed at values taken from the ΔNNLO_GO(450) interaction [34]. Specifically, we consider deuteron observables and use seven active model parameters: c _1,2,3,4, C ̃ 3S1 , C _3S1 , C _E1. Our Gaussian likelihood contains three data wih independent errors: the deuteron ground state energy E, its point-proton radius R _p and one-body quadrupole moment Q with experimental and theoretical targets from Refs. [35–37]. Note that the target point-proton radii were transformed from experimental charge radii using the same relation as in Ref. [38]. For the target Q we use the theoretical result obtained by the CD-Bonn [37] model. With these simplified conditions, we perform S/IR as well as MCMC sampling to study 1) the posterior PDF for the LECs and 2) posterior predictive distributions (PPDs) for selected few-body observables. This application therefore allows a straightforward comparison of the two different sampling methods in a realistic setting. We note that the inclusion of all 17 LECs as active parameters would have required more careful tuning of the MCMC sampling algorithm and corresponding convergence studies.

It is the computation of observables, e.g., for likelihood evaluations, which is usually the major, time-consuming bottleneck in Bayesian analyses using MCMC methods. In this application, the statistical analysis is enabled by the use of emulators which mimic the outputs of many-body solvers but are faster by orders of magnitude. The emulators employed here for the ground-state observables of the deuteron, and later for few-body observables, are based on eigenvector continuation [39–41]. These emulators allow to reduce the computation time from seconds to milliseconds while keeping the relative error (compared with full no-core shell model calculation) within 0.001%. Unfortunately, emulators are not yet available for all nuclear observables. The MCMC sampling of posterior PDFs, or the evaluation of expectation integrals such as Eq. 2, will typically not work for models with observables that require heavy calculations.

The experimental target values and error assignments for the calibration observables used to condition the posterior PDF are listed in the upper half of Table 1. In this study we assume a normally-distributed likelihood, and consider different sources of error when calibrating the model predictions with experimental data. The errors are assumed to be independent. They include experimental, ɛ _exp, model (χEFT truncation) discrepancy, ɛ _model, many-body method, ɛ _method, and emulator, ɛ _em, errors. The χEFT truncation errors are estimated based on order-by-order calculations as in Ref. [33]. More details on the determination of the error scales can be found in Ref. [42].

Furthermore, we take advantage of previous studies and incorporate information about c _1,2,3,4 from a Roy-Steiner analysis of pion-nucleon scattering data [43] and identify a non-implausible domain for C ̃ 3S1 , C _3S1 , C _E1 from a history matching approach in Ref. [42]. With this prior knowledge we set up the prior distribution of the seven LECs as the product of a multivariate Gaussian for c _1,2,3,4 and a uniform distribution for C ̃ 3S1 , C _3S1 , C _E1 ¹ . We note that the use of history matching is very beneficial for both S/IR and MCMC sampling. For S/IR it allows to select a sampling distribution that promises a large overlap with the target distribution and it identifies prior samples that are likely to have large weights in the resampling step. For MCMC, the non-implausible samples from history matching serve as good starting points for the walkers and thereby give faster convergence.

In this paper we reviewed an established sampling method known as S/IR. Specifically, we applied importance resampling using the weighted bootstrap algorithm and sampled the posterior PDF for selected LECs of the ΔNNLO interaction model conditioned on deuteron observables. The resulting PDF and PPD were compared with those obtained from MCMC sampling and a very good agreement was found. We also demonstrated Bayesian updating using S/IR by the addition of ³H and ⁴He observables to the calibration data set. As expected, the predictions of ³H and ⁴He observables were improved, but also the description of the deuteron ground-state energy which could be explained by the correlation structure between E(²H) and R _p (⁴He). Finally, we illustrated some limitations of the S/IR method that were signaled by small numbers of effective samples. We found that such situations occured when the likelihood became too complex for the limited model, or when prior samples failed to resolve a very peaked posterior that resulted from small tolerances. We also argued that prior knowledge of the posterior landscape is very useful to avoid possible failure scenarios that might not be signaled by the number of effective samples.