A Bayesian Approach for Parameter Estimation and Prediction using a Computationally Intensive Model

Bayesian methods have been very successful in quantifying uncertainty in physics-based problems in parameter estimation and prediction. In these cases, physical measurements y are modeled as the best fit of a physics-based model $\eta(\theta)$ where $\theta$ denotes the uncertain, best input setting. Hence the statistical model is of the form $y = \eta(\theta) + \epsilon$, where $\epsilon$ accounts for measurement, and possibly other error sources. When non-linearity is present in $\eta(\cdot)$, the resulting posterior distribution for the unknown parameters in the Bayesian formulation is typically complex and non-standard, requiring computationally demanding computational approaches such as Markov chain Monte Carlo (MCMC) to produce multivariate draws from the posterior. While quite generally applicable, MCMC requires thousands, or even millions of evaluations of the physics model $\eta(\cdot)$. This is problematic if the model takes hours or days to evaluate. To overcome this computational bottleneck, we present an approach adapted from Bayesian model calibration. This approach combines output from an ensemble of computational model runs with physical measurements, within a statistical formulation, to carry out inference. A key component of this approach is a statistical response surface, or emulator, estimated from the ensemble of model runs. We demonstrate this approach with a case study in estimating parameters for a density functional theory (DFT) model, using experimental mass/binding energy measurements from a collection of atomic nuclei. We also demonstrate how this approach produces uncertainties in predictions for recent mass measurements obtained at Argonne National Laboratory (ANL).

💡 Research Summary

The paper addresses the challenge of performing Bayesian parameter estimation and prediction when the underlying physics‑based model is computationally expensive to evaluate. In many scientific problems the data are modeled as y = η(θ) + ε, where η(θ) is a deterministic simulation that may take minutes or hours for a single set of input parameters θ, and ε represents measurement error. When η(·) is nonlinear and the dimension of θ is moderate to high, the posterior distribution π(θ | y) becomes analytically intractable and must be explored with Markov chain Monte Carlo (MCMC). Conventional MCMC requires thousands to millions of calls to η(·), which is infeasible for costly codes such as nuclear density functional theory (DFT).

To overcome this bottleneck the authors adopt a Bayesian model‑calibration framework that replaces the expensive simulator with a statistical surrogate (emulator). They first generate an ensemble of m = 183 DFT runs over a space‑filling Latin hyper‑cube design in the 12‑dimensional parameter space that defines the UNEDF1 Skyrme functional. For each of the 92 nuclei considered (28 spherical, 47 deformed, and 17 newly measured neutron‑rich isotopes) the DFT code outputs a predicted mass. The collection of simulator outputs forms a training set for a Gaussian process (GP) model. The GP provides a probabilistic mapping from the input vector t to the simulated mass η(t), delivering both a mean prediction and a predictive variance for any new t without additional DFT evaluations.

The statistical model for the observed masses incorporates three components: (i) the GP emulator for η(θ), (ii) a discrepancy term δ that captures systematic model‑data mismatch, and (iii) a measurement‑error term ε. The likelihood is constructed by integrating the GP predictive distribution with the variance of δ and ε. Prior distributions for the 12 physical parameters are taken as broad uniform or normal ranges reflecting expert knowledge, while hyper‑priors are placed on the GP covariance parameters and the discrepancy variance.

Because the GP surrogate is analytically tractable, MCMC can now be performed on the posterior π(θ,δ,σ²_GP,σ²_δ | y) using only the GP mean and covariance, eliminating the need for costly DFT calls. The authors employ a Metropolis‑Hastings algorithm, generating thousands of posterior samples in a matter of hours on a standard workstation.

Results show substantial contraction of the posterior for the DFT parameters relative to the original UNEDF1 prior ranges, and the inferred discrepancy term is small but non‑zero, indicating that the DFT model does not perfectly reproduce the experimental masses. Using the calibrated parameters and the GP emulator, the authors predict masses for the 17 newly measured neutron‑rich nuclei from the CARIBU facility at Argonne National Laboratory. Predictive distributions (means and 95 % credible intervals) are compared with the actual measurements; the calibrated model reduces the root‑mean‑square error by roughly 30 % compared with the uncalibrated UNEDF1 predictions, and the credible intervals appropriately capture the experimental uncertainties.

The paper highlights several methodological contributions: (1) a practical workflow for building a GP emulator from a modest number of expensive simulations, (2) a full Bayesian treatment that jointly estimates physical parameters, emulator hyper‑parameters, and model discrepancy, and (3) a demonstration that calibrated predictions with quantified uncertainty can be obtained for data not used in the calibration step. Limitations discussed include the sensitivity of the GP to the design of the training points, challenges in estimating high‑dimensional covariance structures, and the need for careful specification of the discrepancy model.

In conclusion, the authors provide a scalable Bayesian calibration approach that makes it feasible to perform rigorous uncertainty quantification for computationally intensive models. While illustrated on nuclear DFT, the methodology is broadly applicable to any domain where expensive simulations must be combined with experimental data—such as climate modeling, aerospace design, or high‑energy physics—offering a path to reliable predictions without prohibitive computational cost.