BAYESIAN CALIBRATION WITH ADAPTIVE
MODEL DISCREPANCY

Nicolas Leoni, Olivier Le Maitre, Maria-Giovanna Rodio, Pietro Marco Congedo

Congedo

► To cite this version:

Nicolas Leoni, Olivier Le Maitre, Maria-Giovanna Rodio, Pietro Marco Congedo. BAYESIAN CALIBRATION WITH ADAPTIVE MODEL DISCREPANCY. International Journal for Uncertainty Quantification, 2024, 14 (1), pp.19-41. 10.1615/Int.J.UncertaintyQuantification.2023046331. hal-04323574v2

HAL Id: hal-04323574

https://inria.hal.science/hal-04323574v2

Submitted on 5 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# BAYESIAN CALIBRATION WITH ADAPTIVE MODEL DISCREPANCY

Nicolas Leoni,¹,²,* Olivier Le Maître,³ Maria-Giovanna Rodio,² & Pietro Marco Congedo¹

¹Inria, Centre de Mathématiques Appliquées, Ecole polytechnique, IPP, Route de Saclay, 91128 Palaiseau Cedex, France

²Commissariat à l'Energie Atomique et aux Energies Alternatives, ISAS, D36, 91191 Gif-sur-Yvette Cedex, France

³CNRS, Inria, Centre de Mathématiques Appliquées, Ecole polytechnique, IPP, Route de Saclay, 91128 Palaiseau Cedex, France

*Address all correspondence to: Nicolas Leoni, Inria, Centre de Mathématiques Appliquées, Ecole polytechnique, IPP, Route de Saclay, 91128 Palaiseau Cedex, France, E-mail: nicolas.leoni@inria.fr

Original Manuscript Submitted: mm/dd/yyyy; Final Draft Received: mm/dd/yyyy

We investigate a computer model calibration technique inspired by the well-known Bayesian framework of Kennedy and O'Hagan. We tackle the full Bayesian formulation where model parameter and model discrepancy hyperparameters are estimated jointly and reduce the problem dimensionality by introducing a functional relationship that we call the Full Maximum a Posteriori (FMP) method. This method also eliminates the need for a true value of model parameters that caused identifiability issues in the KOH formulation. When the joint posterior is approximated as a mixture of Gaussians, the FMP calibration is proved to avoid some pitfalls of the KOH calibration, namely missing some probability regions and underestimating the posterior variance. We then illustrate two numerical examples where both model error and measurement uncertainty are estimated together. Using the solution to the full Bayesian problem as a reference, we show that the FMP results are accurate, robust and avoid the need for high-dimensional Markov Chains for sampling.

**KEY WORDS:** Uncertainty Quantification, Bayesian calibration, Model error, Model discrepancy, Identifiability
1

# 1. INTRODUCTION

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

Models are widely employed to understand, interpret, and predict phenomena. As the phenomenon of interest is increasingly complex, the models must be more elaborated and, eventually, have to be solved numerically using computers. The software used to solve a scientific model using numerical techniques is called a computer model or a computer code. Computer codes have become central to scientific work, thanks to the increasing computing power available and their cost, usually lower than full-scale experiments. The accuracy and reliability of computer code predictions are consequently of paramount importance today [1].

Models usually involve parameters that are not perfectly known. Calibration procedures aim to estimate these parameters using experimental observations. Bayesian techniques are well suited to the calibration task, as they offer a natural framework for plausible reasoning [2] and enable the incorporation of expert knowledge. This knowledge can take many forms, like favoring particular values of the parameters using an informative prior distribution or performing posterior predictive checks to ensure the statistical assumptions led to physically acceptable predictions.

All models are imperfect representations of reality; they involve *model error*, or *model discrepancy*, affecting the quality of predictions. The pioneering article of Kennedy and O'Hagan (KOH) [3] was the first to introduce model discrepancy in the calibration process. Their Bayesian framework represents the model discrepancy with Gaussian processes. It has been widely applied in various fields such as aerodynamics [4–7], fluid mechanics [8,9] and solid mechanics ([10,11]). Compared to standard calibration, inferring the model discrepancy jointly with the parameters yields more objective and less informative parameter posteriors. The work of Kennedy and O'Hagan [3] has triggered a whole body of calibration literature dealing with computer models with high-dimensional outputs [12–14] and with computer model validation [15,16]. The nature of model discrepancy and its treatment in the calibration process was also extensively discussed in [17–19]. For instance, different calibration methods with model discrepancy were recently proposed in [20], [21] and [22].

Although the calibration with model discrepancy has encountered tremendous successes, multiple questions remain unanswered. One question concerns the choice of representation of the model error. A Gaussian process is added to the model output in KOH to represent the discrepancy. Other works use constant bias, random walks [23], or deterministic functions [24]. For Gaussian process models, the user must provide the prior mean, the covariance function structure, and the prior of its hyperparameter. Setting these quantities is not an obvious task and may significantly impact the calibration results. Another approach to calibration in the presence of model error embeds representations to add uncertainties in the model parameters [25,26]. The frequentist point of view is taken in [27,28] permitting to cure a potential *L*₂-*inconsistent calibration of the original KOH method*.
In twenty years of applications of the Bayesian calibration framework, two classes of methods have emerged: fully Bayesian methods and modular methods. The fully Bayesian methods aim at solving the complete Bayesian formulation of the problem accounting for all uncertainties related to the data, the surrogate for the computer model if necessary, and the estimation of the model discrepancy. Although theoretically superior to others, the resolution of the full Bayesian problems is seldom practical due to their complexity and computational cost [12], and numerical approximation techniques are necessary to solve these problems in practice [13]. The modular methods tackle the complexity of the calibration by introducing simplified formulations of the fully Bayesian problem. Typically, they rely on sequential pointwise estimations of hyperparameters. These strategies significantly reduce the sample spaces' dimensionality, improving the whole computational efficiency of the calibration. The authors of [29], who coined the name, also claimed a theoretical advantage of modular approaches, as the separated sequential estimations avoid confounding issues.

The KOH framework is a modular approach based on a pointwise estimation of the model discrepancy before estimating the model parameters. As a result, the method underestimates the uncertainty related to the model discrepancy. Although generally considered a "second-order uncertainty," recent studies [19] show that underestimating uncertainties in the model discrepancy might impact the posterior predictions, which are of utmost importance to assess the calibration quality and predict new values with objective confidence intervals. Therefore, there is a need for new modular techniques that approximate the fully Bayesian calibration problem better while retaining the advantage of reduced sample spaces. In this paper, we propose a new calibration method, named Full Maximum a Posteriori (FMP), which involves a parametric estimation of the model discrepancy hyperparameters. Specifically, the model error is explicitly dependent on model parameter values in our calibration framework. The resulting predictive uncertainty split into two contributions: the calibrated model error and the residual uncertainty. The FMP method leads to simple approximations of the full Bayesian parameters and predictive posterior distributions, for reduced cost compared to the full Bayesian inference. However, the FMP method corrects the defects of the KOH estimations when the posterior of the hyper-parameters has several modes.

The organization of the paper is as follows. In section 2 we present the Bayesian calibration problem, introducing the general framework with model discrepancy in section 2.1. Sections 2.2 and 2.3 discuss the fully Bayesian and modular methods, respectively. In section 3, we introduce the FMP method, discussing its underlying hypotheses and providing the expressions of the resulting parameters' posterior and posterior predictive density. We show that the FMP estimation of the parameter has higher plausibility than any other modular approach. In section 3.3, we contrast the calibration results for the KOH and FMP for analytic examples consisting of a Gaussian or a mixture of Gaussians
1 for the full Bayesian posterior. The superiority of the FMP method, compared to the KOH approach, is demonstrated
2 in these analytical examples. In section 4, we consider more complex examples. The first one uses a simple model
3 having two sets of parameter values yielding a comparable explanation of the observations. The second example
4 concerns the calibration of a boiling model on experimental measurements. The FMP method proves superior to the
5 KOH method in approaching the full Bayesian solution in the two examples.

## 2. BAYESIAN CALIBRATION OF COMPUTER CODES

### 2.1 General framework and reference solution

8 The first requirement for parameter calibration is to specify a statistical model to explain the observations of a quantity
9 y. We restrict ourselves to a scalar observed quantity, $y \in \mathbb{R}$, of which we have $n$ observations collected in the
10 observations vector $\mathbf{y}_{\text{obs}} \in \mathbb{R}^n$. The observations consist of measurements for different experimental conditions; a
11 particular experimental condition is described by a vector $\mathbf{x} \in \mathbb{R}^d$. For instance, in an experiment with fluids flowing
12 through a pipe, one may measure the pressure loss along the pipe for different conditions $\mathbf{x}$, consisting of the flow
13 rate, the fluid density, temperature,.... We denote $\{\mathbf{x}_i\}_{1\le i\le n}$ the set of experimental conditions corresponding to the
14 observations in $\mathbf{y}_{\text{obs}}$.

15 Let $f$ be a computer model that provides predictions of $y$. The inputs for $f$ includes the conditions $\mathbf{x}$ and other
16 model parameters $\theta \in \Theta \subset \mathbb{R}^p$. The latter are usually imperfectly known, and the goal of the calibration procedure
17 is to learn their value from the observation $\mathbf{y}_{\text{obs}}$. In the example, $\theta$ might consists of the coefficients of an empirical
18 friction model.

19 In practice, the observations involve a measurement error $\epsilon(\mathbf{x})$. The simplest observation model relates the model
20 predictions to the observations through the additive form

$$ \mathbf{y}_{\text{obs}}(\mathbf{x}) = f(\mathbf{x}, \theta) + \epsilon(\mathbf{x}). \quad (1) $$

21 A statistical model for the measurement error is needed. For simplicity, we shall assume that the measurement error
22 does not depend on the condition $\mathbf{x}$ and that the errors of different observations are independent. For a Gaussian
23 distribution of the error we have $\epsilon(\mathbf{x}_i) \sim N(0, \sigma_e^2)$, where $\sigma_e^2$ is the mean-squared measurement error. These settings
24 are standard and used in many calibration studies, such as [30].

25 As discussed below, a Bayesian procedure can be employed to learn the parameters $\theta$ with the observation model
26 in (1). The procedure updates the prior distribution of $\theta$ to obtain its posterior distribution and, subsequently, posterior
predictions using the computer code. The Bayesian procedure was successfully applied to numerous calibration problems. However, there are situations where it yields inconsistent results, with residuals (differences between posterior predictions and observations) that do not follow the assumed measurement error distribution. Such inconsistencies may have two causes: an incorrect measurement error distribution and an observation model that is too poor. This work considers situations with adequate measurement error distribution and focuses on the observations model. For instance, the residuals may be highly correlated with significant mean values when the measurement errors are known to be unbiased and independent for different $x_i$. These inconsistencies are interpreted in the work of Kennedy and O'Hagan [3] as a model error describing the discrepancy between the observed data and the model predictions from even the best-fitting parameter values. The authors proposed to account for the model error by completing the observation model with a model discrepancy term $z(\mathbf{x})$, leading to the following calibration equation:

$$y_{\text{obs}}(\mathbf{x}) = f(\mathbf{x}, \theta) + z(\mathbf{x}) + \epsilon(\mathbf{x}). \qquad (2)$$

Note that the original formulation in [3] includes a multiplicative coefficient $\rho$ on the computer model prediction $f$, but the expression in (2) is generally adopted in the calibration literature.

More statistical assumptions are needed to proceed with the calibration of the parameters. First, measurement error and model discrepancy are usually supposed independent: $z(\mathbf{x}_i) \perp \epsilon(\mathbf{x}_i)$. Second, the model discrepancy is taken as a Gaussian process: $z|\psi_z \sim GP(\mu(\cdot), c(\cdot, \cdot))$, where $\mu$ and $c$ are the prior mean and covariance of the process. The covariance of $z$ depends on some hyperparameters $\psi_z$. For the prior mean, we set $\mu = 0$ as advocated in [12,31] for better identifiability. However, $\mu(\mathbf{x})$ can be expressed as a linear combination of basis functions in the context of Universal Kriging [32].

We denote $\psi$ the vector of hyperparameters of the entire error model, so that $\psi = (\psi_z, \sigma_\epsilon)$ or $\psi = \sigma_\epsilon$ depending if model error is considered or not in the calibration. The space of hyperparameters is $\Psi \subset \mathbb{R}^h$.

In a Bayesian framework, the posterior distribution $p(\theta, \psi|y_{\text{obs}})$ provides all information about parameters and hyperparameters after seeing the observations. Its derivation requires the specification of a prior distribution $p(\theta, \psi)$ reflecting our knowledge about their values before seeing the observations and of a likelihood function $p(y_{\text{obs}}|\theta, \psi)$. Then, the posterior follows the Bayes' Rule:

$$p(\theta, \psi|y_{\text{obs}}) \propto p(\theta, \psi)p(y_{\text{obs}}|\theta, \psi). \qquad (3)$$
1 The log-likelihood function corresponding to the calibration equation without model error (1) is the probability density function of a multivariate normal law with diagonal covariance matrix:

$$ \log p(\mathbf{y}_{\text{obs}}|\theta, \sigma_{\epsilon}) = -\frac{n}{2}\log 2\pi - \frac{n}{2}\log \sigma_{\epsilon}^2 - \frac{1}{2\sigma_{\epsilon}^2}\|\mathbf{y}_{\text{obs}} - \mathbf{f}_{\theta}\|^2, \quad (4) $$

3 where $\mathbf{f}_{\theta} = (f(\theta, \mathbf{x}_1), f(\theta, \mathbf{x}_2), \dots, f(\theta, \mathbf{x}_n))^T$ is the vector of evaluations of the computer code at the observation points.

5 The log-likelihood function corresponding to the calibration equation with model discrepancy (2) writes:

$$ \log p(\mathbf{y}_{\text{obs}}|\theta, \mathbf{\psi}_z, \sigma_{\epsilon}) = -\frac{n}{2}\log 2\pi - \frac{1}{2}\log \det(\Sigma_{\psi} + \sigma_{\epsilon}^2 \mathbf{I}_n) - \frac{1}{2}(\mathbf{y}_{\text{obs}} - \mathbf{f}_{\theta})^T (\Sigma_{\psi} + \sigma_{\epsilon}^2 \mathbf{I}_n)^{-1} (\mathbf{y}_{\text{obs}} - \mathbf{f}_{\theta}), \quad (5) $$

6 where $\mathbf{I}_n$ is the identity matrix of size $n$ and $\Sigma_{\psi}$ is the prior covariance matrix of $z$: $(\Sigma_{\psi})_{i,j} = c_{\psi_z}(\mathbf{x}_i, \mathbf{x}_j)$.

7 The choice of prior distributions in Bayesian statistics is not straightforward and affects the calibration results. 8 Model parameters and hyperparameters are generally considered independent a priori: $p(\theta, \psi) = p(\theta)p(\psi)$. 9 parameters might be physical quantities with a range of plausible values in the literature. Consequently, it is often 10 possible to formulate informative priors for them. For the hyperparameters of the model discrepancy, classical choices 11 are uniform priors and Jeffrey's prior (e.g., $p(\sigma) \propto 1/\sigma$). Inverse-gamma priors for correlation lengths $l$ have the 12 advantage of ruling-out the cases of vanishing $l$, inducing a clear separation in the estimation of model error and 13 measurement error (assumed uncorrelated). 14

14 The marginal posterior distribution of $\theta$ is obtained by integration over the hyperparameters $\psi$:

$$ p_{\text{Bayes}}(\theta) := p(\theta|\mathbf{y}_{\text{obs}}) = \int_{\Psi} p(\theta, \psi|\mathbf{y}_{\text{obs}}) \, d\psi. \quad (6) $$

15 The posterior predictive density is used to make predictions of the true process at any experimental condition $\mathbf{x}_*$, 16 possibly not observed:

$$ p(y(\mathbf{x}_*)|\mathbf{y}_{\text{obs}}) = \int_{\Theta} \int_{\Psi} p(y(\mathbf{x}_*)|\theta, \psi, \mathbf{y}_{\text{obs}})p(\theta, \psi|\mathbf{y}_{\text{obs}}) \, d\theta \, d\psi, \quad (7) $$

17 The distribution $p(y(\mathbf{x}_*)|\theta, \psi, \mathbf{y}_{\text{obs}})$ is normal, obtained from the Gaussian Process predictive equations, which
1 writes [33]:

$$y(\mathbf{x}_*)|\Psi, \Theta, \mathbf{y}_{\text{obs}} \sim \mathcal{N}(\mu_{\text{pred}}, \sigma_{\text{pred}}^2),$$
$$\mu_{\text{pred}} = f(\mathbf{x}_*, \Theta) + \mathbf{k}_*^T (\mathbf{\Sigma}_{\Psi} + \sigma_e^2 \mathbf{I}_n)^{-1} (\mathbf{y}_{\text{obs}} - \mathbf{f}_{\Theta}), \quad (8)$$
$$\sigma_{\text{pred}}^2 = c_{\Psi}(\mathbf{x}_*, \mathbf{x}_*) - \mathbf{k}_*^T (\mathbf{\Sigma}_{\Psi} + \sigma_e^2 \mathbf{I}_n)^{-1} \mathbf{k}_*.$$

2 with $\mathbf{k}_* = (c_{\Psi}(\mathbf{x}_*, \mathbf{x}_1), c_{\Psi}(\mathbf{x}_*, \mathbf{x}_2), \dots, c_{\Psi}(\mathbf{x}_*, \mathbf{x}_n))^T$.

3 In the following, we refer to the general calibration framework with model error presented above as the *full Bayesian* framework. The distributions defined in equations (6) and (7) constitute the exact solution to the calibration problem, as was proposed for example in [12,13]. The interest in the full Bayesian calibration comes from its rigorous derivation and the characterization of the posterior uncertainties in the model parameters and model discrepancy. 4 5 6 7 Further, the framework is flexible and can adapt to various situations. For instance, in the case of an expensive 8 computer code requiring its substitution by a surrogate model, the framework can be extended introducing a surrogate 9 error model with new hyperparameters.

10 11 12 13 14 15 16 17 18 Despite its apparent simplicity, this calibration can not be solved precisely in most applications. The reason is that it requires an accurate estimation of the posterior density, which might exhibit substantial variations and many local maxima over the joint space $\Theta \times \Psi$. Further, for calibrations involving multiple quantities and experimental configurations, the number of hyperparameters might become significant, making the estimation in high dimensions very complex. Monte-Carlo techniques [34] and other advanced sampling techniques have been employed to generate samples from the high dimensional posterior distributions and approximate the required integrals. The following section discusses some works proposed in the literature on direct full Bayesian calibration. However, the so-called modular approaches ([29]) have gained popularity as alternatives to full Bayesian calibration. These approaches, discussed in section 2.3, reduce the dimensionality of the sample space by estimating point values for some hyperparameters.

## 2.2 Fully Bayesian approaches

20 21 22 23 24 An example of application of full Bayesian method is in [12], where the full joint density is sampled using the Metropolis-Hastings algorithm in dimension six (one parameter and five hyperparameters), allowing a complete representation of the posterior uncertainty. In [13], the authors deal with a high-dimensional output code and a sample space with a dimension greater than 30. Besides, more than 18,000 observations are considered making one evaluation of the likelihood function costly. Consequently, the authors relied on reduced basis representations, block matrix
1

inversion formulas, and adapted informative priors to make the problem tractable and ease the posterior sampling.

2

3

4

5

6

7

8

9

10

11

Another example of a fully Bayesian technique is presented in [15], where the authors chose a covariance function for *z* that involves variance and correlation length hyperparameters. Before the sampling, correlation lengths are estimated with **Maximum Likelihood (ML)**, as well as any eventual prior mean for *z*. ML estimates of variance hyperparameters serve to construct a prior for them. The MCMC sampling itself is simplified by allowing only variance hyperparameters to vary (as well as other hyperparameters involved in the surrogate for the computer model), using Gibbs updates as a conjugate prior was chosen. The authors declare that combining ML estimates with sampling achieves an efficient approximation to the Full Bayesian calibration. Note that they claim that ML estimates are a sufficient approximation to Full Bayesian, but only concerning hyperparameters for the computer model surrogate; they recognize that the variations of model discrepancy hyperparameters must be considered. Conjugate priors are also central to the sampling with Gibbs steps proposed in [17].

12

13

14

15

Yuan and Ng [35] approximated the full Bayesian solution using an Expectation-Maximization algorithm, which can be interpreted as a repetition of a modular method of the first class (see next section). In [36], the full posterior distribution was sampled with a Gibbs sampler. In [37], the authors proposed a weighted normal approximation to the posterior to reduce the overall cost of evaluation and accommodate large numbers of hyperparameters.

16

17

18

19

20

21

22

23

24

25

26

We now discuss two calibration frameworks that account for the uncertainty in model error through its depen-
dence on the model parameters. In [20], a "true value" of model parameters is defined as the minimizer of a loss
function. It leads to an orthogonality condition for the model discrepancy and the gradient of the computer model.
For each $\Theta$, a model discrepancy is estimated using a covariance function that verifies the orthogonality condition.
This calibration aims to find the "true parameters value", a questionable objective when the model predictions can
not satisfactorily represent the observations. Nonetheless, using an adaptive representation of model error instead of
a distribution is a way to include multiple model discrepancies that might be very different in structure. In [38], the
calibration was recast as an optimization problem, where a single value of $\Theta$ is sought to maximize the likelihood
of the corresponding model discrepancy term. These applications show that the full Bayesian calibration remains of
practical interest despite being less popular than modular approaches due to its cost and the difficulties involved in
sampling.

27

**2.3 Modular approaches**

28

29

The principle underlying modular approaches is to separate into groups (modules) the treatment of parameters and hyperparameters during the calibration and obtain a problem of sequential estimations. In [29], three modules are
considered in the calibration of computer models: the surrogate of $f$, the measurement error in the observations, and the model discrepancy. Each module is then separately estimated. This approach breaks the high-dimensional calibration problem into smaller ones. One expects that these separated problems are more accessible to sample than the original problem. The modular approach can also improve the identifiability of the hyperparameters belonging to different modules.

Since the calibration is sequential, the calibration problem for the modules must be adapted, and, in particular, assumptions regarding subsequent modules not yet calibrated are needed. Since the posterior distribution of the parameters $\theta$ is arguably the goal of the calibration, the last module estimates $p(\theta|\mathbf{y}_{\text{obs}}, \mathbf{\psi} = \hat{\mathbf{\psi}})$ where the estimator $\hat{\mathbf{\psi}}$ of $\mathbf{\psi}$ have been estimated in previous modules such that $p(\theta|\mathbf{y}_{\text{obs}}, \mathbf{\psi} = \hat{\mathbf{\psi}}) \approx p_{\text{Bayes}}(\theta)$. We now distinguish between two classes of methods proposed in the calibration literature.

In the first class of methods, one starts by determining a value $\theta_0$ such that $f(\mathbf{x}, \theta_0)$ is an accurate prediction of $y(\mathbf{x})$. Then, the hyperparameters are estimated, usually by ML, considering $\theta = \theta_0$. With our notations, this amounts to the following estimator:

$$ \hat{\mathbf{\psi}}_1 = \arg\max_{\mathbf{\Psi}} p(\mathbf{y}_{\text{obs}}|\mathbf{\psi}, \theta_0). \qquad (9) $$

The methods in the first class differ by the definition of $\theta_0$: some works advocate for elicited values such as the prior mean [29,39]. Other common approaches estimate $\theta_0$ from the observation, either minimizing the residuals between the model prediction and the observations [40], or solving the calibration problem without model error (1) to propose a plausible value [11,19,41,42]. These approaches inherently assess the relevance of introducing a model discrepancy term by judging the quality of the model predictions using $\theta = \theta_0$. In [43], it is demonstrated that $f(\mathbf{x}, \theta_0) + z(\mathbf{x})$ can be an accurate prediction of $y(\mathbf{x})$, even for a poor value of $\theta_0$, because the model discrepancy term $z$ can compensate for a wide range of residuals.

In the second class of methods, one averages the posterior over the model parameters value to obtain the distribution $p(\mathbf{y}_{\text{obs}}|\mathbf{\psi})$, and then the hyperparameters are estimated using:

$$ \hat{\mathbf{\psi}}_2 = \arg\max_{\mathbf{\Psi}} p(\mathbf{y}_{\text{obs}}|\mathbf{\psi}). \qquad (10) $$

This approach is used for instance in [14,31,44]. It was originally proposed by Kennedy and O'Hagan [3,45], with also the consideration of the hyperparameter prior, leading to the Maximum A Posteriori (MAP) estimate:

$$ \hat{\mathbf{\psi}}_{\text{KOH}} = \arg\max_{\mathbf{\Psi} \in \Psi} p(\mathbf{\psi}|\mathbf{y}_{\text{obs}}) = \arg\max_{\mathbf{\Psi} \in \Psi} p(\mathbf{\psi}) \int_{\Theta} p(\theta)p(\mathbf{y}_{\text{obs}}|\theta, \mathbf{\psi}) d\theta, \qquad (11) $$
1 and the following estimation for the parameters posterior distribution:

$$p_{\text{KOH}}(\theta) \propto p(\theta)p(y_{\text{obs}}|\theta, \hat{\psi}_{\text{KOH}}). \qquad (12)$$

2 Note also that an approximation of $p(y_{\text{obs}}|\psi)$ is sometimes used to reduce the computational cost; Kennedy and O'Hagan propose to use a Gaussian approximation in [45]. Compared to the first class of methods, the computation of $p(y_{\text{obs}}|\psi)$ is generally more costly than estimating a single parameter value $\theta_0$. Compared to the first class, the second class of methods benefits from the complete account of the parameter uncertainty when estimating the model discrepancy.

7 The two classes of modular methods are compared in [19], where the second class estimates higher model discrepancies, showing the benefit of including the parameter uncertainty. In general, the modular approaches provide fairly accurate parameter posteriors. Nonetheless, due to the single point estimate for $\psi$, they do not account for the complete uncertainty about the model discrepancy term. In what follows, we propose a calibration technique where the uncertainty about the values of $\psi$ is estimated to include all relevant uncertainties in the calibration framework.

## 3. THE FULL MAXIMUM A POSTERIORI METHOD

13 This section introduces the FMP approximation. It starts with its founding assumptions and the resulting approximations to the distributions of interest, presented in section (3.1). The FMP method is shown to produce marginal likelihoods greater than any other similar estimation technique in section (3.2). Then, we derive the exact FMP, KOH and Bayes calibration solutions, in the case of a Gaussian joint posterior in section (3.3.1), and a mixture of Gaussians with well-separated modes in section (3.3.2).

### 3.1 Calibration framework

#### 3.1.1 Fundamental hypotheses

20 We follow the idea that the plausibility of parameter values should be compared by considering their most favorable model discrepancy distributions. This is done by first recognizing that the model discrepancy term depends on the model parameter:

$$z_{\theta}(\mathbf{x}) = y(\mathbf{x}) - f(\mathbf{x}, \theta). \qquad (13)$$
1 Therefore, the calibration equation becomes:

$$y_{\text{obs}}(\mathbf{x}) = f(\mathbf{x}, \theta) + z_{\theta}(\mathbf{x}) + \epsilon(\mathbf{x}). \quad (14)$$

2 This equation appeared previously in the work of [20], but we do not modify the covariance of $z_{\theta}$ as they do. Instead,

3 we perform a joint estimation of $\theta$ and $\psi$ like in the full Bayesian framework of [12]. To make the estimation cheaper

4 we introduce a functional relationship between the parameters and the "optimal hyperparameters" defined as follows

$$\psi_{\text{FMP}}(\theta) = \arg\max_{\psi} p(\psi)p(y_{\text{obs}}|\theta, \psi), \quad (15)$$

5 where FMP stands for "Full Maximum A Posteriori", a combination between full Bayesian analysis and MAP estimation. Two additional assumptions are introduced to derive the approximation of the posterior.

7 1. For all $\theta \in \Theta$, the distribution $p(\psi|\theta, y_{\text{obs}})$ can be approximated as a point mass distribution at its mode $\hat{\psi}_{\text{FMP}}(\theta)$. This assumption is standard in calibration contexts. For instance, the underlying assumption of KOH's framework is a point mass approximation of the distribution $p(\psi|y_{\text{obs}})$ at its mode. Here we rely on a point mass approximation of the distributions $p(\psi|\theta, y_{\text{obs}})$. In this way, we do not consider one value of $\psi$ that fits the error for all model predictions (when $\theta$ varies). Instead, the errors for each model prediction can be fitted well by individual values of $\psi$.

13 2. We assume that, for all $\theta \in \Theta$ with significant posterior,

$$\int_{\Theta} p(\theta'|y_{\text{obs}})p(y_{\text{obs}}|\theta, \hat{\psi}_{\text{FMP}}(\theta')) d\theta' \approx p(y_{\text{obs}}|\theta, \hat{\psi}_{\text{FMP}}(\theta)). \quad (16)$$

14 This assumption is exploited below to simplify the expression and computation of the FMP posterior. In practice, it holds true for multi-modal posteriors that are not well-separated in both $\theta$ and $\psi$. Otherwise, the approximation in 16 introduces a bias in the relative strengths of the modes, but retain the global multi-modal structure of the $\theta$'s posterior, in contrast with the KOH method. This is illustrated in the example of section 3.3.2.
### 1.3.1.2 Posterior approximation and predictive density

With the previous assumptions, the parameter posterior, named $p_{\text{Bayes}}(\theta|\mathbf{y}_{\text{obs}})$, can be expressed as

$$
\begin{align*}
p_{\text{Bayes}}(\theta|\mathbf{y}_{\text{obs}}) &= \int_{\Psi} p(\theta|\mathbf{\psi}, \mathbf{y}_{\text{obs}})p(\mathbf{\psi}|\mathbf{y}_{\text{obs}}) d\mathbf{\psi} \\
&= \int_{\Psi} p(\theta|\mathbf{\psi}, \mathbf{y}_{\text{obs}}) \left( \int_{\Theta} p(\mathbf{\psi}|\theta', \mathbf{y}_{\text{obs}})p(\theta'|\mathbf{y}_{\text{obs}}) d\theta' \right) d\mathbf{\psi} \\
&= \int_{\Theta} p(\theta'|\mathbf{y}_{\text{obs}}) \left( \int_{\Psi} p(\theta|\mathbf{\psi}, \mathbf{y}_{\text{obs}})p(\mathbf{\psi}|\theta', \mathbf{y}_{\text{obs}}) d\mathbf{\psi} \right) d\theta' \\
&\approx \int_{\Theta} p(\theta'|\mathbf{y}_{\text{obs}})p(\theta|\mathbf{y}_{\text{obs}}, \hat{\mathbf{\psi}}_{\text{FMP}}(\theta')) d\theta' \tag*{(Assumption 1, Sec. 3.1.1)} \\
&\approx \int_{\Theta} p(\theta'|\mathbf{y}_{\text{obs}})p(\theta)p(\mathbf{y}_{\text{obs}}|\theta, \hat{\mathbf{\psi}}_{\text{FMP}}(\theta')) d\theta' \tag*{(Bayes Theorem)} \\
&\approx p(\theta)p(\mathbf{y}_{\text{obs}}|\theta, \hat{\mathbf{\psi}}_{\text{FMP}}(\theta)) \doteq p_{\text{FMP}}(\theta|\mathbf{y}_{\text{obs}}) \tag*{(Assumption 2, Sec. 3.1.1).} \quad (17)
\end{align*}
$$

2 We denote with $p_{\text{FMP}}$ the last term of Eq. (17), i.e. the probability density proportional to $p(\theta)p(\mathbf{y}_{\text{obs}}|\theta, \hat{\mathbf{\psi}}_{\text{FMP}}(\theta))$.

3 From the previous equation, $p_{\text{FMP}}$ is a reasonable approximation of $p_{\text{Bayes}}$.

4 Note that for most covariance functions, $c_{\psi}(x, x')$ is a continuous function of $\psi$ for all $(x, x')$. By the Maximum 5 Theorem [46,47] it can be shown that $p_{\text{FMP}}$ is a continuous function of $\theta$. With the additional assumption that for 6 each value of $\theta$ there is a unique value of optimal hyperparameters (which is not true in general), we have that $\hat{\psi}_{\text{FMP}}$ 7 is a continuous function of $\theta$. These continuity properties pave the way for calibration techniques based on surrogate 8 modelling of the FMP quantities that will be explored in subsequent work.

The predictive density for the true process at an unobserved experimental condition $x^*$ is given by:

$$
\begin{align}
p(y(\mathbf{x}^*)|\mathbf{y}_{\text{obs}}) &= \int_{\Psi} \int_{\Theta} p(y(\mathbf{x}^*)|\theta, \mathbf{\psi}, \mathbf{y}_{\text{obs}})p(\mathbf{\psi}|\theta, \mathbf{y}_{\text{obs}})p(\theta|\mathbf{y}_{\text{obs}}) d\theta d\mathbf{\psi} \nonumber \\
&= \int_{\Theta} p(\theta|\mathbf{y}_{\text{obs}}) \left( \int_{\Psi} p(y(\mathbf{x}^*)|\theta, \mathbf{\psi}, \mathbf{y}_{\text{obs}})p(\mathbf{\psi}|\theta, \mathbf{y}_{\text{obs}}) d\mathbf{\psi} \right) d\theta \nonumber \\
&\approx \int_{\Theta} p(\theta|\mathbf{y}_{\text{obs}})p(y(\mathbf{x}^*)|\theta, \mathbf{\psi}) = \hat{\psi}_{\text{FMP}}(\theta), \mathbf{y}_{\text{obs}}) d\theta. \tag{18}
\end{align}
$$

9 The mean and variance of the predictive density can be computed explicitly. It comes

$$
\begin{align}
\mathbb{E}[y(\mathbf{x}_*)|\mathbf{y}_{\text{obs}}] &= \mathbb{E}_{\theta}[\mathbb{E}[y(\mathbf{x}_*)|\mathbf{y}_{\text{obs}}, \theta]] \nonumber \\
&= \underbrace{\mathbb{E}_{\theta}[f(\mathbf{x}_*, \theta)]}_{\text{averaged model prediction at } \mathbf{x}_*} + \underbrace{\mathbb{E}_{\theta}[\mathbf{k}_*^T (\mathbf{\Sigma} + \sigma_{\epsilon}^2 \mathbf{I}_n)^{-1} (\mathbf{y}_{\text{obs}} - \mathbf{f}_{\theta})],}_{\text{averaged model discrepancy at } \mathbf{x}_*} \tag{19}
\end{align}
$$
1 and

$$
\begin{align}
\mathrm{Var}[y(\mathbf{x}_*)|\mathbf{y}_{\mathrm{obs}}] &= \mathrm{Var}_{\theta}[\mathbb{E}[y(\mathbf{x}_*)|\mathbf{y}_{\mathrm{obs}}, \boldsymbol{\theta}]] + \mathbb{E}_{\theta}[\mathrm{Var}[y(\mathbf{x}_*)|\mathbf{y}_{\mathrm{obs}}, \boldsymbol{\theta}]] \nonumber \\
&= \mathrm{Var}_{\theta}[f(\mathbf{x}_*, \boldsymbol{\theta}) + \mathbf{k}_*^T (\boldsymbol{\Sigma}_{\hat{\psi}_{\mathrm{FMP}}(\boldsymbol{\theta})} + \sigma_e^2 \mathbf{I}_n)^{-1} (\mathbf{y}_{\mathrm{obs}} - \mathbf{f}_{\theta})] \nonumber \\
&\qquad \underbrace{\text{uncertainty in the corrected model}} \tag{20} \\
&\quad + \mathbb{E}_{\theta}[c_{\hat{\psi}_{\mathrm{FMP}}(\boldsymbol{\theta})}(\mathbf{x}_*, \mathbf{x}_*) - \mathbf{k}_*^T (\boldsymbol{\Sigma}_{\hat{\psi}_{\mathrm{FMP}}(\boldsymbol{\theta})} + \sigma_e^2 \mathbf{I}_n)^{-1} \mathbf{k}_*] \nonumber \\
&\qquad \underbrace{\text{residual uncertainty}}
\end{align}
$$

2 where $\mathbf{k}_* = (c_{\hat{\psi}_{\mathrm{FMP}}(\boldsymbol{\theta})}(\mathbf{x}_*, \mathbf{x}_1), c_{\hat{\psi}_{\mathrm{FMP}}(\boldsymbol{\theta})}(\mathbf{x}_*, \mathbf{x}_2), \dots, c_{\hat{\psi}_{\mathrm{FMP}}(\boldsymbol{\theta})}(\mathbf{x}_*, \mathbf{x}_n))^T$.

Equation (20) provides a decomposition of the predictive error into two contributions. The first contribution, named "uncertainty in the corrected model", is the variability of the corrected predictions. It is significant when the corrected predictions vary with $\boldsymbol{\theta}$, indicating that the model discrepancy term cannot adequately correct the predictions. The second term, named "residual uncertainty," is the part of the posterior uncertainty that the corrected model cannot explain. It corresponds to the remaining uncertainty in calibration using a finite amount of noisy observations.

The FMP framework avoids ambiguity in defining the "true value" of parameters. Nevertheless, other types of identifiability issues might still plague the calibration process. If two distinct values of $\boldsymbol{\theta}$ with the same prior value lead to the same computer model prediction, one might deem them unidentifiable. Likewise, if the prior distribution of $z$ is not suited to the actual model discrepancy, optimal hyperparameters might be meaningless and posterior samples of $z$ would be inaccurate. These issues can be addressed by testing several prior distributions for $\boldsymbol{\theta}$, $z$, and $\psi$. The multivariate Gaussian example in section 3.3.2 below gives a graphical representation of these identifiability cases.

## 3.2 Marginal likelihood of the FMP approximation

We now examine the relevance of the FMP method in the scope of Bayesian Model Comparison (BMC). For a model $M$ and an observations set $\mathbf{y}_{\mathrm{obs}}$, the marginal likelihood $p(\mathbf{y}_{\mathrm{obs}}|M)$ is the probability of observing $\mathbf{y}_{\mathrm{obs}}$ given that its generative model is $M$. Comparing the marginal likelihood values for different models allows for determining which model is the most likely to have generated the observations. In our case, consider weak prior information on the hyperparameters, the maximum a posteriori estimator simplifies to maximum likelihood: $\hat{\mathbf{\psi}}_{\mathrm{FMP}}(\boldsymbol{\theta}) = \arg\max_{\psi} p(\psi)p(\mathbf{y}_{\mathrm{obs}}|\boldsymbol{\theta}, \psi) \approx \arg\max_{\psi} p(\mathbf{y}_{\mathrm{obs}}|\boldsymbol{\theta}, \psi)$.

Assume that, in an alternative approximation method, the functional relationship $\hat{\mathbf{\psi}}_h(\boldsymbol{\theta}) = h(\boldsymbol{\theta})$ is proposed, with $h$ a generic function. The posterior of the model parameters would become $p_h(\boldsymbol{\theta}|\mathbf{y}_{\mathrm{obs}}) \propto p(\boldsymbol{\theta})p(\mathbf{y}_{\mathrm{obs}}|\boldsymbol{\theta}, \psi =$
$\psi_h(\theta)$). This setting also encompasses single-point estimation methods, such as the KOH or cross-validation estimator of the hyperparameters. We have:

$$ 
\begin{align*}
& p(\mathbf{y}_{\text{obs}}|\theta, \mathbf{\psi} = \hat{\mathbf{\psi}}_{\text{FMP}}(\theta)) \geq p(\mathbf{y}_{\text{obs}}|\theta, \mathbf{\psi} = \psi_h(\theta)), \quad \forall \theta \in \Theta \\
\Rightarrow \quad & \int_{\Theta} p(\theta)p(\mathbf{y}_{\text{obs}}|\theta, \mathbf{\psi} = \hat{\mathbf{\psi}}_{\text{FMP}}(\theta)) d\theta \geq \int_{\Theta} p(\theta)p(\mathbf{y}_{\text{obs}}|\theta, \mathbf{\psi} = \psi_h(\theta)) d\theta \\
\Leftrightarrow \quad & p(\mathbf{y}_{\text{obs}}|\mathbf{\psi} = \hat{\mathbf{\psi}}_{\text{FMP}}) \geq p(\mathbf{y}_{\text{obs}}|\mathbf{\psi} = \psi_h).
\end{align*}
$$

1 Thus, the likelihood of the observations when making the FMP approximation will always be higher than any other functional relationship between hyperparameters and parameters. The BMC will then systematically favor the FMP method against any other approach based on such a relationship or single-point estimation.

## 3.3 Elementary examples

5 To illustrate the behavior of the FMP calibration, we consider two examples. In the first one, the (Bayesian) joint posterior distribution of $(\theta, \psi)$ is Gaussian. In the second one, we assume it is a multimodal mixture of Gaussians. 6 The Bayes, KOH, and FMP calibration methods feature analytical solutions in these two cases.

8 A Gaussian posterior is unlikely in practice because priors on hyperparameters are rarely Gaussian but classically follow uniform or inverse Gamma distributions. However, if a single value for $\theta$ and $\psi$ perfectly fits the data, then the posterior is unimodal and a Gaussian approximation might be reasonable. However, the Gaussian case is very instructive, and posterior approximations by Gaussian distribution are standard approximation when the posteriors become increasingly sharp around their maximum as the number of observations increases. When several, well-separated, values of $\theta$ and $\psi$ explain the data equally well, the multimodal Gaussian assumption is more appropriate. 9 10 11 12 13 14 15

### 3.3.1 Unimodal Gaussian posterior

17 Let us assume that the Bayesian procedure leads to a joint Gaussian posterior for $\theta$ and $\psi$. We recall that $\dim(\theta) = p$ and $\dim(\psi) = h$. Let $\gamma = (\theta, \psi)$, $\mu = (\mu_\theta, \mu_\psi)$ is the posterior mean and $\Lambda = \begin{pmatrix} V_\theta & C_{\theta,\psi}^T \\ C_{\theta,\psi} & V_\psi \end{pmatrix}$ is the posterior
1 covariance, symmetric and positive definite. The joint posterior distribution writes:

$$p(\gamma|\mathbf{y}_{\text{obs}}) = \frac{1}{(2\pi)^{(p+h)/2}|\mathbf{\Lambda}|^{1/2}} \exp\left(-\frac{1}{2}(\gamma - \boldsymbol{\mu})^T \mathbf{\Lambda}^{-1} (\gamma - \boldsymbol{\mu})\right). \quad (21)$$

2 Following the properties of the multivariate Gaussian distribution, the marginal distributions are also Gaussian:

$$\theta|\mathbf{y}_{\text{obs}} \sim \mathcal{N}(\boldsymbol{\mu}_{\theta}, \mathbf{V}_{\theta}), \quad \text{and} \quad \psi|\mathbf{y}_{\text{obs}} \sim \mathcal{N}(\boldsymbol{\mu}_{\psi}, \mathbf{V}_{\psi}), \quad (22)$$

3 which directly gives:

$$p_{\text{Bayes}}(\theta) = \mathcal{N}(\boldsymbol{\mu}_{\theta}, \mathbf{V}_{\theta}). \quad (23)$$

4 Another property is that the conditional distributions are also Gaussian. We can thus obtain the distribution of parameters $\theta$ conditioned to a value of hyperparameters $\psi$:

5

$$\begin{align}
\theta|\psi, \mathbf{y}_{\text{obs}} &\sim \mathcal{N}(\boldsymbol{\mu}_{\theta|\psi}, \mathbf{V}_{\theta|\psi}), \nonumber \\
\boldsymbol{\mu}_{\theta|\psi} &= \boldsymbol{\mu}_{\theta} + \mathbf{C}_{\theta, \psi}^{T} \mathbf{V}_{\psi}^{-1} (\psi - \boldsymbol{\mu}_{\psi}), \tag{24} \\
\mathbf{V}_{\theta|\psi} &= \mathbf{V}_{\theta} - \mathbf{C}_{\theta, \psi}^{T} \mathbf{V}_{\psi}^{-1} \mathbf{C}_{\theta, \psi}. \nonumber
\end{align}$$

6 Accordingly, the distribution of hyperparameters $\psi$ conditioned to a value of parameters $\theta$ has for expression

$$\begin{align}
\psi|\theta, \mathbf{y}_{\text{obs}} &\sim \mathcal{N}(\boldsymbol{\mu}_{\psi|\theta}, \mathbf{V}_{\psi|\theta}), \nonumber \\
\boldsymbol{\mu}_{\psi|\theta} &= \boldsymbol{\mu}_{\psi} + \mathbf{C}_{\theta, \psi} \mathbf{V}_{\theta}^{-1} (\theta - \boldsymbol{\mu}_{\theta}), \tag{25} \\
\mathbf{V}_{\psi|\theta} &= \mathbf{V}_{\psi} - \mathbf{C}_{\theta, \psi} \mathbf{V}_{\theta}^{-1} \mathbf{C}_{\theta, \psi}^{T}. \nonumber
\end{align}$$

To perform the KOH estimation, the first step is to compute the marginal of hyperparameters and retrieve its maximum. From equation (22) we obtain

$$\hat{\boldsymbol{\psi}}_{\text{KOH}} = \boldsymbol{\mu}_{\psi}.$$

7 The approximation of the parameters' marginal is then obtained by conditioning $\theta$ on the estimated value $\hat{\boldsymbol{\psi}}_{\text{KOH}}$.

8 Thus,

$$p_{\text{KOH}}(\theta|\mathbf{y}_{\text{obs}}) = p(\theta|\hat{\boldsymbol{\psi}}_{\text{KOH}}, \mathbf{y}_{\text{obs}}) = p(\theta|\boldsymbol{\mu}_{\psi}, \mathbf{y}_{\text{obs}}) = \mathcal{N}(\boldsymbol{\mu}_{\theta|\boldsymbol{\mu}_{\psi}}, \mathbf{V}_{\theta|\psi}) = \mathcal{N}(\boldsymbol{\mu}_{\theta}, \mathbf{V}_{\theta|\psi}). \quad (26)$$
1 Note the reduction of variance property: because the matrix $C_{\theta,\psi}^T V_{\psi}^{-1} C_{\theta,\psi}$ is symmetric positive definite, under the
2 Loewner order we have $V_{\theta|\psi} \le V_{\theta}$. This order can be understood as the following: let $\{\lambda_i\}_{1\le i\le p}$ be the eigenvalues
3 of $V_{\theta|\psi}$ in a decreasing order and $\{\lambda'_i\}_{1\le i\le p}$ be the eigenvalues of $V_{\theta}$ in a decreasing order. Then, for all $i \le p$, we
4 have $\lambda_i \le \lambda'_i$. It is seen that $p_{\text{KOH}}$ remains Gaussian, with correct mean $\mu_{\theta}$, but a reduced covariance compared to
5 $p_{\text{Bayes}}$.

To proceed with the FMP estimation, the first step is to estimate the optimised hyperparameters for each param-
eter value $\theta$. We directly get, using equation (24),

$$
\hat{\psi}_{\text{FMP}}(\theta) = \arg \max_{\psi} p(\psi|\theta, y_{\text{obs}}) = \mu_{\psi|\theta}(\theta) = \mu_{\psi} + C_{\theta,\psi} V_{\theta}^{-1}(\theta - \mu_{\theta}).
$$

Thus $\hat{\psi}_{\text{FMP}}$ is an affine function of $\theta$. To evaluate the approximate posterior distribution of parameters we need to
express the inverse of the covariance matrix $\Lambda$, using the Matrix Block Inversion Lemma [33],

$$
\Lambda^{-1} = \begin{pmatrix}
V_{\theta}^{-1} + V_{\theta}^{-1} C_{\theta, \psi}^{T} V_{\psi|\theta}^{-1} C_{\theta, \psi} V_{\theta}^{-1} & -V_{\theta}^{-1} C_{\theta, \psi}^{T} V_{\psi|\theta}^{-1} \\
-V_{\psi|\theta}^{-1} C_{\theta, \psi} V_{\theta}^{-1} & V_{\psi|\theta}^{-1}
\end{pmatrix},
$$

6 we obtain

$$
\begin{align}
p_{\text{FMP}}(\theta|\mathbf{y}_{\text{obs}}) &\propto p(\theta, \mathbf{\psi} = \hat{\mathbf{\psi}}_{\text{FMP}}(\theta)|\mathbf{y}_{\text{obs}}) \nonumber \\
&\propto \exp\left(-\frac{1}{2}\left(\begin{pmatrix} \theta - \boldsymbol{\mu}_{\theta} \\ \hat{\mathbf{\psi}}_{\text{FMP}}(\theta) - \boldsymbol{\mu}_{\psi} \end{pmatrix}^{\mathrm{T}} \boldsymbol{\Lambda}^{-1} \begin{pmatrix} \theta - \boldsymbol{\mu}_{\theta} \\ \hat{\mathbf{\psi}}_{\text{FMP}}(\theta) - \boldsymbol{\mu}_{\psi} \end{pmatrix}\right)\right) \nonumber \\
&\propto \exp\left(-\frac{1}{2}\left(\begin{pmatrix} \theta - \boldsymbol{\mu}_{\theta} \\ C_{\theta,\psi} \mathbf{V}_{\theta}^{-1}(\theta - \boldsymbol{\mu}_{\theta}) \end{pmatrix}^{\mathrm{T}} \boldsymbol{\Lambda}^{-1} \begin{pmatrix} \theta - \boldsymbol{\mu}_{\theta} \\ C_{\theta,\psi} \mathbf{V}_{\theta}^{-1}(\theta - \boldsymbol{\mu}_{\theta}) \end{pmatrix}\right)\right) \tag{27} \\
&\propto \exp\left(-\frac{1}{2}(\theta - \boldsymbol{\mu}_{\theta})^{\mathrm{T}}\begin{pmatrix} 1_{p \times p} \\ C_{\theta,\psi} \mathbf{V}_{\theta}^{-1} \end{pmatrix}^{\mathrm{T}} \boldsymbol{\Lambda}^{-1} \begin{pmatrix} 1_{p \times p} \\ C_{\theta,\psi} \mathbf{V}_{\theta}^{-1} \end{pmatrix}(\theta - \boldsymbol{\mu}_{\theta})\right) \nonumber \\
&\propto \exp\left(-\frac{1}{2}(\theta - \boldsymbol{\mu}_{\theta})^{\mathrm{T}}\mathbf{V}_{\theta}^{-1}(\theta - \boldsymbol{\mu}_{\theta})\right). \nonumber
\end{align}
$$

7 The density $p_{\text{FMP}}$ is proportional to the exponential of a quadratic form in $\theta$, so it is Gaussian, and its mean and
8 covariance matrix correspond to the true marginal distribution. In this case, the FMP approximation is exact.

9 The behaviour of the methods is illustrated in figure 1 for $\dim(\theta) = \dim(\psi) = 1$. In this example, the FMP
method matches exactly the Bayesian solution. The KOH solution has a reduced variance due to the conditioning on a fixed value of the hyperparameters, which might lead to the "false certitude" effect. This higher certitude (lower uncertainty) is simply an artifact of the KOH estimation method.

<figure><img src="image_2.png" /><figcaption>Figure 1: Joint posteriors of parameters and hyperparameters: Gaussian case.</figcaption></figure>

### 3.3.2 Mixture of Gaussians

We now consider the case where the joint posterior $p(\gamma|\mathbf{y}_{\text{obs}})$ is a mixture of Gaussians, with modes that are well-separated, with similar orders of magnitude. The joint density is a mixture of $m$ Gaussians, with weights $(w_i)_{i \le m}$, such that $\sum_{i=1}^m w_i = 1$. Their respective means are $\boldsymbol{\mu}_i = \begin{pmatrix} \mu_{\theta,i} \\ \mu_{\psi,i} \end{pmatrix}$ and their respective covariance matrices: $\boldsymbol{\Lambda}_i = \begin{pmatrix} \mathbf{V}_{i,\theta} & \mathbf{C}_{i,\theta,\psi}^T \\ \mathbf{C}_{i,\theta,\psi} & \mathbf{V}_{i,\psi} \end{pmatrix}$. The joint density is

$$ \begin{aligned} p(\gamma|\mathbf{y}_{\text{obs}}) &= \sum_{i=1}^{m} w_i p(\gamma|\mathbf{y}_{\text{obs}}, \boldsymbol{\mu}_i, \boldsymbol{\Lambda}_i) \\ &= \frac{1}{(2\pi)^{(p+h)/2}} \sum_i \frac{w_i}{\sqrt{|\boldsymbol{\Lambda}_i|}} \exp\left(-\frac{1}{2}(\gamma - \boldsymbol{\mu}_i)^T \boldsymbol{\Lambda}_i^{-1} (\gamma - \boldsymbol{\mu}_i)\right). \end{aligned} \quad (28) $$

The well-separated hypothesis concerns the projections of the modes on the spaces $\Theta$ and $\Psi$. Geometrically, the well-separated hypothesis on $\Theta$ means that, when projecting the Gaussians on the first $p$ coordinates, the 95% confidence ellipses of each mode do not intersect each other. It can be written as:

> "The intervals $\Theta_j = \{\Theta \text{ s. t. } (\Theta - \boldsymbol{\mu}_j)^T \mathbf{V}_{\theta,j}^{-1} (\Theta - \boldsymbol{\mu}_j) \le t_{95}(p)\}$ are disjoint for $1 \le j \le m,$"
where we have noted $t_{95}(p)$ the 95% quantile of the $\chi_2$ law with $p$ degrees of freedom. The equivalent condition for the projections over $\Psi$ is also supposed to be true. Another assumption is that the weights $\{w_i\}_{1\le i\le m}$ have the same order of magnitude, so that, for $1 \le i \le m$, if $(\theta, \psi)$ is close to $\mu_i$ then $p(\theta, \psi|\mathbf{y}_{\text{obs}}) \approx w_i p(\theta, \psi|\mathbf{y}_{\text{obs}}, \mu_i, \Lambda_i)$. By linearity, the true marginal density of the parameters is the linear combination of marginals:

$$ \begin{aligned} p_{\text{Bayes}}(\theta|\mathbf{y}_{\text{obs}}) &= \int_{\Psi} p(\gamma|\mathbf{y}_{\text{obs}}) d\psi = \sum_{i=1}^{m} w_i \int_{\Psi} p(\gamma|\mathbf{y}_{\text{obs}}, \mu_i, \Lambda_i) d\psi \\ &= \frac{1}{(2\pi)^{p/2}} \sum_{i=1}^{m} \frac{w_i}{\sqrt{|\mathbf{V}_{i,\theta}|}} \exp \left( -\frac{1}{2} (\theta - \mu_{\theta,i})^T \mathbf{V}_{i,\theta}^{-1} (\theta - \mu_{\theta,i}) \right). \end{aligned} \quad (29) $$

The KOH estimation first computes the hyperparameter marginal distribution:

$$ p(\psi|\mathbf{y}_{\text{obs}}) = \frac{1}{(2\pi)^{h/2}} \sum_{i=1}^{m} \frac{w_i}{\sqrt{|\mathbf{V}_{i,\psi}|}} \exp \left( -\frac{1}{2} (\psi - \mu_{\psi,i})^T \mathbf{V}_{i,\psi}^{-1} (\psi - \mu_{\psi,i}) \right). $$

According to the separation hypothesis over $\Psi$, the maximum of $p(\psi|\mathbf{y}_{\text{obs}})$ is:

$$ \begin{aligned} \hat{\psi}_{\text{KOH}} &= \arg\max_{\psi} p(\psi|\mathbf{y}_{\text{obs}}) = \mu_{i_{\text{KOH}},\psi}, \\ i_{\text{KOH}} &= \arg\max_{i \le m} \frac{w_i}{\sqrt{V_{i,\psi}}}. \end{aligned} \quad (30) $$

So, the KOH posterior is:

$$ p_{\text{KOH}}(\theta|\mathbf{y}_{\text{obs}}) = p(\theta|\mathbf{y}_{\text{obs}}, \psi = \mu_{i_{\text{KOH}},\psi}) = \mathcal{N}(\mu_{\theta,i_{\text{KOH}}}, \mathbf{V}_{i_{\text{KOH}},\theta|\psi}), \quad (31) $$

with

$$ \mathbf{V}_{i_{\text{KOH}},\theta|\psi} = \mathbf{V}_{i_{\text{KOH}},\theta} - \mathbf{C}_{i_{\text{KOH}},\theta,\psi}^{T} \mathbf{V}_{i_{\text{KOH}},\psi} \mathbf{C}_{i_{\text{KOH}},\theta,\psi}. $$

The solution of the KOH estimation is Gaussian with reduced variance matrix. Besides, the selection of the mode is driven by the criteria defining $i_{\text{KOH}}$, which does not necessarily correspond to the true maximum of the Bayesian solution.

Optimal hyperparameters in the FMP method are given by the solution to the optimization problem:

$$ \hat{\psi}_{\text{FMP}}(\theta) = \arg\max_{\psi} p(\gamma|\mathbf{y}_{\text{obs}}) = \arg\max_{\psi} \sum_{i=1}^{m} \frac{w_i}{\sqrt{|\Lambda_i|}} \exp \left( -\frac{1}{2} (\gamma - \mu_i)^T \Lambda_i^{-1} (\gamma - \mu_i) \right) $$
According to the separation hypothesis, the confidence intervals {$\Theta_i$}$_{i \le m}$ are disjoint so that

$$
\text{for } \theta \in \Theta_i, \quad p(\gamma|\mathbf{y}_{\text{obs}}) \approx w_i p(\gamma|\mathbf{y}_{\text{obs}}, \boldsymbol{\mu}_i, \boldsymbol{\Lambda}_i),
$$

$$
\text{and } \hat{\psi}_{\text{FMP}}(\theta) = \hat{\psi}_{i,\text{FMP}}(\theta) := \mu_{i,\psi} + C_{i,\theta,\psi} \mathbf{V}_{i,\theta}^{-1}(\theta - \mu_{i,\theta}).
$$

1 We note $p_{\text{FMP}}^*(\theta) = p(\theta, \psi = \hat{\psi}_{\text{FMP}}(\theta)|\mathbf{y}_{\text{obs}}) = \sum_{j=1}^m w_j p(\theta, \hat{\psi}_{\text{FMP}}(\theta)|\mathbf{y}_{\text{obs}}, \boldsymbol{\mu}_j, \boldsymbol{\Lambda}_j)$ the unnormalized FMP approximation of the parameter posterior. The following two properties are true:

3. $\forall \theta \in \Theta_i$, $w_i p(\theta, \hat{\psi}_{i, \text{FMP}}(\theta) | \mathbf{y}_{\text{obs}}, \boldsymbol{\mu}_i, \boldsymbol{\Lambda}_i) \gg \sum_{j \neq i} w_j p(\theta, \hat{\psi}_{i, \text{FMP}}(\theta) | \mathbf{y}_{\text{obs}}, \boldsymbol{\mu}_j, \boldsymbol{\Lambda}_j)$. This is the application of the separation property in $\Psi$ and the low discrepancy in the weights {$w_i$}$_{1 \le i \le m}$.

5. $w_i p(\theta, \hat{\psi}_{i, \text{FMP}}(\theta) | \mathbf{y}_{\text{obs}}, \boldsymbol{\mu}_i, \boldsymbol{\Lambda}_i) = \frac{w_i}{(2\pi)^{(p+h)/2} \sqrt{|\boldsymbol{\Lambda}_i|}} \exp \left( -\frac{1}{2} (\theta - \boldsymbol{\mu}_{i,\theta})^T \mathbf{V}_{i,\theta}^{-1} (\theta - \boldsymbol{\mu}_{i,\theta}) \right)$. This result was obtained in the unimodal case (see equation (27)).

Thus, a natural approximation for $p_{\text{FMP}}^{*}$, for $\theta \in \Theta$, is

$$
P_{\text{FMP}}^{*}(\theta) \approx \sum_{j=1}^{m} w_j p(\theta, \hat{\psi}_{j, \text{FMP}}(\theta) | \mathbf{y}_{\text{obs}}, \mu_j, \Lambda_j),
$$

and if $\theta \in \Theta_i$ it simplifies to

$$
P_{\text{FMP}}^{*}(\theta) \approx \frac{w_i}{(2\pi)^{(p+h)/2} \sqrt{|\Lambda_i|}} \exp \left( -\frac{1}{2} (\theta - \mu_{i,\theta})^T \mathbf{V}_{i,\theta}^{-1} (\theta - \mu_{i,\theta}) \right).
$$

7. Since most of the probability mass of $p_{\text{FMP}}^{*}(\theta|\mathbf{y}_{\text{obs}})$ is contained within the intervals $\Theta_i$, we have

$$
P_{\text{FMP}}(\theta) = \frac{1}{K} \sum_{i=1}^{m} \frac{w_i}{\sqrt{|\Lambda_i|}} \exp \left( -\frac{1}{2} (\theta - \mu_{i,\theta})^T \mathbf{V}_{i,\theta}^{-1} (\theta - \mu_{i,\theta}) \right) \quad \text{for } \theta \in \Theta, \quad (32)
$$

with the normalizing constant

$$
\begin{align*}
K &= \int_{\theta \in \Theta} \sum_{i=1}^{m} \frac{w_i}{\sqrt{|\Lambda_i|}} \exp \left( -\frac{1}{2} (\theta - \mu_{i,\theta})^T \mathbf{V}_{i,\theta}^{-1} (\theta - \mu_{i,\theta}) \right) d\theta \\
&\approx (2\pi)^{(p/2)} \sum_{i=1}^{m} w_i \frac{\sqrt{|\mathbf{V}_{i,\theta}|}}{\sqrt{|\Lambda_i|}}.
\end{align*}
$$

8. The distribution $p_{\text{FMP}}$ is a linear combination of the $m$ marginals associated with the original Gaussian mixture. The

9. variances of the different peaks are correctly estimated, but the estimation introduces a bias on the weights of the
1 Gaussian mixture: the FMP estimates the original weights $w_i$ to be $w'_i = \frac{1}{k} w_i \frac{\sqrt{|V_{i,o}|}}{\sqrt{|\Lambda_i|}}$ where $k$ is a normalizing constant such that $\sum_{i=1}^m w'_i = 1$. Specifically, the FMP method emphasizes the importance of the Gaussians with low hyperparameter variance in the mixture.

4 An illustration of the methods is given in figure 2. In this example, the KOH estimation finds only one peak that does not correspond to the true maximum of the posterior, so the conclusions about $\Theta$ are potentially misleading. 5 In the FMP estimation, all the peaks are found, yet their relative weights might be wrongly estimated, leading to an 6 inversion of the importance of the peaks so that the maximum a posteriori estimator is incorrect. This feature will be 7 found in any method that uses estimated hyperparameters instead of marginalization (see assumption 2 in equation 16) 8 because volume effects cannot be seen. We argue that it is still favorable to perform the FMP estimation in this case 9 so that no possible explanation of the observations is missed. 10

<figure><img src="image_3.png" /><figcaption>**Figure 2:** Joint posterior of parameters and hyperparameters: well separated Gaussian mixture case.</figcaption></figure>

11 We now discuss the case where the two Gaussians of the mixture are not well-separated. In figure 3a, the $\Psi$-projections of the two modes overlaps. The corresponding estimations of the KOH and FMP methods are plotted in the figure. In this example, the KOH criterion selects the right mode; contrary to the well-separated case, the other (left) mode contributes to the parameter posterior because part of its probability mass is captured when conditioning on $\hat{\psi}_{\text{KOH}}$. However, KOH badly estimates the second mode, contrary to the FMP approximation, which retrieves the two modes correctly. This situation occurs when a single model discrepancy adequately applies to the whole parameters domain.

18 Figure 3b shows the case where the $\Theta$-projections of the Gaussians overlap. This situation is challenging for the two approximation methods since they rely on a point-mass approximation of the marginals of $\psi$, which is not verified here. In the example shown, the FMP approximation performs slightly better with a posterior variance of $\Theta$ underestimated but closer to the reference than for KOH. This situation occurs when the distribution of the model discrepancy is at odds with the observations, and no hyperparameters value stands out; it calls for a change in the
1 prior of $z_{\theta}$ (e.g., selecting another covariance structure).

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Comparison of KOH and FMP estimations: case of Gaussian mixtures without modes separation.</figcaption></figure>

2 **4. APPLICATIONS**

3 In this section, we present two examples. The first one in Section 4.1 is an analytical model in which predictions
4 have different shapes depending on the parameter value, yielding two explanations of the observations. The posterior
5 is then bimodal, challenging the KOH estimation. The elementary nature of this example also allows studying the
6 normality of the posterior by fitting a Gaussian Mixture Model.

7 The second application in Section 4.2 deals with the calibration using actual experimental measurements of a
8 boiling model. Physical insights are required to formulate the statistical assumptions. The model error and the mea-
9 surement uncertainty are both significant in this problem, and the FMP calibration correctly attributes the uncertainty
10 of each source, whereas the KOH calibration sees only one.

11 **4.1 Calibration of an inadequate model**

12 We apply the FMP calibration to an elementary problem where model inadequacy is present, and the structure of
13 the model predictions is sensitive to parameter variations. We also perform the KOH and full Bayesian calibrations
14 on this example. The Bayesian calibration shows that the joint posterior can be approximated by a mixture of two
15 Gaussian distributions, illustrating the previous theoretical results.
### 1 4.1.1 Problem formulation

2 In this example, the true function is $y(x) = x$ and the computer model is given by

$$f(x, \theta) = x \sin(2\theta x) + (x + 0.15)(1 - \theta). \quad (33)$$

3 The input variable $x$ and the model parameter $\theta$ lie within the restricted range $(x, \theta) \in [0, 1] \times [-0.5, 1.5]$. We use 4 observations of the true function at 8 points uniformly spaced in the interval $[0, 1]$. The observations are noisy with a 5 centered Gaussian noise with standard deviation $\sigma_e = 0.1$. The observation noise is not known a priori but is learned 6 along with the other hyperparameters of the model error.

7 The prior of the model discrepancy is a Gaussian Process with zero mean and a Squared Exponential covariance

$$c_{\psi}(x, x') = \sigma^2 \exp \left( - \frac{(x - x')^2}{2l^2} \right). \quad (34)$$

8 Thus, a total of three hyperparameters $\psi = (\sigma, l, \sigma_e)$ are involved in the calibration. The prior for $\theta$ is uniform over 9 $[-0.5, 1.5]$. The priors for $\sigma$ and $\sigma_e$ are uniform on the interval $[1e-5, 1]$. The prior on $l$ is also uniform, truncated 10 on $[1e-3, 5]$.

11 Samples of the posterior distributions are obtained using Monte-Carlo Markov Chain methods with the Metropolis-Hastings algorithm [48]. For the KOH calibration, we first estimate the hyperparameters using (11). The MCMC 12 method is then run for the target density $p_{\text{KOH}}(\theta)$ (see (12)). In the FMP calibration, the target density $p_{\text{FMP}}(\theta)$ is 13 also of dimension 1, but each step of the MCMC requires the determination of the optimal hyperparameters. The cost 14 of these optimizations is not prohibitive, thanks to the available gradients [33, Chapter 5], and because the optimal 15 hyperparameters of the previous step are generally a good starting point. For the Bayes calibration, the target density 16 $p(\theta, \psi|y_{\text{obs}})$ is in dimension four and requires significantly more steps. The FMP and KOH chains have 5e5 steps, 17 and the Bayes chain has 2.5e6 steps, all with a burn-in of 10%. From each chain, we extract a posterior sample of size 18 19 5000 by taking regularly-spaced visited states. This subsampling rate is large enough compared to the self-correlation 20 length of the chains to consider the resulting samples independent. For more details about this procedure, see [49].

### 2 4.1.2 Calibration Results

22 We apply a Kernel Density Estimation method to the MCMC samples to estimate the posterior marginal densities 23 of parameters and hyperparameters. These marginals are shown in figure 4. The parameter's posterior marginals are
bimodal with a first mode located at $\theta = -0.025$ and the second at $\theta = 1.02$. For the reference Bayes solutions, the relative importance of the two modes is in the ratio 4.3:1, computed as the relative frequency of $\theta < 0.5$. The ratio is 2.75:1 for the FMP method and 69.4:1 for the KOH method. The FMP method overestimates the second mode's importance and puts more probability mass between the two modes. The KOH estimation almost misses the second mode focusing on a single interpretation of the observations.

Hyperparameters' posterior marginals reveal that the FMP method improves the calibration results. All three hyperparameters feature a variance a posteriori that the KOH estimation cannot capture. The KOH hyperparameters, by definition, are at the Bayes posterior maxima; they match the maxima of the marginals of $\sigma_e$ and $\sigma$, although with a slight offset for $l$. The FMP marginals have modes located at these maxima and have smaller variances than the Bayes solution due to the reduction in their optimal values. Note that the posterior marginal of $l$ accumulates mass around the maximal value a priori. We have repeated the calibration with high maximal values up to $l = 20$ without finding a notable difference, as the Bayes posterior remains flat. Note that for the inferred $\sigma_e$, the Bayes solution is the only one to correctly estimate the true value of 0.1 when other methods infer smaller values $\sigma_e \approx 0.06$.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Prior and posterior distributions for the model parameter (top left) and hyperparameters. In the KOH calibration, hyperparameters are estimated with point mass distributions, represented with vertical lines. The Gaussian Mixture Model, fitted on the Bayes posterior, is represented with the dashed curves. In the bottom row, some probability mass lies outside the support of the densities due to the Kernel Density Estimation. All densities are normalized over their support for proper comparison.</figcaption></figure>

The bimodality of the parameter posterior reflects in the posterior model predictions shown in figure 5. The left mode, the only one found by KOH, corresponds to quasi-linear predictions with a constant offset from the true
process $y(x)$. The second modes correspond to predictions close to the true process but with an oscillation in $x$. The FMP method accounts for the two types of predictions present in the Bayes solution. We also see a clear separation between the two types of predictions in the Bayes solutions, while the FMP method produces some predictions with intermediate structures.

<figure><img src="image_6.png" /><figcaption>Figure 5: Posterior samples of the model predictions for all three calibration methods. 50 samples are represented on each figure, and observations are shown with $2\sigma_e$ confidence intervals. Model predictions that are quasi-linear correspond to $\theta \approx 0$, and the ones with slow oscillations correspond to $\theta \approx 1$.</figcaption></figure>

### 4.1.3 Gaussian Mixture Model fit on the Bayes solution

As a further step, we fit a mixture of Gaussians model on the MCMC samples of the Bayes solution. To do so, we employed a hard clustering algorithm based on Expectation-Maximization *. The selected form of the mixture is a weighted sum of two Gaussians

$$p(\theta, \psi|\mathbf{y}_{\text{obs}}) = w_1 N(\mu_1, \mathbf{C}_1) + w_2 N(\mu_2, \mathbf{C}_2). \qquad (35)$$

The estimated coefficients of the GMM are presented in Table 1.

The two estimated modes reveal that, even though they are distinct in $\theta$, they are not in $\psi$. We are thus in a situation without clear separation featuring the intersection of projections in $\psi$ space. The marginals of the GMM are shown in Figure 4. A good fit of the Bayesian and GMM marginals is reported for $\theta$. The fit is much less satisfactory for the hyperparameters: the normality assumption is inappropriate due to the truncated supports. Nevertheless, the GMM illustrates some results of section 3.3.2 as discussed below:

* The KOH estimation selects a mode according to the criterion in equation (30). We check that $\frac{w_1}{\sqrt{|\mathbf{C}_1^{[3,3]}|}} = 70.8$ is greater than $\frac{w_2}{\sqrt{|\mathbf{C}_2^{[3,3]}|}} = 10.0$, which is the reason why it selects the first mode.

<footnote>*https://perso.telecom-paristech.fr/bonald/documents/gmm.pdf</footnote>
<table><caption>**Table 1:** Estimated coefficients of the Gaussian Mixture Model for the posterior distribution of $\mu = (\theta, \sigma_{\text{mes},\Delta T}, \sigma, l)$.</caption><thead><tr><th>Coefficient</th><th>Estimation</th></tr></thead><tbody><tr><td>w<sub>1</sub></td><td>0.81</td></tr><tr><td>w<sub>2</sub></td><td>0.19</td></tr><tr><td>μ<sub>1</sub></td><td>(-0.03, 0.09, 0.40, 2.93)<sup>T</sup></td></tr><tr><td>μ<sub>2</sub></td><td>(1.02, 0.12, 0.29, 2.56)<sup>T</sup></td></tr><tr><td>C<sub>1</sub></td><td>(<table><tbody><tr><td>2.9e-2</td><td>-1.9e-4</td><td>1.3e-3</td><td>-1.6e-3</td></tr><tr><td>-1.9e-4</td><td>1.5e-3</td><td>-8.6e-5</td><td>-1.4e-3</td></tr><tr><td>1.3e-3</td><td>-8.6e-5</td><td>5.7e-2</td><td>3.0e-2</td></tr><tr><td>-1.6e-3</td><td>-1.4e-3</td><td>3.0e-2</td><td>1.6</td></tr></tbody></table>)</td></tr><tr><td>C<sub>2</sub></td><td>(<table><tbody><tr><td>2.1e-2</td><td>-3.0e-4</td><td>-8.0e-4</td><td>5.1e-2</td></tr><tr><td>-3.0e-4</td><td>2.9e-3</td><td>5.7e-4</td><td>7.8e-3</td></tr><tr><td>-8.0e-3</td><td>5.7e-4</td><td>6.2e-2</td><td>4.9e-2</td></tr><tr><td>5.1e-2</td><td>7.8e-3</td><td>4.9e-2</td><td>2.1</td></tr></tbody></table>)</td></tr></tbody></table>

• The posterior correlation between $\theta$ and $\psi$ is low, so the variance at the first mode is roughly unchanged after projection: $V_{\theta} = C_{1}(0,0) = 2.90e-2$, and $V_{\theta|\psi} = C_{1}(0,0) - C_{1,\theta,\psi}^{T}(C_{1}^{[3,3]})^{-1}C_{1,\theta,\psi} = 2.89e-2$, with $C_{1,\theta,\psi} = (-1.9e-4, 1.3e3, -1.6e-3)^T$. The KOH method underestimates the posterior variance due to selecting only one mode.

• The FMP weights can be computed using (32). We have $w'_{i} \propto w_{i} \frac{\sqrt{|C_{i}|}}{\sqrt{C_{i}(0,0)}}$, which gives $w'_{1} = 0.73$ and $w'_{2} = 0.27$, in a proportion 2.71:1 which corresponds to the ratio found in the FMP results.

Thus, even if the assumption of normal joint posterior is questionable, this example shows the key role of the posterior covariance in understanding the behaviour of the calibration methods. However, this covariance can not be accessed before the calibration.

To further contrast the quality of the FMP and KOH distributions, we measure their respective similarities with the Bayes distribution using the Jensen-Shannon divergence ($D_{JSD}$). The JS-Divergence between two distributions with density $p_{*}$ and the Bayes' density $p_{\text{Bayes}}$ is given by

$$D_{JS}(p_{*} \| p_{\text{Bayes}}) = \frac{1}{2} D_{\text{KL}}(p_{*}, p_{1/2}) + \frac{1}{2} D_{\text{KL}}(p_{\text{Bayes}}, p_{1/2}),$$

where $p_{1/2} = (p_{*} + p_{\text{Bayes}})/2$ and $D_{\text{KL}}(p_{1} \| p_{2})$ is the standard Kullbach-Leibler divergence. The JS-Divergence is symmetric and is zero when the two distributions compared are equal. A larger value of $D_{JS}$ means a higher dissimilarity. The present work estimates the divergences using the MCMC sample sets. In the first stage, we use a KDE method to estimate the densities. In the second stage, the densities are evaluated at the sample set's points to estimate the averages in the expression of the KL divergences. Compared to the KL divergence, the JS divergence is
1 much less susceptible to sampling noise and more stable when the densities are dissimilar. In the present analysis, 2 resampling with bootstrapping has shown that the estimates are sufficiently converged to characterize the differences 3 between the two methods.

4 The divergence for the posterior distributions of the parameter $\theta$ and the model's prediction $f(\cdot, \theta)$ at some $x$ are 5 reported in Table 2. These results show that the FMP distribution has a much lower JS-Divergence with the Bayesian 6 distribution than the KOH distribution. Specifically, in this example, the divergence from the Bayesian distribution of 7 the FMP method is less than 7.5% than that of the FMP method.

<table><caption>**Table 2:** JS-Divergences between the Bayes and FMP and the Bayes and KOH distributions. The divergences are computed for the posterior distributions of the parameter ($\theta$) and model's predictions ($f(\cdot, \theta)$) at different values of $x$ as indicated. The last column provides the ratio of the Bayes-FMP and Bayes-KOH divergences.</caption><thead><tr><th rowspan="2">Quantity</th><th colspan="2">JS Divergences</th><th rowspan="2">Relative (%)</th></tr><tr><th>Bayes-FMP</th><th>Bayes-KOH</th></tr></thead><tbody><tr><td>$\theta$</td><td>0.0072</td><td>0.1416</td><td>5.1</td></tr><tr><td>$f(x = 0.56, \theta)$</td><td>0.0052</td><td>0.1222</td><td>4.2</td></tr><tr><td>$f(x = 0.69, \theta)$</td><td>0.0074</td><td>0.1334</td><td>5.5</td></tr><tr><td>$f(x = 0.94, \theta)$</td><td>0.0097</td><td>0.1295</td><td>7.5</td></tr></tbody></table>

## 4.2 Calibration of a boiling model

9 We now tackle the calibration of a thermal model using experimental observations. Three parameters are selected for 10 the calibration for their strongly non-linear impact on the model predictions. The model discrepancy distribution is 11 formulated using physical knowledge.

### 4.2.1 Problem formulation

13 The MIT Boiling model [50] simulates the boiling of fluids in contact with a heated wall. The inputs of the model 14 are experimental variables (fluid pressure, velocity, temperature, and configuration geometry) and the wall superheat 15 $\Delta T_{sup}$, defined as the difference between the wall temperature and the saturation temperature of the fluid. The model 16 output is the heat flux $\phi$ from the wall to the liquid. The model is semi-empirical and based on experimental correla- 17 tions, with low evaluation cost. The experimental observations used in this study come from the boiling experiment of 18 Kennel [51], precisely case number 6 (one of the observation sets used by Kommajosyula to validate the MIT Boiling 19 model). It consists of 8 joint measurements of heat flux and wall temperature in different boiling regimes, ranging 20 from no-boiling to fully developed nucleate boiling.

21 Three calibration parameters ($\theta_1, \theta_2, \theta_3$) are considered; the motivations for this choice and the significance of
each parameter is discussed in [52]. The measurement uncertainty lies primarily on $\Delta T_{sup}$, which is an input model quantity, not an output. This specificity requires an adaptation of the framework. Following [53], the sensitivity (derivatives) of $y_{obs}$ at the uncertain inputs is estimated from a polynomial fit over the observations and used to include this uncertainty in the framework (see below).

We base the formulation of the model discrepancy on two considerations. First, because $\phi$ is a smooth increasing function of $\Delta T_{sup}$, the model discrepancy should be smooth too. The second consideration concerns the magnitude of the discrepancy. At low values of $\Delta T_{sup}$, there is no boiling, and the model predictions are accurate (linear tendency). At high values of $\Delta T_{sup}$, the emergence of nucleate boiling induces an exponential increase of the heat flux that might be incorrectly represented by the model, requiring a higher variability of the discrepancy term. Consequently, the selected model discrepancy is a Gaussian Process with mean zero and covariance function given by:

$$c_{\psi}(\Delta T_{sup}, \Delta T'_{sup}) = \sigma^2 \Delta T_{sup} \Delta T'_{sup} \exp \left( - \frac{|\Delta T_{sup} - \Delta T'_{sup}|^2}{2l^2} \right). \quad (36)$$

The vector of observations has a normal distribution, with covariance matrix $\Sigma_{\psi} + \sigma_{mes, \Delta T}^2 D$, where D is the diagonal matrix with diagonal coefficients equal to the derivatives of $y_{obs}$ at the observations.

Table 3 reports the prior distributions considered for the calibration. The priors are truncated normals and uniform distributions for the parameters and hyperparameters, respectively. The priors' ranges are selected to ensure physically meaningful value (parameters) and using a priori considerations on the magnitude of the model and measurement errors (hyperparameters).

As in the first example, we sample the posteriors with the Metropolis-Hastings algorithm. The KOH and FMP chains have an acceptance rate of ≈ 11% and the Bayes chain ≈ 3%. Figure 6 shows the self-correlation of the MCMC chains. The mixing lengths of the chains are $\tau = 28$ for KOH, $\tau = 30$ for FMP, and $\tau = 470$ for Bayes. To ensure sample sets of similar quality, we generate chains with respective lengths 5e5, 5e5, and 1e7 for the KOH, FMP, and Bayes methods and extract 1,000 samples from the chains by uniform subsampling.

A polynomial fit of the model is used to estimate derivatives. The FMP optimizations use the C++ library NLopt†, with the global algorithm MLSL, using the local optimizer SBPLX. A maximal time of 10e-4s is allocated to individual optimization. For the KOH estimation of hyperparameters, the integral over $\Theta$ is computed by quadrature, using a QMC grid of size 300 over the parameters domain.

<footnote>†https://nlopt.readthedocs.io/en/latest/</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 6: Self-correlation of the chains for the three methods applied to the MIT Boiling model calibration.</figcaption></figure>

## 1 4.2.2 Results

Table 3 summarizes the results of the inference problem for the three methods. The most noticeable result is that the KOH method estimates a zero model error $\sigma$ and attributes all the discrepancies to measurement uncertainty. Consistently, the KOH method yields an erroneous estimation of $\theta_1$ with an underestimated variance. On the contrary, the FMP method identifies a model error comparable to the full Bayesian calibration and, despite underestimating the hyperparameters' variances, provides a globally correct inference of the parameters and hyperparameters.

<table><caption>**Table 3:** Prior distributions and posterior summaries obtained for each calibration technique. Notations are $\mathcal{U}$ for uniform distributions, and $N_{t,I_k}(\mu, \sigma^2)$ refers to a truncated normal distribution over the interval $I_k$, with $I_1 = [-0.5, 1.75]$ and $I_2 = [0, 1.5]$.</caption><thead><tr><th>Calibration</th><th>θ<sub>1</sub></th><th>θ<sub>2</sub></th><th>θ<sub>3</sub></th><th>σ</th><th>l</th><th>σ<sub>mes,ΔT</sub></th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-weight: bold;">Prior distributions</td></tr><tr><td>-</td><td>N<sub>t,I<sub>1</sub></sub>(0.5, 0.3<sup>2</sup>)</td><td>N<sub>t,I<sub>2</sub></sub>(0.5, 0.3<sup>2</sup>)</td><td>N<sub>t,I<sub>2</sub></sub>(0.5, 0.3<sup>2</sup>)</td><td>U([0, 1e6])</td><td>U([1, 30])</td><td>U([0.1, 0.8])</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Posterior summaries (mean ± std)</td></tr><tr><td>KOH</td><td>0.45 ± 0.22</td><td>0.50 ± 0.28</td><td>0.37 ± 0.12</td><td>0 ± 0</td><td>15.0 ± 0</td><td>0.48 ± 0</td></tr><tr><td>FMP</td><td>0.28 ± 0.29</td><td>0.53 ± 0.28</td><td>0.41 ± 0.14</td><td>5.7e3 ± 8.6e3</td><td>13.2 ± 5.5</td><td>0.32 ± 0.14</td></tr><tr><td>Bayes</td><td>0.27 ± 0.29</td><td>0.48 ± 0.28</td><td>0.43 ± 0.18</td><td>5.6e4 ± 1.0e5</td><td>17.2 ± 7.8</td><td>0.37 ± 0.17</td></tr></tbody></table>

Figure 7 shows the posterior of the corrected model prediction ($f + z_{\theta}$) of the true process $y$. The credible intervals of the corrected prediction in FMP and Bayes are generally larger than for KOH, especially at low values of $\Delta T_{sup}$ where no observations are available. Note that the uncertainty is zero at $\Delta T_{sup}$ because of the structure of the model discrepancy covariance. The predictive variance is split into the model and residual contributions, as described in Section 3.1.2. For KOH, the variance a posteriori of $z_{\theta}$ is uniformly zero, as $\sigma$ is inferred to be zero. This result is unsatisfactory, especially at low values of $\Delta T_{sup}$ where few observations are available. All approaches agree on the dominance of the model predictions variance for the high range of $\Delta T_{sup}$ values.

More insight can be gained by plotting the projection of the posterior samples in the $\sigma_{mes,\Delta T}$ and $l$ plane; see
<figure><img src="image_8.png" /><figcaption>**Figure 7:** Posterior predictions of the corrected model ($f + z_{\theta}$) with uncertainty. Observations are shown with inferred error bars of length $2\mathbb{E}_{\theta}[\sigma_{mes, \Delta T}]$. The bottom plots show the decomposition of the predictive variance in residual uncertainty (stripped region) and uncertainty in the corrected model (grey filled region) as explained in Section 3.1.2.</figcaption></figure>

1. the left plot of Fig. 8. The Bayes samples cover the full support of the prior; the FMP samples fall primarily in two regions corresponding to two distinct interpretations of the observations: high measurement error with zero model error or a combination of both. The KOH estimator falls into the former. Since KOH estimates $\sigma = 0$, all lengths $l$ are equally likely, and the reported value $l = 15.0$ comes from the initialization of the optimization problem. The second and third plots of Fig. 8 report the FMP model predictions corresponding to these two regions. In the interpretation with measurement error only, plausible model predictions pass through all observations and have tight dispersions.
7. In the interpretation with non-zero model error, some model predictions come close to most observations, but others are further away, following the observations' trend. From a practical perspective, acknowledging the possibility of low measurement and non-zero model errors is crucial to properly assess the uncertainty in the model predictions of non-observed quantities that can not be corrected. The FMP and Bayes methods achieve this goal of considering alternative interpretations of the observations.
9.
10.
11.

<figure><img src="image_9.png" /><figcaption>**Figure 8:** Samples in the $\sigma_e$ and $l$ plane for the three methods (left plot). Samples of the FMP model predictions conditioned on a zero (middle) and non-zero model error (right).</figcaption></figure>
1

# 5. CONCLUSIONS AND DISCUSSION

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

We have proposed an approach to estimate both model parameters and model discrepancy. For the first time in literature, to our knowledge, the model error is explicitly dependent on model parameter values in the calibration framework. We showed that it solves the identifiability problem between model error and parameter uncertainty. The resulting predictive uncertainty splits into two natural contributions: the error in the calibrated model and the residual uncertainty. The FMP approximation has a reduced cost compared to the full Bayesian inference while correcting the defects of the KOH estimations when the posterior of the hyper-parameters has several modes.

In two applications, we found that the FMP calibration performs better than the KOH calibration in terms of uncertainty estimation, values of parameters, and posterior predictions. Allowing some variance in the hyperparameters of the model discrepancy helps avoid pitfalls such as the false certainty effect or missing entire probability regions. The FMP approach proves to be an accurate approximation of the full Bayesian calibration and significantly reduces the dimension of the MCMC sample space. Thus, the FMP method is well suited for situations requiring complex model error terms with many hyperparameters or when performing calibration on a collection of experiments with independent model discrepancy terms.

The reduction of the sample space in the FMP method is achieved by using an optimized value for the hyper-parameters, rather than sampling its conditional distribution. The FMP method then solves an optimization problem for the pointwise estimation of the hyperparameters, at each step of the MCMC, which induces a significant computational overhead compared to using a global pointwise estimates. Specifically, we solve the optimization problem in equation (15), where each evaluation of the likelihood requires inverting the covariance matrix of the experimental data. This is the most costly step of the FMP procedure since this operation is repeated for each iteration of the MCMC algorithm, especially when the number of observations is high or there are many model error hyperparameters. Because of the relatively low dimensionality of the two application problems presented in this paper, we have not observed significant computational benefits of using the FMP method compared to the full Bayesian approach in our numerical experiments. We are currently developing adaptive methods that construct surrogate models of the FMP hyperparameters $\hat{\psi}_{\text{FMP}}(\theta)$ in an offline stage. These surrogates bypass the optimization steps with considerable computational savings compared to the full Bayesian method, without significantly altering the quality of posterior approximation. With the hyperparameters surrogates, the cost of sampling the FMP posterior is comparable to the cost of the KOH method. These results will be presented in an article in preparation.

Our calibration equation offers new challenges to propose appropriate priors. Consistently with their definition of model error, previous works used model discrepancy hyperparameters a priori independent of model parameters.
1 We also adopt this hypothesis in the present work, but one could imagine more complex prior. We advocate for
2 using model discrepancy distributions that are more complex and more physics-informed to provide more accurate
3 predictions.

4 **SUPPLEMENTARY MATERIAL**

5 The implementation of the MIT Boiling model, as well as the Kennel measurements data set used in section 4.2, are
6 available online at https://github.com/nleoni95/MITB.

7 **ACKNOWLEDGMENTS**

8 This work was partly funded by the French Innovation Agency for Defense (Agence de l'Innovation de Défense,
9 AID). N. Leoni would also like to thank Ravikishore Kommajosyula for providing the boiling model, and Kathryn
10 Maupin for helpful discussions.

11 This project was partly funded from the Clean Sky 2 Joint Undertaking (JU) under grant agreement No 101008257.

12 **References**

13 1. Oberkampf, W.L. and Roy, C.J., *Verification and Validation in Scientific Computing*, Cambridge University Press, Cambridge, 2010.

15 2. Jaynes, E.T. and Bretthorst, G.L., *Probability Theory: The Logic of Science*, Cambridge University Press, Cambridge, UK; New York, NY, 2003.

17 3. Kennedy, M.C. and O'Hagan, A., Bayesian calibration of computer models, *Journal of the Royal Statistical Society: Series B* (*Statistical Methodology*), 63(3):425-464, 2001.

19 4. Edeling, W., Cinnella, P., Dwight, R., and Bijl, H., Bayesian estimates of parameter variability in the k-ε turbulence model, *Journal of Computational Physics*, 258:73-94, February 2014.

21 5. Xiao, H. and Cinnella, P., Quantification of model uncertainty in RANS simulations: A review, *Progress in Aerospace Sciences*, 108:1-31, July 2019.

23 6. Nadiga, B., Jiang, C., and Livescu, D., Leveraging Bayesian analysis to improve accuracy of approximate models, *Journal of Computational Physics*, 394:280-297, October 2019.

25 7. Doronina, O.A., Murman, S.M., and Hamlington, P.E., Parameter Estimation for RANS Models Using Approximate Bayesian Computation, arXiv:2011.01231 [physics], November 2020.
1

8. Damblin, G., Keller, M., Barbillon, P., Pasanisi, A., and Parent, É., Bayesian Model Selection for the Validation of Computer Codes, *Quality and Reliability Engineering International*, 32(6):2043-2054, 2016.

2

3

9. del Val, A., Le Maître, O.P., Magin, T.E., Chazot, O., and Congedo, P.M., A surrogate-based optimal likelihood function for the bayesian calibration of catalytic recombination in atmospheric entry protection materials, *Applied Mathematical Modelling*, 101:791-810, 2022.

4

5

6

10. Gally, T., Groche, P., Hoppe, F., Kuttich, A., Matei, A., Pfetsch, M.E., Rakowitsch, M., and Ulbrich, S., Identification of Model Uncertainty via Optimal Design of Experiments applied to a Mechanical Press, arXiv:1910.08408 [cs, math, stat], February 2020.

7

8

9

11. Maupin, K.A. and Swiler, L.P., Model Discrepancy Calibration Across Experimental Settings, *Reliability Engineering & System Safety*, p. 106818, January 2020.

10

11

12. Higdon, D., Kennedy, M., Cavendish, J.C., Cafeo, J.A., and Ryne, R.D., Combining Field Data and Computer Simulations for Calibration and Prediction, *SIAM Journal on Scientific Computing*, 26(2):448-466, January 2004.

12

13

13. Higdon, D., Gattiker, J., Williams, B., and Rightley, M., Computer Model Calibration Using High-Dimensional Output, *Journal of the American Statistical Association*, 103(482):570-583, 2008.

14

15

14. Arendt, P.D., Apley, D.W., Chen, W., Lamb, D., and Gorsich, D., Improving Identifiability in Model Calibration Using Multiple Responses, *Journal of Mechanical Design*, 134(10), October 2012.

16

17

15. Bayarri, M.J., Berger, J.O., Paulo, R., Sacks, J., Cafeo, J.A., Cavendish, J., Lin, C.H., and Tu, J., A Framework for Validation of Computer Models, *Technometrics*, 49(2):138-154, May 2007.

18

19

16. Bayarri, M.J., Berger, J.O., Cafeo, J., Garcia-Donato, G., Liu, F., Palomo, J., Parthasarathy, R.J., Paulo, R., Sacks, J., and Walsh, D., Computer model validation with functional output, *The Annals of Statistics*, 35(5):1874-1906, October 2007.

20

21

17. Brynjarsdóttir, J. and O'Hagan, A., Learning about physical parameters: The importance of model discrepancy, *Inverse Problems*, 30(11):114007, October 2014.

22

23

18. Ling, Y., Mullins, J., and Mahadevan, S., Selection of model discrepancy priors in Bayesian calibration, *Journal of Computational Physics*, 276:665-680, November 2014.

24

25

19. Gardner, P., Rogers, T.J., Lord, C., and Barthorpe, R.J., Learning model discrepancy: A Gaussian process and sampling-based approach, *Mechanical Systems and Signal Processing*, 152:107381, May 2021.

26

27

20. Plumlee, M., Bayesian Calibration of Inexact Computer Models, *Journal of the American Statistical Association*, 112(519):1274-1285, July 2017.

28

29

21. Gu, M. and Wang, L., Scaled Gaussian Stochastic Process for Computer Model Calibration and Prediction, *SIAM/ASA Journal on Uncertainty Quantification*, 6(4):1555-1583, January 2018.

30

31

22. Xie, F. and Xu, Y., Bayesian Projected Calibration of Computer Models, arXiv:1803.01231 [math, stat], February 2019.
1

2

3

Arhonditsis, G.B., Papantou, D., Zhang, W., Perhar, G., Massos, E., and Shi, M., Bayesian calibration of mechanistic aquatic biogeochemical models and benefits for environmental management, *Journal of Marine Systems*, 73(1-2):8–30, September 2008.

4

5

6

DeCarlo, E.C., Mahadevan, S., and Smarslok, B.P., Bayesian Calibration of Aerothermal Models for Hypersonic Air Vehicles, In *54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference*, Boston, Massachusetts, April 2013. American Institute of Aeronautics and Astronautics.

7

8

He, Y. and Xiu, D., Numerical strategy for model correction using physical constraints, *Journal of Computational Physics*, 313:617–634, May 2016.

9

10

Sargsyan, K., Huan, X., and Najm, H.N., Embedded Model Error Representation for Bayesian Model Calibration, *International Journal for Uncertainty Quantification*, 9(4), 2019.

11

12

Tuo, R. and Wu, C.F.J., A theoretical framework for calibration in computer models: Parametrization, estimation and convergence properties, *arXiv:1508.07155 [stat]*, August 2015.

13

14

Tuo, R. and Wu, C.F.J., Prediction based on the Kennedy-O'Hagan calibration model: Asymptotic consistency and other properties, *arXiv:1703.01326 [math, stat]*, March 2017.

15

16

Liu, F., Bayarri, M.J., and Berger, J.O., Modularization in Bayesian analysis, with emphasis on analysis of computer models, *Bayesian Analysis*, 4(1):119–150, March 2009.

17

18

Cox, D.D., Park, J.S., and Singer, C.E., A statistical method for tuning a computer code to a data base, *Computational Statistics & Data Analysis*, 37(1):77–92, July 2001.

19

20

Carmassi, M., Barbillon, P., Keller, M., Parent, E., and Chiodetti, M., Bayesian calibration of a numerical code for prediction, *arXiv:1801.01810 [stat]*, March 2019.

21

Stein, M.L., *Interpolation of Spatial Data*, Springer Series in Statistics, Springer New York, New York, NY, 1999.

22

23

Rasmussen, C.E. and Williams, C.K.I., *Gaussian Processes for Machine Learning*, Adaptive Computation and Machine Learning, MIT Press, Cambridge, Mass, 2006.

24

25

Robert, C.P. and Casella, G., *Monte Carlo Statistical Methods*, Springer Texts in Statistics, Springer, New York, NY, 2. ed., softcover reprint of the hardcover 2. ed. 2004 edition, 2010.

26

27

Yuan, J. and Ng, S.H., A sequential approach for stochastic computer model calibration and prediction, *Reliability Engineering & System Safety*, 111:273–286, March 2013.

28

29

Yuan, J., Nian, V., Su, B., and Meng, Q., A simultaneous calibration and parameter ranking method for building energy models, *Applied Energy*, 206:657–666, November 2017.

30

31

Yuan, J. and Ng, S.H., An Integrated Method for Simultaneous Calibration and Parameter Selection in Computer Models, *ACM Transactions on Modeling and Computer Simulation*, 30(1):7:1–7:23, February 2020.
1 38. Gramacy, R.B., Bingham, D., Holloway, J.P., Grosskopf, M.J., Kuranz, C.C., Rutter, E., Trantham, M., and Drake, R.P., Calibrating a large computer experiment simulating radiative shock hydrodynamics, *The Annals of Applied Statistics*, 9(3):1141–1168, September 2015.

4 39. Wu, X., Kozlowski, T., Meidani, H., and Shirvan, K., Inverse uncertainty quantification using the modular Bayesian approach based on Gaussian Process, Part 2: Application to TRACE, *Nuclear Engineering and Design*, 335:417–431, August 2018.

6 40. Joseph, V.R. and Melkote, S.N., Statistical Adjustments to Engineering Models, *Journal of Quality Technology*, 41(4):362–375, October 2009.

8 41. Wong, R.K.W., Storlie, C.B., and Lee, T.C.M., A frequentist approach to computer model calibration, *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 79(2):635–648, 2017.

10 42. Rumsey, K., Huerta, G., Brown, J., and Hund, L., Dealing with Measurement Uncertainties as Nuisance Parameters in Bayesian Model Calibration, *SIAM/ASA Journal on Uncertainty Quantification*, pp. 1287–1309, January 2020.

12 43. Loeppky, J.L., Bingham, D., and Welch, W.J., Computer Model Calibration or Tuning in Practice, Tech. Rep., 2006.

13 44. Arendt, P.D., Apley, D.W., and Chen, W., Quantification of Model Uncertainty: Calibration, Model Discrepancy, and Identifiability, *Journal of Mechanical Design*, 134(10):100908, October 2012.

15 45. Kennedy, M.C. and O'Hagan, A. Supplementary details on Bayesian Calibration of Computer Models. 2001.

16 46. Berge, C., *Topological Spaces: Including a Treatment of Multi-Valued Functions, Vector Spaces, and Convexity*, Macmillan Co., New York, 1963.

18 47. Tian, G. and Zhou, J., The maximum Theorem and the existence of Nash equilibrium of (generalized) games without lower semicontinuities, *Journal of Mathematical Analysis and Applications*, 166(2):351–364, May 1992.

20 48. Andrieu, C. and Thoms, J., A tutorial on adaptive MCMC, *Statistics and Computing*, 18(4):343–373, December 2008.

21 49. Gamerman, D. and Lopes, H.F., *Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference*, number 68 in Texts in Statistical Science Series, Taylor & Francis, Boca Raton, 2nd ed edition, 2006.

23 50. Kommajosyula, R., Development and assessment of a physics-based model for subcooled flow boiling with application to CFD, PhD thesis, Massachusetts Institute of Technology, 2020.

25 51. Kennel, W.E., Local boiling of water and superheating of high pressure steam in annuli, Thesis, Massachusetts Institute of Technology, 1949.

27 52. Leoni, N., Bayesian inference of model error for the calibration of two-phase CFD codes, PhD thesis, Institut Polytechnique de Paris, 2022.

29 53. Mchutchon, A. and Rasmussen, C., Gaussian Process Training with Input Noise, In *Advances in Neural Information Processing Systems*, Vol. 24. Curran Associates, Inc., 2011.