# Efficient treatment of the model error in the calibration of computer codes: the Complete Maximum a Posteriori method

Omar Kahol, Pietro Marco Congedo, Olivier Le Maitre, Enora Denimal Goy

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

## To cite this version:

Omar Kahol, Pietro Marco Congedo, Olivier Le Maitre, Enora Denimal Goy. Efficient treatment of the model error in the calibration of computer codes: the Complete Maximum a Posteriori method. International Journal for Uncertainty Quantification, 2025, 10.1615/Int.J.UncertaintyQuantification.2025056317. [hal-05090880](https://hal.science/hal-05090880v1)

HAL Id: hal-05090880

https://hal.science/hal-05090880v1

Submitted on 30 May 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Efficient treatment of the model error in the calibration
of computer codes: the Complete Maximum a Posteriori
method

Omar Kahol<sup>a,b</sup>, Pietro Marco Congedo<sup>a</sup>, Olivier Le Maître<sup>c</sup>, Enora Denimal
Goy<sup>a</sup>

<sup>a</sup>Inria, Centre de Mathématiques Appliquées, Ecole polytechnique, IPP, Route de Saclay,
91128 Palaiseau Cedex, France

<sup>b</sup>Dipartimento di Scienze e Tecnologie Aerospaziali, Politecnico di Milano, Via La Masa
34, 20156 Milano, Italy

<sup>c</sup>CNRS, Inria, Centre de Mathématiques Appliquées, Ecole polytechnique, IPP, Route de Saclay, 91128 Palaiseau Cedex, France

Abstract

Computer models are widely used for the prediction of complex physical phe-
nomena. Based on observations of these physical phenomena, it is possible
to calibrate the model parameters. In most cases, such computer models are
mis-specified, and the calibration process must be improved by including a
model error term. The model error hyperparameters are, however, rarely
learned jointly with the model parameters to reduce the dimensionality of
the problem. Sequential and non-sequential approaches have been intro-
duced to estimate the hyperparameters. The former, such as the Kennedy
and O'Hagan (KOH) framework, estimates the model error hyperparameters
before calibrating the model parameters. The latter, such as the Full Max-
imum a Posteriori (FMP), introduces a functional dependence between the
model parameters and the model error hyperparameters. Despite being more
reliable in some cases (bimodality e.g.), the FMP method still fails to esti-
mate correctly the posterior distribution shape. This work proposes a new
methodology for treating the model error term in computer code calibration.
It builds upon the KOH and FMP framework. Called the Complete Maxi-
mum a Posteriori (CMP) method, it provides a closed-form expression for the
marginalization integral over the model error hyperparameters, significantly

Email address: omar.kahol@inria.fr (Omar Kahol)

Preprint submitted to IJUQ

April 5, 2025
reducing the dimensionality of the calibration problem. Such expression relies on a set of assumptions that are more general and less stringent than the ones usually employed. The CMP method is applied to four examples of increasing complexity, from elementary to real fluid dynamics problems, including or not bimodality. Compared to the true reference solution and unlike the KOH and FMP, the CMP method correctly captures the shape of the posterior distribution, including all modes and their weights. Moreover, it provides an accurate estimate of the distribution tails.

Keywords: Uncertainty Quantification, Model Calibration, Bayesian
Method, Model Error

# 1. Introduction

A model is a tool used by scientists and engineers to understand and predict the world by mathematically describing different phenomena. The intrinsic mathematical details can be ignored here; however, it is useful to point out that a prediction generally depends on some parameters, $\theta$. These parameters can be observable quantities, such as dimensional constants, or non-observable variables which influence the prediction. The process of learning the value of such parameters, from available experimental data, is called *model calibration* and belongs to the class of inverse problems [1].

In this paper, we adopt the Bayesian perspective to solve the calibration problem. According to this perspective, the parameters, $\theta$, are random variables and the calibration process is done by updating their probability distribution using available experimental data. The *prior distribution* is the probability distribution of the parameters before the experimental data is observed, and the *posterior distribution* is the updated prior distribution after the experimental data is observed. The two are related by Bayes' formula; see [1–3] or §2 for more details.

The mathematical formulation of the model to be calibrated might present some internal inconsistencies or rely on simplifying assumptions; moreover, the numerical tools and approximations used to compute the model prediction introduce errors in the output. For these reasons, it is reasonable to assume that in most cases the model is an approximate representation of the true physical process. If this is the case, the solution of the calibration problem might yield biased and over-confident posteriors distributions [4, 5].
25

The seminal work of Kennedy and O'Hagan [6] was the first to consider an additional model error term in the Bayesian calibration procedure. The model error term is a random function that accounts for the discrepancy between the model and the true physical process and could be expressed in a finite-dimensional form (such as polynomial expansions) or with an infinite-dimensional form (such as Gaussian Processes) [1]. The model error term usually requires an additional parametrization; we call its parameters *hyperparameters* and denote them with $\psi$ to distinguish them from the model parameters, $\theta$.

30

35

The model error term has been applied to the calibration of complex models in many fields, such as aerodynamics [7], solid mechanics [8], energy [9] and climate [10, 11].

40

In a full Bayesian approach, the hyperparameters, $\psi$, are random variables and must be learned jointly with the parameters, $\theta$. The introduction of the model error term therefore increases the dimensionality and complexity of the calibration problem, calling for simplified solutions [12, 13]. These solutions usually target sampling from a lower-dimensional distribution.

45

First, sequential techniques propose to find a deterministic estimate of the hyperparameters, $\psi$, and then condition the joint posterior distribution on that deterministic estimate. The original calibration framework proposed by Kennedy and O'Hagan [6, 12] uses a sequential approach that is called the Kennedy O'Hagan (KOH) method.

50

Non-sequential methods, like the Full Maximum a Posteriori (FMP) method [13], approximate the marginal distribution of the parameters by introducing a functional relationship between the hyperparameters, $\psi$ and the parameters, $\theta$, and conditioning the joint posterior distribution on that relationship. The marginal posterior of the parameters is the integral of the joint distribution over all possible values of the hyperparameters, $\psi$. For this reason, approximating the marginal posterior is beneficial and can potentially capture more features of the posterior distribution. In [13], the authors show that in some cases, for example, bimodal distribution, the FMP method outperforms the KOH method as the latter sometimes fails to capture all the modes. In the same work it is pointed out that, although able to correctly capture all possible modes, the FMP method still fails to correctly represent their weights, incorrectly biasing the calibration results. These failures are mostly caused by the very stringent assumptions on which the FMP method relies. In particular, the FMP method assumes that when conditioning the posterior distribution of the hyperparameters on the parameters, the result is a point mass distribution. This assumption is very strict and, very rarely it is satisfied. The FMP method, therefore, misrepresents the shape of the

55

60
hyperparameters' conditional distribution. Moreover, the FMP method does not work when the prior of the hyperparameters is non-uniform (informative prior) as it is not present in its final expression.

In this paper, we revisit the FMP method and propose several improvements. In particular, we propose a new non-sequential method, the Complete Maximum a Posteriori (CMP) method, which provides a closed-form expression for the marginal distribution of the model parameters that is exact in some cases. The CMP method builds upon the FMP method and improves it by relaxing some assumptions and by providing a more accurate approximation of the posterior distribution. The name CMP refers to the fact that the new method incorporates a complete representation of the uncertainty in the model error term, from different values of the hyperparameters $\psi$ to potentially different shapes of their distribution.

The paper is organized in the following way. Section 2 presents a self-contained description of the calibration problem using the Bayesian framework. Section 3 is dedicated to the presentation of the CMP method. In §4, we test the CMP method on two elementary examples which, because of their simplicity, have a closed-form solution. In §5, we test the CMP method on two calibration problems and compare it to the KOH and FMP methods. Finally, §6 concludes the discussion and outlines possible future investigations.

## 2. The Bayesian Calibration Framework

This section briefly introduces the calibration problem with model error. The basic mathematical formulation and the notation are presented in §2.1. Section 2.2 focuses on the solution of the calibration problem discussing the different methods that can be used to solve the problem. In §2.3 we discuss the limitations of the current methods and outline the motivations for introducing the CMP method.

### 2.1. General Framework

We consider a mathematical model that describes a physical system. The response is an abstract quantity $u \in \mathcal{U}$ that represents the system's state given a value of some measurable control variables, $x \in \mathcal{X}$. The mathematical model, $\mathcal{M}$, is represented as an operator that operates on the response and is parameterized by some parameters, $\theta \in \Theta \subset \mathbb{R}^{d_{\theta}}$. If $u$ is the solution of the model, the following equation holds:

$$ \mathcal{M}(u | x, \theta) = 0. \qquad (1) $$
The solution, $\mathbf{u}$, is usually not directly measurable as it is an abstract representation of the response of the physical system. One usually has access to some observables, $\{y^{(i)}\}_{i=1...N_h}$: $y^{(i)} \in \mathcal{Y}^{(i)}$, which are scalar physical quantities that can be measured in an experiment. To simplify the discussion, we initially deal with the case in which only a single scalar observable, $y$, is available. The extension to the case of multiple observables is straightforward and will be discussed later. We model the observable, $y$, as the sum between an observation operator, $\mathcal{G}$, which operates on the solution of the model, $\mathbf{u}$, and a residual term, $r$,

$$y | \mathbf{x} = \mathcal{G}(\mathbf{u} | \mathbf{x}, \boldsymbol{\theta}) + r(\mathbf{x} | \boldsymbol{\theta}). \qquad (2)$$

The presence of the residual term, $r$, is due to the fact that the model is an approximate representation of the true physical system. The residual term is usually unknown but depends on the model parameters, $\boldsymbol{\theta}$, since a different choice of $\boldsymbol{\theta}$ will lead to a different residual term [13]. We call a model *well-specified* if there exists a value of the parameters, $\boldsymbol{\theta}^*$, that makes the residual term, $r$, zero: $r(\mathbf{x} | \boldsymbol{\theta}^*) = 0$. We call a model *miss-specified* if otherwise. In this work, we consider the case of miss-specified models.

The term *model calibration* refers to the process of learning the value of the parameters $\boldsymbol{\theta}$ given some experimental data, $\mathcal{D} = \{(\mathbf{x}_{\text{obs}}^i, y_{\text{obs}}^i)\}_{i=1...N}$, where $\mathbf{x}_{\text{obs}}^i$ represents the value of the control variables and $y_{\text{obs}}^i$ is the experimental measurement of the observable, $y$. In this work, we consider the case in which the experimental observations, $y_{\text{obs}}$, are affected by measurement error. To proceed, we need to define a statistical model (and some assumptions) that explains the observations. We model the observations as a vector of $N$ random variables, $\mathbf{Y}$, which is the sum of the outputs of three random functions:

$$\mathbf{Y} = \mathcal{G}_{\boldsymbol{\theta}} + \mathbf{z} + \boldsymbol{\epsilon}. \qquad (3)$$

The random vector $\mathcal{G}_{\boldsymbol{\theta}}$ contains the output of the observation operator at each control variable, $\mathcal{G}_{\boldsymbol{\theta}} = (\mathcal{G}(\mathbf{u} | \mathbf{x}_{\text{obs}}^i, \boldsymbol{\theta}))_{i=1...N}$. Following the Bayesian framework, $\boldsymbol{\theta}$ is a random variable and it is distributed according to the prior distribution, $\pi(\boldsymbol{\theta})$.

The random vector $\mathbf{z} = (z(\mathbf{x}_{\text{obs}}^i | \psi_z))_{i=1...N}$ is the model of the residual, $r(\mathbf{x} | \boldsymbol{\theta})$, and is called *model error* [6]. In this work, the model error term is considered to be a Gaussian Process (GP) with mean $\mu_{\psi_z}(\mathbf{x})$ and covariance $c_{\psi_z}(\mathbf{x}, \mathbf{x}')$. The GP is assumed to have zero mean, leading to better identifiability [14]. The parameters $\psi_z$ are the hyperparameters that parameterize the model error term. Similarly to the model parameters, the
hyperparameters are random variables and are distributed according to the prior distribution $\pi(\psi_z)$.

The random vector $\epsilon$ contains the measurement error at each control point, $\epsilon = (\epsilon_i)_{i=1...N}$. In this work, we assume that the $\epsilon_i$ are i.i.d. centered Gaussian random variables with equal variance $\sigma_e^2$ and call $\pi(\sigma_e)$ its prior distribution. The hyperparameters $\psi = (\psi_z, \sigma_e) \in \Psi \subset \mathbb{R}^{d_{\psi}}$ are the vector of all the hyperparameters that parameterize the model and experimental errors.

The prior distributions are updated using the available experimental data to obtain the posterior distribution, $p(\theta, \psi | \mathcal{D})$. The posterior distribution is related to the prior distributions by Bayes' formula [1]:

$$p(\theta, \psi | \mathcal{D}) = \frac{\mathcal{L}(\mathcal{D} | \theta, \psi) \pi(\theta) \pi(\psi)}{p(\mathcal{D})} . \quad (4)$$

The term $\mathcal{L}(\mathcal{D} | \theta, \psi)$ is the likelihood that measures the probability of observing the experimental data given a particular value of the parameters and hyperparameters. Following our assumptions on Eq. 3, the likelihood is a multivariate normal distribution:

$$\mathcal{L}(\mathcal{D} | \theta, \psi) = \frac{1}{\sqrt{(2\pi)^N \det(\mathbf{K}_{\psi_z} + \sigma_e^2 \mathbf{I}_N)}} \exp\left(-\frac{1}{2} \mathbf{r}_{\theta}^T (\mathbf{K}_{\psi_z} + \sigma_e^2 \mathbf{I}_N)^{-1} \mathbf{r}_{\theta}\right). \quad (5)$$

The vector $\mathbf{r}_{\theta}$ contains the residuals, $\mathbf{r}_{\theta} = (y_{\text{obs}}^i - \mathcal{G}(\mathbf{u} | \theta, \mathbf{x}_{\text{obs}}^i))_{i=1...N}$, and $\mathbf{K}_{\psi_z}$ is the covariance matrix of the model error term evaluated at the observation points $(\mathbf{K}_{\psi_z})_{ij} = c_{\psi_z}(\mathbf{x}_{\text{obs}}^i, \mathbf{x}_{\text{obs}}^j)$. The notation $\det(\cdot)$ denotes the determinant of a matrix and $\mathbf{I}_N$ is the identity matrix of size $N$. In the literature, one can find different expressions of Eq. 5 that can account also for data inconsistency [15], multiplicative experimental error [16], and model errors with non-zero mean [1].

When multiple observables are available, the likelihood is the product of the likelihoods of each observable,

$$\mathcal{L}(\mathcal{D}^{(1)}, \dots, \mathcal{D}^{(N_h)} | \theta, \psi^{(1)}, \dots, \psi^{(N_h)}) = \prod_{i=1}^{N_h} \mathcal{L}(\mathcal{D}^{(i)} | \theta, \psi^{(i)}), \quad (6)$$

where $\psi^{(i)}$ is the hyperparameters vector that parameterize the model and experimental error terms for the $i$-th observable and $\mathcal{D}^{(i)}$ is the experimental data concerning the $i$-th observable.

Finally, the denominator in Eq. 4 is the *model evidence* given by:
$$ p(\mathcal{D}) = \int_{\Theta} \int_{\Psi} \mathcal{L}(\mathcal{D} | \theta, \psi) \pi(\theta) \pi(\psi) \, d\psi \, d\theta. \quad (7) $$

Integrating the joint posterior in Eq. 4 over the hyperparameters yields the parameters' marginal posterior,

$$ p(\theta | \mathcal{D}) = \int_{\Psi} p(\theta, \psi | \mathcal{D}) \, d\psi. \quad (8) $$

The corrected model prediction of an observable at a new coordinate, $x^*$, denoted with $y_{corr}^* | x^* = \mathcal{G}(u | x^*, \theta) + z(x^* | \psi_z)$, is the sum of the model and model error prediction. Its marginal distribution can be computed using

$$ p(y_{corr}^* | \mathcal{D}) = \int_{\Theta} \int_{\Psi_z} p(y_{corr}^* | \theta, \psi_z, \mathcal{D}) p(\theta, \psi_z | \mathcal{D}) \, d\theta \, d\psi_z, \quad (9) $$

where the distribution $p(y_{corr}^* | \theta, \psi_z, \mathcal{D})$ is normal with mean and variance specified by the Gaussian process predictive equations [17]:

$$ \begin{aligned} p(y_{corr}^* | \theta, \psi_z, \mathcal{D}) = & \mathcal{N}(y_{corr}^* | \mu_{pred}, \sigma_{pred}^2), \\ \mu_{pred} = & \mathcal{G}(u | x^*, \theta) + k_*^T K_{\psi_z}^{-1} r_{\theta}, \\ \sigma_{pred}^2 = & c_{\psi_z}(x^*, x^*) - k_*^T K_{\psi_z}^{-1} k_* \end{aligned} \quad (10) $$

where $(k_*)_i = c_{\psi_z}(x^*, x_{\text{obs}}^i)$.

This framework, together with the above-stated assumptions, constitutes the Bayesian Calibration Framework and provides the tools and methodology to calibrate mis-specified models. The expression of the posterior distribution, Eq. 4, is the starting point for the solution of the calibration problem but one rarely has access to a closed-form expression.

## 2.2. Solving the Bayesian Calibration Problem

A popular way to solve the calibration problem is to use sampling techniques, like Markov Chain Monte Carlo (MCMC) methods [1, 2, 18]. MCMC methods construct a Markov Chain whose steady-state distribution is the target distribution and extract samples when convergence is reached [19].

To extract samples from the parameters' posterior distribution, Eq. 8, one can use the so-called full Bayesian approach. In this case, one samples from the joint posterior distribution, $p(\theta, \psi | \mathcal{D})$, and then marginalizes over the hyperparameters to obtain the parameters' posterior distribution, $p(\theta | \mathcal{D})$. Examples of this technique can be found in [13, 14, 20]. The full Bayesian
approach is the most general but it is also the most computationally expensive as its efficiency decreases as the dimensionality of the problem increases [18].

Modular approaches, on the other hand, reduce the complexity of the problem by sampling from a lower dimensional distribution. Sequential approximations fix the hyperparameters to a constant value, $\bar{\psi}$, and sample from the parameters' posterior distribution conditioned on that value,

$$ p(\theta | \mathcal{D}, \bar{\psi}) \propto \mathcal{L}(\mathcal{D} | \theta, \bar{\psi}) \pi(\theta). \qquad (11) $$

Different authors argue for different choices for the deterministic estimate of the hyperparameters, $\bar{\psi}$, and discuss the implications of such choices [21–24]. Of these, we report the KOH method [6, 12], which fixes the hyperparameters to the maximizer of the parameter-averaged posterior distribution,

$$ \begin{aligned} \bar{\psi} = \psi_{\text{KOH}} &= \arg\max_{\psi \in \Psi} \int_{\Theta} p(\theta, \psi | \mathcal{D}) \, d\theta, \\ &= \arg\max_{\psi \in \Psi} p(\psi | \mathcal{D}). \end{aligned} \qquad (12) $$

Sequential methods are potentially cheaper than the full Bayesian approach but they might miss important features of the posterior distribution. The FMP method [13] is an example of a non-sequential method that approximates the marginal distribution of the parameters by introducing a functional relationship between the hyperparameters and the parameters,

$$ \psi_{\text{MAP}}(\theta) = \arg\max_{\psi \in \Psi} p(\psi, \theta | \mathcal{D}). \qquad (13) $$

The method then proceeds by conditioning the likelihood on the functional relationship, $\psi_{\text{MAP}}(\theta)$, and sampling from the following distribution:

$$ P_{\text{FMP}}(\theta | \mathcal{D}) \propto \mathcal{L}(\mathcal{D} | \theta, \psi_{\text{MAP}}(\theta)) \pi(\theta). \qquad (14) $$

The authors of the FMP method claim that the method is more accurate than the KOH method and that it approximates the marginal distribution of the parameters Eq. 8.

Note that $\psi_{\text{KOH}}$ and $\psi_{\text{MAP}}(\theta)$ can be related. This is done, initially, by using the marginalization formula,

$$ p(\theta, \psi | \mathcal{D}) = p(\psi | \theta, \mathcal{D}) p(\theta | \mathcal{D}). \qquad (15) $$

Substituting Eq. 15 in Eq. 12 yields
$$
\begin{align}
\psi_{\text{KOH}} &= \arg\max_{\psi \in \Psi} \mathbb{E}_{\theta \sim \text{p}(\theta|\mathcal{D})} [\text{p}(\psi | \theta, \mathcal{D})] \nonumber \\
&= \mathbb{E}_{\theta \sim \text{p}(\theta|\mathcal{D})} \left[ \arg\max_{\psi \in \Psi} \text{p}(\psi | \theta, \mathcal{D}) \right] \tag{16} \\
&= \mathbb{E}_{\theta \sim \text{p}(\theta|\mathcal{D})} [\psi_{\text{MAP}}(\theta)] . \nonumber
\end{align}
$$

The KOH value of the hyperparameters is the posterior-averaged value of
the functional relationship between the hyperparameters and the parameters
introduced by the FMP method.

2.3. Limitations

In the previous subsection, several approximations of the full Bayesian cal-
ibration problem were introduced. On the one hand, these approaches are
different both from the methodological point of view and from the theoretical
point of view. On the other hand, they share some interesting similarities as
their sampling distribution is a conditioning of the true posterior distribution
on a curve, $\hat{\psi}(\theta)$. In the case of sequential methods, this curve is a straight
line, $\hat{\psi}(\theta) = \bar{\psi}$, while in the case of the FMP method, the curve is more com-
plex, $\hat{\psi}(\theta) = \psi_{\text{MAP}}(\theta)$. The FMP method is therefore superior to sequential
methods as it can potentially explore more regions of the posterior.

The exact representation of the marginal posterior of the parameters is,
however, the result of the integration of Eq. 8. The FMP method claims
to approximate the parameters’ marginal distribution, but there are two
main problems. First and foremost, the FMP method assumes that the
hyperparameters conditional distribution is a point mass p(ψ|θ, D, M) ≈
δ(ψ − ψMAP(θ)). This assumption turns out to be very strict and does not
consider the true shape of the hyperparameters’ conditional distribution. The
second problem is that the FMP parameters’ posterior distribution, Eq. 14,
evaluates only the likelihood on ψMAP(θ) and not the prior of the hyper-
parameters. This can cause problems when the prior distribution of the
hyperparameters is not uniform.

To lift these limitations, we propose in this work a new non-sequential ap-
proach, the CMP method, which builds upon the FMP method and im-
proves it by relaxing its assumptions. Under these new assumptions, the
CMP method is an exact representation of the marginal posterior distribu-
tion of the parameters, increasing the efficiency and reliability of the model
calibration.
### 3. The Complete Maximum a Posteriori Method

This section introduces the CMP method. Section 3.1 formulates the assumptions on which the method relies. Sections 3.2 and 3.3 lay down the CMP approximation of the marginal and predictive posterior distributions showing that, if the assumptions are satisfied, the resulting distribution is equal to the full Bayesian one. Section 3.4 and Appendix A discuss some computational aspects.

#### 3.1. Assumptions of the CMP Method

We assume that the hyperparameters' conditional distribution, $p(\psi | \theta, \mathcal{D})$, has a unique global maximum value, $\psi_{MAP}(\theta) \in \Psi$, for each possible value of the parameters, $\theta$. This maximum is defined in the FMP approximation by Eq. 13 and [13].

Furthermore, we assume that the conditional distribution $p(\psi | \theta, \mathcal{D})$ can be approximated by a fixed distribution, $f: \Phi \to \mathbb{R}^{+}$, which we call shape function. The conditional distribution and the shape function are connected via an affine transformation

$$ \phi = S_{\theta}(\psi - \psi_{MAP}(\theta)). \qquad (17) $$

The transformation matrix $S_{\theta}$ must be invertible. Therefore,

$$ p(\psi | \theta, \mathcal{D}) \approx f(\phi) | \det(S_{\theta}) | . \qquad (18) $$

These hypotheses correspond to requiring that the hyperparameters' conditional distribution, as $\theta$ changes, changes its shape only by some scaling and shifting.

<figure><img src="image_3.png" /><figcaption>Figure 1: Example (sketch) of a consistent hyperparameters' conditional distribution for two values of the parameters (blue and orange curves).</figcaption></figure>
In Fig. 1 we report a sketch of a consistent hyperparameters’ conditional distribution for two values of the parameters (blue and orange curves). Here the shape function is $f(\phi) \propto \text{sinc}^2(\phi)$, and we have chosen this example to stress the fact that, although the function can have a multiple local maxima, the global maximum should be unique. As visible in the figure, the affine transformation correspond to some shifting and scaling of the shape function. We remind that if the shape function is a point mass distribution, $f(\phi) = \delta(\phi)$, the CMP method reduces to the FMP method. Moreover, if $f(\phi) = \frac{1}{\sqrt{2\pi}} \exp -\frac{\phi^T \phi}{2}$, we recover the Laplace approximation.

In any case, computing the shape function is not necessary for the computation of the CMP approximation of the marginal posterior distribution, and we only need to postulate its existence. We anticipate, however, that for the computation of the predictive posterior distribution, the shape function must be known, and we are going to present a further approximation to avoid computing it in §3.3.

## 3.2. Computation of the Marginal Posterior

The parameters’ posterior distribution must be computed by integrating out the hyperparameters as in Eq. 8. To simplify the integral we note that, under the CMP assumptions, we have that

$$ p(\psi_{\text{MAP}}(\theta) | \theta, \mathcal{D}) \frac{1}{|\det(\mathcal{S}_{\theta})|} \approx f(0). \qquad (19) $$

We can therefore multiply the left and right-hand sides of Eq. 19 by the parameters’ marginal posterior distribution, $p(\theta | \mathcal{D}, \mathcal{M})$, and obtain

$$ \begin{aligned} p(\theta | \mathcal{D}) f(0) &\approx p(\theta | \mathcal{D}) p(\psi_{\text{MAP}}(\theta) | \theta, \mathcal{D}) \frac{1}{|\det(\mathcal{S}_{\theta})|} \qquad (20) \\ &= p(\theta, \psi_{\text{MAP}}(\theta) | \mathcal{D}) \frac{1}{|\det(\mathcal{S}_{\theta})|}. \end{aligned} $$

This means that marginalizing the joint distribution over the hyperparameters, a procedure that previously required the evaluation of a possibly high-dimensional integral, can now be done by evaluating the joint distribution along a curve, $c(\theta) = (\theta, \psi_{\text{MAP}}(\theta))$, and dividing by a correction factor, $|\det(\mathcal{S}_{\theta})|$. Inserting the expression of the joint posterior distribution, Eq. 4, in Eq. 20 we get the expression of the CMP marginal posterior distribution,

$$ p_{\text{CMP}}(\theta | \mathcal{D}) \propto \mathcal{L}(\mathcal{D} | \theta, \psi_{\text{MAP}}(\theta)) \pi(\theta) \pi(\psi_{\text{MAP}}(\theta)) \frac{1}{|\det(\mathcal{S}_{\theta})|}. \qquad (21) $$
The expression still relies on the availability of $|\det(S_\theta)|$ whose expression is discussed in §3.4. One can immediately see that, contrary to the FMP method, the CMP method evaluates both the likelihood and the prior distribution of the hyperparameters on $\psi_{\text{MAP}}(\theta)$. This is a crucial difference as the prior distribution of the hyperparameters can be informative and can influence the parameters' marginal posterior distribution. Furthermore, dividing by the correction factor, $|\det(S_\theta)|$, ensures that the CMP distribution correctly recovers the parameters' marginal distribution, Eq. 8, if the assumptions of the method are satisfied.

290

295

### 3.3. Approximation of the Predictive Posterior

We now turn to the question of approximating the posterior distribution of the corrected model, $p(y_{\text{corr}}^* | \mathcal{D})$ defined in Eq. 9, i.e. the distribution of the observable, $y_{\text{corr}}^* = \mathcal{G}(u | \theta, x^*) + z(x^* | \psi_z)$, at a new coordinate, $x^*$. To simplify the notation we remove the explicit dependence on the data,

$$p(y_{\text{corr}}^*) = \int_{\Theta} \int_{\Psi} p(y_{\text{corr}}^* | \theta, \psi) p(\theta, \psi) \, d\psi \, d\theta$$

Using Eq. 15,

$$= \int_{\Theta} p(\theta) \int_{\Psi} p(y_{\text{corr}}^* | \theta, \psi) p(\psi | \theta) \, d\psi \, d\theta$$

From Eq. 18,

(22)

$$\propto \int_{\Theta} p(\theta) \int_{\Psi} p(y_{\text{corr}}^* | \theta, \psi) f(\phi) |\det(S_{\theta})| \, d\psi \, d\theta$$

Using Eq. 17

$$= \int_{\Theta} p(\theta) \, d\theta \int_{\Phi} p(y_{\text{corr}}^* | \theta, \psi_{\text{MAP}}(\theta) + S_{\theta}^{-1}\phi) f(\phi) \, d\phi.$$

300

We introduce the following abbreviation for the argument of the integral,

$$\begin{align}
P_{\text{CMP}} (y_{\text{corr}}^* | \theta) &= \int_{\Phi} p (y_{\text{corr}}^* | \theta, \psi_{\text{MAP}} (\theta) + S_{\theta}^{-1} \phi) f(\phi) \, d\phi \tag{23} \\
&= \mathbb{E}_{\phi \sim f(\phi)} [ p (y_{\text{corr}}^* | \theta, \psi_{\text{MAP}} (\theta) + S_{\theta}^{-1} \phi) ] . \nonumber
\end{align}$$

This integral cannot be evaluated without an explicit expression for the shape function, $f(\phi)$, preventing direct sampling from $p(y_{\text{corr}}^* | \mathcal{D})$. To address this issue, we propose to return to a point-wise approximation, by assuming that the shape function is a point mass distribution, $f(\phi) = \delta(\phi)$,

$$\hat{P}_{\text{CMP}} (y_{\text{corr}}^* | \theta, \mathcal{D}) = p (y_{\text{corr}}^* | \theta, \psi_{\text{MAP}} (\theta), \mathcal{D}) . \quad (24)$$
The proposed approximation, Eq. 24, can be evaluated using the GPs predictive equations [17], reported in §2.1. The result can then be averaged over the parameters’ marginal posterior distribution, Eq. 21, to obtain the predictive distribution of the corrected model,

$$ p(y_{\text{corr}}^* | \mathcal{D}) \approx \int_{\Theta} \hat{p}_{\text{CMP}}(y_{\text{corr}}^* | \boldsymbol{\theta}, \mathcal{D}) p_{\text{CMP}}(\boldsymbol{\theta} | \mathcal{D}) d\boldsymbol{\theta}. \quad (25) $$

On the one hand, introducing the point mass distribution assumption neglects the uncertainty of the hyperparameters and makes Eq. 25 overconfident. On the other hand, the uncertainty of the hyperparameters is already accounted in the parameters’ marginal posterior distribution which is used in Eq. 25 to average Eq. 24. This means that the CMP predictive distribution, although not exact, is still an improvement over the FMP predictive distribution.

### 3.4. Evaluation of the Correction Factor

The CMP method relies on the availability of the correction factor, $|\det(S_\theta)|$, which, as of now, was not specified. We can however show that this quantity is connected to the Hessian of the posterior distribution at its maximum value in the hyperparameters. This is firstly done by substituting the hypothesis of the CMP method, Eq. 18, in the log of Eq. 15:

$$ \nabla_{\psi}^{2} \log p(\boldsymbol{\theta}, \psi_{\text{MAP}}(\boldsymbol{\theta}) | \mathcal{D}, \mathcal{M}) \approx \nabla_{\psi}^{2} \log f(\mathbf{0}). \quad (26) $$

Using the chain rule we obtain

$$ \nabla_{\psi}^{2} \log f(\mathbf{0}) = S_{\theta}^{T} H_{0} S_{\theta}, \quad (27) $$

where $H_0 = \nabla_{\psi}^2 \log f(\mathbf{0})$ is a constant matrix equal to the Hessian of the log of the shape function at the origin. Because of the rules of the determinant,

$$ \sqrt{|\det(\nabla_{\psi}^2 \log p(\boldsymbol{\theta}, \psi_{\text{MAP}}(\boldsymbol{\theta}) | \mathcal{D}))|} \propto |\det(S_{\theta})|. \quad (28) $$

The Hessian of the posterior distribution is available in closed form provided that the derivatives of the model error’s covariance function and the log-prior are available. Appendix A reports a closed-form expression of the Hessian of the log-likelihood with respect to the hyperparameters for the case considered in this paper.
### 4. Elementary Examples

In this section, we present two elementary examples to illustrate the application of the CMP method and show its advantages over the KOH and FMP methods. Both examples are chosen because they can be solved analytically, giving valuable insights into the method. The first example, in §4.1, is a calibration problem in which the model is well-specified. The second one, in §4.2, is a calibration problem in which the true posterior distribution is a mixture of Gaussian distributions with well-separated modes.

#### 4.1. Example 1: Calibration of a Well-Specified Model

In this subsection, we start from Eq. 5 and simplify it by assuming that the model error term is zero, hence the model is well-specified. The only hyper-parameter left is the experimental error noise, $\sigma_e$. Under this assumption, Eq. 5 simplifies to,

$$ \mathcal{L}(\mathcal{D} | \theta, \sigma_e) = \frac{1}{\sqrt{(2\pi\sigma_e^2)^N}} \exp\left(-\frac{\mathbf{r}_{\theta}^T \mathbf{r}_{\theta}}{2\sigma_e^2}\right). \qquad (29) $$

Concerning the prior distribution of the experimental error noise, $\pi(\sigma_e)$, we choose an improper power law prior distribution with parameter $a > 0$,

$$ \pi(\sigma_e) = \frac{1}{\sigma_e^a}. \qquad (30) $$

If $a = 0$ we recover a uniform prior distribution, if $a = 1$ we recover a log-uniform prior distribution and if $a = 2$ we recover Jeffrey's prior. The solution is reported as a function of the parameter $a$. Using Eq. 4 we can compute the unnormalized posterior distribution. The hyperparameters' conditional distribution is, up to a normalization constant,

$$ \begin{aligned} p(\sigma_e | \theta, \mathcal{D}) &= \frac{p(\theta, \sigma_e | \mathcal{D})}{\int_0^\infty p(\theta, \sigma_e | \mathcal{D}) \, d\sigma_e} \\ &\propto \frac{1}{\sqrt{\mathbf{r}_\theta^T \mathbf{r}_\theta}} e^{-\frac{\mathbf{r}_\theta^T \mathbf{r}_\theta}{2\sigma_e^2}} \left( \frac{\sqrt{\mathbf{r}_\theta^T \mathbf{r}_\theta}}{\sigma_e} \right)^{N+a}. \end{aligned} \qquad (31) $$

We can compute $\sigma_e^{\text{MAP}}$ and $|\det(S_\theta)|$ by following Eq. 13 and Eq. 28 respectively. We obtain
$$ \sigma_e^{\text{MAP}}(\boldsymbol{\theta}) = \sqrt{\frac{\boldsymbol{r}_{\boldsymbol{\theta}}^T \boldsymbol{r}_{\boldsymbol{\theta}}}{N+a}}, \qquad (32) \\ |\det(S_{\boldsymbol{\theta}})| = \sqrt{2 \frac{a+N}{\boldsymbol{r}_{\boldsymbol{\theta}}^T \boldsymbol{r}_{\boldsymbol{\theta}}}} \propto \sqrt{\frac{1}{\boldsymbol{r}_{\boldsymbol{\theta}}^T \boldsymbol{r}_{\boldsymbol{\theta}}}}. $$

Calling $t = \frac{\sigma_e}{\sqrt{\boldsymbol{r}_{\boldsymbol{\theta}}^T \boldsymbol{r}_{\boldsymbol{\theta}}}}$ and by using Eq. 32, we can rewrite Eq. 31 as

$$ p(\sigma_e | \boldsymbol{\theta}, \mathcal{D}) \propto |\det(S_{\boldsymbol{\theta}})| e^{-\frac{1}{2t^2}} \left(\frac{1}{t}\right)^{N+a}. \qquad (33) $$

We can recover the affine transformation by noting that

$$ \phi = S_{\boldsymbol{\theta}}(\sigma_e - \sigma_e^{\text{MAP}}) = t - \sqrt{\frac{1}{N+a}} = t - \phi_0. \qquad (34) $$

Finally, substituting Eq. 34 in Eq. 33 it yields

$$ p(\sigma_e | \boldsymbol{\theta}, \mathcal{D}) \propto |\det(S_{\boldsymbol{\theta}})| e^{-\frac{1}{2(\phi + \phi_0)^2}} \left(\frac{1}{\phi + \phi_0}\right)^{N+a}. \qquad (35) $$

As it is possible to see, Eq. 35 respects exactly the hypothesis of the CMP method. It can be rewritten as described by Eq. 18 with

$$ f(\phi) : (-\phi_0, \infty) \rightarrow \mathbb{R}^+, \\ \propto e^{-\frac{1}{2(\phi + \phi_0)^2}} \left(\frac{1}{\phi + \phi_0}\right)^{N+a}. \qquad (36) $$

<figure><img src="image_4.png" /><figcaption>Figure 2: Shape function.</figcaption></figure>
Fig. 2 depicts the shape function for the CMP calibration problem with no model error choosing $N + a = 10$. In this case, the CMP approximation will yield the exact result.

The correction factor, given in Eq. 32, is crucial in this case. Omitting it would lead to a parameter's posterior distribution that overemphasizes areas where the residuals are small and underemphasizes areas where the residuals are large. Since large residuals are typically found in the tails of the distributions, the correction factor is essential to prevent the false certainty effect.

### 4.2. Example 2: Mixture of Gaussians

In this subsection, we assess the CMP approximation of the parameters' marginal posterior distribution in the case where the posterior distribution is a mixture of Gaussian distributions with well-separated modes. This example was already treated by [13], which reported the KOH and FMP approximations. They found that the KOH method can capture only a single mode and that the FMP approximation, while being able to capture all the modes, misrepresents their weights. We shall therefore perform similar computations to show that the CMP method can capture both the modes and their weights, giving rise to an exact representation of the posterior distribution.

We consider the case in which the joint posterior, $p(\gamma = (\theta, \psi) | \mathcal{D})$, is a mixture of $m$ Gaussian distributions with weights $w_i: \sum_{i=1}^{m} w_i = 1$, means $\mu_i = (\mu_{\theta,i} \mu_{\psi,i})$, and covariance matrices $K_i = [K_{\theta,i} \quad K_{\theta,\psi,i} \quad K_{\psi,i}]$. The posterior distribution is given by

$$ p(\theta, \psi | \mathcal{D}) \propto \sum_{i=1}^{m} w_i \frac{1}{\sqrt{(2\pi)^{d_{\theta}+d_{\psi}} \det(K_i)}} \exp \left( -\frac{1}{2}(\gamma - \mu_i)^T K_i^{-1} (\gamma - \mu_i) \right). \quad (37) $$

We assume that the modes are well-separated. This translates into requiring that “the intervals: $\Theta_i = \{\theta : (\mu_{i,\theta} - \theta)^T K_{\theta,i}^{-1} (\mu_{i,\theta} - \theta) \le t_{95\%}^{d_{\theta}}\}_{i=1...m}$ (with $t_{95\%}^{d_{\theta}}$ equal to the 95% quantile of the $\chi_2$ law with $d_{\theta}$ degrees of freedom) are disjoint” [13]. An equivalent condition is supposed to hold for the hyperparameters. Under this hypothesis, the hyperparameters' conditional distribution can be well approximated by
$$ \begin{aligned} \mathrm{p}(\psi | \theta, \mathcal{D}) &\approx \sum_{i=1}^{m} \mathcal{I}(\theta \in \Theta_i) \frac{1}{\sqrt{(2\pi)^{d_{\psi}} \det(\mathbf{K}_{\psi,i|\theta})}} \\ &\exp\left(-\frac{1}{2}(\psi - \boldsymbol{\mu}_{\psi,i|\theta})^T \mathbf{K}_{\psi,i|\theta}^{-1} (\psi - \boldsymbol{\mu}_{\psi,i|\theta})\right). \end{aligned} \quad (38) $$

The conditional means can be expressed as $\mu_{\psi,i|\theta} = \mu_{\psi,i} + K_{\theta,\psi,i}^T K_{\theta,i}^{-1}(\theta - \mu_{\theta,i})$ and the conditional covariances with $K_{\psi,i|\theta} = K_{\psi,i} - K_{\theta,\psi,i}^T K_{\theta,i}^{-1} K_{\theta,\psi,i}$. The function $\mathcal{I}(\theta \in \Theta_i)$ is the indicator function that is equal to 1 if $\theta \in \Theta_i$ and 0 otherwise. Note that Eq. 38 holds because we assumed that the modes are well-separated and hence we can treat each mode as a separated Gaussian distribution. Furthermore, conditioning a multivariate Gaussian distribution on a subset of its variables results in a Gaussian distribution [17].

According to [13] and Eq. 28, the optimal hyperparameters, $\psi_{\mathrm{MAP}}(\theta)$, and the correction factor, $|\det(\mathbf{S}_{\theta})|$, are well approximated by the following piecewise constant functions

$$ \begin{aligned} \psi_{\mathrm{MAP}}(\theta) &\approx \sum_{i=1}^{m} \mu_{\psi,i|\theta} \mathcal{I}(\theta \in \Theta_i), \\ |\det(\mathbf{S}_{\theta})| &\approx \sum_{i=1}^{m} \frac{1}{\sqrt{\det(\mathbf{K}_{\psi,i|\theta})}} \mathcal{I}(\theta \in \Theta_i). \end{aligned} \quad (39) $$

<figure><img src="image_5.png" /><figcaption>Figure 3: Example of a mixture of Gaussian distribution. The red line represents the MAP value of the hyperparameters, $\psi_{\mathrm{MAP}}(\theta)$, and the blue line represents the KOH value.</figcaption></figure>

Fig. 3 shows an example of a mixture of Gaussian distribution with two modes centered in $(-1, 0)$ and $(1, -1)$ respectively. The covariance matrices are diagonal with elements $(0.1, 0.01)$ and $(0.1, 0.1)$ respectively. The red line
represents the MAP value of the hyperparameters, $\psi_{\text{MAP}}(\theta)$, and the blue line represents the KOH value whose expression is reported in [13]. Using the definition of the CMP method, Eq. 21, with Eq. 39 and the hypothesis of well-separated modes, we can compute the CMP approximation of the parameters’ marginal posterior distribution,

$$ \hat{p}_{\text{CMP}}(\boldsymbol{\theta} | \mathcal{D}) \propto \sum_{i=1}^{m} w_i \frac{\sqrt{\det(\mathbf{K}_{\psi,i|\boldsymbol{\theta}})}}{\sqrt{(2\pi)^{d_{\theta}+d_{\psi}} \det(\mathbf{K}_i)}} \exp\left(-\frac{1}{2}\left(\begin{array}{cc} \boldsymbol{\theta}-\boldsymbol{\mu}_{\theta,i} & \boldsymbol{\mathcal{K}}_{i}^{-1} \\ \boldsymbol{\mu}_{\psi,i|\boldsymbol{\theta}}-\boldsymbol{\mu}_{\psi,i} & \boldsymbol{\theta}-\boldsymbol{\mu}_{\theta,i} \end{array}\right)^{\mathrm{T}} \mathbf{K}_{i}^{-1}\left(\begin{array}{cc} \boldsymbol{\theta}-\boldsymbol{\mu}_{\theta,i} & \boldsymbol{\mu}_{\psi,i|\boldsymbol{\theta}}-\boldsymbol{\mu}_{\psi,i} \end{array}\right)\right). \quad (40) $$

We can exploit the block structure of the covariance matrix, $\mathbf{K}_i$, to express its determinant as $\det(\mathbf{K}_i) = \det(\mathbf{K}_{\theta,i}) \det(\mathbf{K}_{\psi,i|\theta})$ [25]. Also, the argument of the exponential simplifies to $-\frac{1}{2}(\boldsymbol{\theta}-\boldsymbol{\mu}_{\theta,i})^T \mathbf{K}_{\theta,i}^{-1} (\boldsymbol{\theta}-\boldsymbol{\mu}_{\theta,i})$, see [13] for the details. Introducing these simplifications in Eq. 40 and normalizing the distribution yields

$$ \hat{p}_{\text{CMP}}(\boldsymbol{\theta} | \mathcal{D}) \approx \sum_{i=1}^{m} w_i \frac{1}{\sqrt{(2\pi)^{d_{\theta}} \det(\mathbf{K}_{\theta,i})}} \exp\left(-\frac{1}{2}(\boldsymbol{\theta} - \boldsymbol{\mu}_{\theta,i})^T \mathbf{K}_{\theta,i}^{-1} (\boldsymbol{\theta} - \boldsymbol{\mu}_{\theta,i})\right) = p(\boldsymbol{\theta} | \mathcal{D}). \quad (41) $$

The expression in Eq. 41 is the exact representation of the parameters’ marginal posterior distribution [13].

<figure><img src="image_6.png" /><figcaption>Figure 4: Parameters’ marginal posterior distribution.</figcaption></figure>
Fig. 4 shows the resulting parameters' marginal posterior distribution for the mixture of Gaussian distribution reported in Fig. 3. The CMP approximation is exact and can correctly capture the modes and their weights. This is in contrast with the KOH method which can capture only a single mode and the FMP method which can capture all the modes but misrepresents their weights. The reader is referred to [13] for the derivation of the KOH and FMP approximations.

415

It is important to note that the CMP approximation is exact due to the well-separated nature of the modes. When the modes are not well-separated, the CMP method offers an approximate marginal posterior distribution of the parameters that may not be exact. Nonetheless, we contend that the CMP method remains preferable, as it is likely to yield a more accurate approximation compared to the KOH and FMP methods.

420

## 5. Examples

The first example in §5.1 addresses the calibration of an inadequate model, inspired by [13], featuring a bimodal posterior. The results demonstrate that the CMP method can capture both modes and their relative weights, unlike the FMP and KOH methods. The second example, in §5.2, involves the calibration of a problem that depends on three parameters and three hyperparameters, resulting in an unimodal posterior distribution. The novel method accurately captures the variability of this mode, leading to more robust predictions and correcting the false certitude effect characteristic of modular methods.

425

430

### 5.1. Example 1: Bimodal Posterior

In this example, the true function is $y(x) = x$. The experimental data consists of 11 synthetically generated observations, equally spaced within the interval $[0,1]$, based on the true model with the addition of white noise, $\epsilon \sim N(0, 0.01^2)$.

435

The computer model, which is an inadequate representation of the real process, depends on a single parameter called $\theta$. The statistical model is defined in Eq. 3, where the measurement error is modeled as a white noise process $\epsilon|_{\sigma_e} \sim N(0, \sigma_e^2)$ and the model error is represented by a Gaussian process with zero mean and a squared exponential kernel as the covariance function:

440

$$ \begin{array}{ll} \text{Computer model:} & f(x, \theta) = (1 - \theta)(x + 0.15) + x \sin(2\theta x) \\ \text{Model error covariance:} & c_{\psi_z}(x, x') = \sigma^2 \exp\left(-\frac{1}{2}\left(\frac{x - x'}{l}\right)^2\right) \end{array} \quad (42) $$
The prior for $\theta$ is considered uniform, while a set of inverse gamma functions, IG($x, \alpha, \beta$), is used for the hyperparameters. The scale and shape parameters of these distributions are the same as the ones used in [13].

The posterior distributions computed using the KOH, FMP, and CMP methods all have a dimension of 1, and they are evaluated using the analytical formula. In contrast, the full Bayesian posterior, with a dimension of 4, was sampled using an MCMC sampler. These samples were then used to construct a kernel density estimation (KDE) of the true posterior distribution, employing a rectangular kernel on the same grid as that used for the quadrature of the KOH, FMP, and CMP posteriors.

<figure><img src="image_7.png" /><figcaption>Figure 5: Slice of the un-normalized posterior distribution with $\sigma_e = 0.01$ and $l = 0.1$ (MAP value of the hyperparameters in red and the KOH value in blue). The light red background is proportional to the magnitude of the inverse of the correction factor.</figcaption></figure>

Fig. 5 shows a slice of the resulting un-normalized posterior distribution, obtained by fixing two hyperparameters: the measurement noise $\sigma_e = 0.01$ and the correlation length $l = 0.1$. The red line shows the MAP value of the hyperparameters, with variations representative of the scale of the inverse of the correction factor $|\det(S_{\theta})|$, and the blue line shows the KOH value. This figure illustrates how different modular approaches simplify the resulting distribution by evaluating it along a single line in the case of the KOH method, or along a curve for the CMP and FMP methods. The
posterior has two modes corresponding to $\theta \approx -0.1, 0.9$ with the second mode being less significant. The KOH method misses this second mode because its evaluation line does not pass through it. While other sequential approaches might provide better approximations, they are not guaranteed to capture all explanations accurately.

As expected, not including the correction factor underestimates the effect of the tails, leading to a false certitude effect. This can be inferred from Fig. 5 by observing the relative scale of the correction factor, $|\det(S_{\theta})|$, which increases at the tails.

<figure><img src="image_8.png" /><figcaption>Figure 6: Posterior distribution computed using different approximation techniques. The histogram corresponds to samples from the full-Bayesian posterior.</figcaption></figure>

Fig. 6 presents the results of the calibration. It shows that the posterior has two modes at $\theta \approx -0.1$ and $\theta \approx 0.9$, with the second mode being more significant. The KOH approximation completely misses the second mode and exhibits the false certitude effect previously mentioned. The FMP method captures both modes but significantly overestimates the weight of the first mode and also displays a similar false certitude effect. In contrast, the CMP approximation provides an excellent approximation of the true posterior.
<figure><img src="image_9.png" /><figcaption>Figure 7: Hyperparameters and correction factor used by the modular approaches.</figcaption></figure>

Fig. 7 reports the value of the hyperparameters and the Hessian used by the CMP method. The hyperparameters, see Fig. 7a, used by the KOH method correspond to the average value, computed with respect to the final posterior distribution (see Eq. 16), of the ones used by the FMP/CMP methods. Their value will be therefore influenced by the presence of a strong peak around $\theta \approx -0.1$. The magnitude of the correction factor, depicted in Fig. 7b, is significant. However, it is important to focus on relative variations rather than the absolute scale, as the CMP method’s approximation is valid up to a multiplicative constant. It is nevertheless interesting to observe that not including the Hessian would lead to a significant overestimation of the second mode.

## 5.2. Example 2: Higher Dimensional Unimodal Distribution

In this section, we will calibrate a model that relates the drag coefficient, $C_d$, of a smooth sphere moving through a fluid to its Reynolds number, Re. The first comprehensive characterization of this relationship was provided by Lapple and Shepherd in 1940 [26], and is commonly referred to as the standard drag curve. The current literature offers a wide range of tabulated experimental data and models, most of which are derived from empirical correlations. For a recent review of these resources, the interested reader is referred to [27].

For simplicity, this work will utilize only the experimental data from [26], which pertains to Reynolds numbers below 200,000, prior to the transition from laminar to fully turbulent flow. The model to be calibrated is based on the modified Stokes law (Eq. 8 in [27]) and depends on three parameters:
$$C_d = \frac{A}{\mathrm{Re}^B} + C. \qquad (43)$$

For each tested method, the posterior distribution was sampled using an MCMC algorithm. The parameter space has a dimensionality of 3, and we included a model error term with zero mean and a squared exponential kernel as the covariance function, along with an experimental error term. This results in a total of 3 hyperparameters, bringing the overall dimensionality to 6.

The priors for the model parameters are uniform, while a set of inverse-gamma functions, detailed in Tab. 1, is used for the hyperparameters.

<table><caption>Table 1: Parameters of the inverse-gamma priors of the hyperparameters and the corresponding mean and mode of the distribution.</caption><thead><tr><th>Variable</th><th>α</th><th>β</th><th>mean</th><th>mode</th></tr></thead><tbody><tr><td>σ<sub>e</sub></td><td>4</td><td>0.15</td><td>0.05</td><td>0.03</td></tr><tr><td>σ</td><td>3</td><td>0.2</td><td>0.1</td><td>0.05</td></tr><tr><td>l</td><td>3</td><td>4</td><td>2</td><td>1</td></tr></tbody></table>

Inverse gamma priors are a popular choice of hyperparameter prior as they are conjugate priors for scale parameters [1] and put zero probability mass on unwanted events, such as zero, infinite or negative values. The mean and mode of the inverse gamma distribution were chosen in order to separate the model and experimental errors' kernel strengths and to ensure that the correlation length is larger than the length scale of the experimental data.

We choose to perform the calibration of the log $C_d$ versus log Re to account for the different orders of magnitude achieved by Eq. 43.

A KDE of the posterior distribution is constructed using 10,000 independent samples. Independence is ensured by increasing the subsampling ratio to make it compatible with the correlation length.

For the FMP method, the MCMC chains for the posterior distribution were unable to reach convergence, in the sense of [19]. This happens because the prior of the hyperparameters is not present in the FMP approximation, making the resulting posterior distribution wrong and hard to sample from (see §2.3 for a detailed explanation). In [13] the authors use only uniform distributions for the prior of the hyperparameters hence this problem does not appear.
<figure><img src="image_10.png" /><figcaption>Figure 8: Comparison between different approximations of the marginal posterior.</figcaption></figure>

<table><caption>Table 2: Mean and variance of the parameters.</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">A</th><th colspan="2">B</th><th colspan="2">C</th></tr><tr><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th></tr></thead><tbody><tr><td>BAYES</td><td>31.2</td><td>6.87</td><td>0.89</td><td>0.054</td><td>0.41</td><td>0.07</td></tr><tr><td>CMP</td><td>30.9</td><td>5.65</td><td>0.89</td><td>0.050</td><td>0.41</td><td>0.06</td></tr><tr><td>KOH</td><td>29.7</td><td>2.09</td><td>0.89</td><td>0.029</td><td>0.41</td><td>0.03</td></tr></tbody></table>

Fig. 8 reports the approximation of the marginal posterior distribution of the parameters computed using the three different methods along with the one computed without the model error term. The latter was included to justify the inclusion of the model error since the resulting distribution is too overconfident and explains all errors with measurement noise. The mean and standard deviations of the marginal posteriors are also reported in Tab. 2. Both the table and the KDE plot show that the KOH method underestimates the variance of the distribution by more than 50 %. The CMP estimation of the variance, on the other hand, is much closer to the full Bayesian one. This
agrees with the previous discussion and is due to two contributions. First and foremost, the variable model error term can consider more explanations and the correction factor increases the weight of the tails, where the residual is larger.

540

<figure><img src="image_11.png" /><figcaption>Figure 9: Joint plot of the kernel strength, $\sigma$, and the correlation length, $l$.</figcaption></figure>

Fig. 9 shows the joint plot of the two hyperparameters that characterize the model error kernel. For the full Bayesian approach, the samples shown correspond to samples extracted from the marginal posterior of the hyperparameters, $p(\psi | \mathcal{D}, \mathcal{M})$. For the CMP case, they are samples of $\psi_{\text{MAP}}(\theta)$ where $\theta \sim \hat{p}_{\text{CMP}}(\theta | \mathcal{D}, \mathcal{M})$ and the blue dot is the KOH value. This shows that non-sequential methods yield a model error term that has a richer structure than the one computed by sequential approaches. An interesting correlation between these two hyperparameters can also be noted: larger model error corrections (larger $\sigma$) necessitate longer correlation lengths.

545
<figure><img src="image_12.png" /><figcaption>Figure 10: Predictive distribution of the calibrated model.</figcaption></figure>

The predictive distribution without considering the model error correction, computed using the different approximations, is reported in Fig. 10 for three values of the Reynolds number. The latter are considered representative of different scales; looking at Eq. 43, the first value of Re is very small so the non-constant term dominates, the second one is an intermediate value, and at the third one, the constant term dominates.

As visible from Fig. 10, not considering the model error or using a sequential approach has a large impact on the robustness of the posterior prediction as the results are too overconfident. Again, the CMP approximation is revealed to be an excellent one.
<figure><img src="image_13.png" /><figcaption>Figure 11: CMP prediction</figcaption></figure>

Finally, Fig. 11 shows the predictive distribution for all Reynolds numbers in the observed interval. The figure contains also the mean value of the corrected model (see Eq. 9 and Eq. 24) and the experimental data with the measurement noise. This shows that the model error can explain the discrepancy between the observed value and the model’s computation; moreover, the standard deviation of the corrected model prediction is very small. The same plot reveals, however, that there might be a confounding issue between the model and experimental error. The inferred value of the standard deviation of the measurement noise is very small, on the order of $10^{-3}$, and the entire discrepancy is explained by the model error. This problem is intrinsic to the data and the model itself since the inferred value of the correction length is large, see Fig. 9, and hence there is no confusion between the structure experimental error (which is un-correlated) and the one of the model error (whose correction length spans more than one order of magnitude). This confounding issue could be solved by using more data, possibly more log $C_d$ measurements for each log Re value.

## 6. Conclusions

In this work, we introduce a novel approach to solving the calibration problem with the inclusion of model error. The new CMP method is based on the
Bayesian calibration framework for estimating model parameters and model error hyperparameters. This method advances the current state of the art by providing a closed-form approximation of the parameters' marginal posterior distribution, effectively eliminating the need to sample the model error hyperparameters. This approximation surpasses those provided by other modular approaches, such as the KOH and FMP methods, as it consistently accounts for the uncertainty of the hyperparameters. Compared to the reference full Bayesian approach, the CMP method offers an excellent approximation.

To achieve this, we formulated a new set of assumptions on the posterior distribution that are less stringent and more general. These assumptions systematically account for the possibility that conditioning the posterior distribution on different parameter values can result in conditional distributions of the hyperparameters with varying means and variances but constant shapes. Given a specific calibration problem, it seems difficult to verify a-priori if the hypotheses of the CMP method are respected. However, an a-posteriori verification can be performed for a few parameter values by sampling the hyperparameters' conditional distribution and assess the assumption that all samples come from the same distribution up to scale-shift transformations. We tested the CMP method with two theoretical examples. In one example, we assumed the posterior distribution could be represented as a mixture of Gaussians with different means and variances. Here, the CMP approximation provided the exact answer, whereas other popular approaches either failed to capture all modes or assigned incorrect weights to them. In the other example, we examined the calibration problem without model error, where experimental noise was the only hyperparameter. Again, the CMP approximation was exact, while popular modular approaches suffered from a false certitude effect, assigning less weight to the distribution tails.

Additionally, we demonstrated the effectiveness of the CMP method with two numerical examples. In the first, targeting a 1D bimodal distribution, the CMP method accurately identified both modes and their weights, similar to the Gaussian mixture case. In the second example, focusing on a higher-dimensional unimodal distribution, the CMP method avoided the false certitude effect and provided an excellent approximation of the MAP value and its uncertainty.

We acknowledge that, at present, the benefits of the reduction in the dimensionality of the posterior are overshadowed by the significant increase in computational cost due to the introduction of an optimization step. On the one hand, this increase in computational cost might not be justified for simple models, where the full Bayesian approach is already computationally feasible, or for models with significantly more parameters than hyperparameters.
ters, in which the reduction of dimensionality is less significant. On the other hand, reducing the dimensionality of the posterior distribution is crucial for models with many hyperparameters and, even when the number of hyper-parameters is small, the CMP distribution, being a marginal distribution, is always easier to sample from than the joint one.

In the future, we plan to address the issue of the computational cost by developing a precise surrogate representation of either the optimal hyper-parameters or the approximation of the posterior using only a few sample points. These surrogates, which are inexpensive to evaluate, will then be used for parameter calibration. We are currently testing various methods, and interested readers can find some preliminary results in [28].

This work, along with the subsequent investigation into the construction of surrogate representations, will be part of a larger methodological framework for the rapid calibration of computer models. We intend to test this framework on more complex test cases.

## 7. Acknowledgments

This project has received funding from the European Union's Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101072551 (TRACES).

## 8. Supplementary material

The library for the calibration of the sphere drag example in §5.2 is available at https://github.com/omarkahol/cmp.

# Appendix A. Derivatives of the log-likelihood

This section provides a formula for the evaluation of the second derivatives of the log-likelihood, defined in Eq. 5, with respect to the hyperparameters. To do so, we start from the formula for the k-th component of the gradient of the log-likelihood [17]:

$$ \partial_{\psi_k} \log \mathcal{L}(\mathbf{D} | \boldsymbol{\theta}, \boldsymbol{\psi}) = \text{trace} \left( \frac{1}{2} (\boldsymbol{\alpha} \boldsymbol{\alpha}^T - \mathbf{K}_{\psi}^{-1}) \frac{\partial \mathbf{K}_{\psi}}{\partial \psi_k} \right), \quad (A.1) $$

with $\boldsymbol{\alpha} = \mathbf{K}_{\psi}^{-1} \mathbf{r}_{\theta}$. We proceed by taking the derivative a second time and by using the fact that the trace operator is linear,
$$ 
\begin{aligned}
\partial_{\psi_l} \partial_{\psi_k} \log \mathcal{L}(\mathbf{D} | \boldsymbol{\theta}, \boldsymbol{\psi}) = & \frac{1}{2} \operatorname{trace} \left( (\boldsymbol{\alpha} \boldsymbol{\alpha}^T - \mathbf{K}_{\psi}^{-1}) \frac{\partial^2 \mathbf{K}_{\psi}}{\partial \psi_l \partial \psi_k} \right) \\
& + \frac{1}{2} \operatorname{trace} \left( \frac{\partial \boldsymbol{\alpha} \boldsymbol{\alpha}^T}{\partial \psi_l} \frac{\partial \mathbf{K}_{\psi}}{\partial \psi_k} \right) \\
& - \frac{1}{2} \operatorname{trace} \left( \frac{\partial \mathbf{K}_{\psi}^{-1}}{\partial \psi_l} \frac{\partial \mathbf{K}_{\psi}}{\partial \psi_k} \right).
\end{aligned}
\quad (A.2) $$

Eq. A.2 shows that the Hessian of the log-likelihood is the sum of three contributions. The first contribution involves the Hessian of the observations' covariance matrix. The second contribution can be simplified in the following way:

$$ 
\begin{aligned}
\frac{\partial \boldsymbol{\alpha} \boldsymbol{\alpha}^T}{\partial \psi_l} &= 2 \operatorname{sym} \left( \frac{\partial \boldsymbol{\alpha}}{\partial \psi_l} \boldsymbol{\alpha}^T \right), \\
\frac{\partial \boldsymbol{\alpha}}{\partial \psi_l} &= \frac{\partial \mathbf{K}_{\psi}^{-1}}{\partial \psi_l} \boldsymbol{r}_{\theta} = - \mathbf{K}_{\psi}^{-1} \frac{\partial \mathbf{K}_{\psi}}{\partial \psi_l} \mathbf{K}_{\psi}^{-1} \boldsymbol{r}_{\theta}.
\end{aligned}
\quad (A.3) $$

Where we have defined with sym the operator that returns the symmetric component of a tensor, sym($X$) = ($X^T$ + $X$)/2, and used the formula for the derivative of the inverse of a matrix [17]. The third contribution can be simplified by using similar arguments. To simplify the notation, we proceed by defining $\mathbf{A}_i = \mathbf{K}_{\psi}^{-1} \frac{\partial \mathbf{K}_{\psi}}{\partial \psi_i}$. Combining Eq. A.3 and Eq. A.2 yields

$$ 
\begin{aligned}
\partial_{\psi_l} \partial_{\psi_k} \log \mathcal{L}(\mathbf{D} | \boldsymbol{\theta}, \boldsymbol{\psi}) = & \frac{1}{2} \operatorname{trace} \left( (\boldsymbol{\alpha} \boldsymbol{\alpha}^T - \mathbf{K}_{\psi}^{-1}) \frac{\partial^2 \mathbf{K}_{\psi}}{\partial \psi_l \partial \psi_k} \right) \\
& - \operatorname{trace} \left( \operatorname{sym}(\mathbf{A}_l \boldsymbol{\alpha} \boldsymbol{\alpha}^T) \frac{\partial \mathbf{K}_{\psi}}{\partial \psi_k} \right) \\
& + \frac{1}{2} \operatorname{trace}(\mathbf{A}_l \mathbf{A}_k).
\end{aligned}
\quad (A.4) $$

The Hessian of the log-likelihood with respect to the hyperparameters can be cheaply evaluated provided that the kernel and its gradients are known and that the inverse of the covariance matrix is available.

## 660 References

[1] Smith, R.C., *Uncertainty Quantification: Theory, Implementation, and Applications*, Society for Industrial and Applied Mathematics, 2013.

[2] Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., and Rubin, D.B., *Bayesian data analysis*, third edition, 2013.
665

[3] Berger, J.O. *Statistical Decision Theory and Bayesian Analysis*. Springer New York, 1985.

[4] Brynjarsdóttir, J. and O Hagan, A., Learning about physical parameters: the importance of model discrepancy, *Inverse Problems*, **30**(11), pp. 114007, 2014.

670

[5] Ling, Y., Mullins, J., and Mahadevan, S., Selection of model discrepancy priors in bayesian calibration, *Journal of Computational Physics*, **276**, pp. 665–680, 2014.

675

[6] Kennedy, M.C. and O'Hagan, A., Bayesian calibration of computer models, *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, **63**, 2001.

[7] Oliver, T.A. and Moser, R.D., Bayesian uncertainty quantification applied to rans turbulence models, Vol. 318, 2011.

680

[8] Rappel, H., Beex, L., Hale, J., Noels, L., and Bordas, S., A tutorial on bayesian inference to identify material parameters in solid mechanics, *Archives of Computational Methods in Engineering*, **27**(2), pp. 361–385, 2020.

[9] Hou, D., Hassan, I., and Wang, L., Review on building energy model calibration by bayesian inference, *Renewable and Sustainable Energy Reviews*, **143**, 2021.

685

[10] Stainforth, D., Allen, M., Tredger, E., and Smith, L., Confidence, uncertainty and decision-support relevance in climate predictions, *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, **365**(1857), pp. 2145–2161, 2007.

690

[11] Sansó, B. and Forest, C., Statistical calibration of climate system properties, *Journal of the Royal Statistical Society: Series C (Applied Statistics)*, **58**(4), pp. 485–503, 2009.

[12] Kennedy, M., Supplementary details on bayesian calibration of computer models, 2001.

695

[13] Leoni, N., Maître, O.L., Rodio, M.G., and Congedo, P.M., Bayesian calibration with adaptive model discrepancy, *International Journal for Uncertainty Quantification*, **14**(1), pp. 19–41, 2024.
[14] 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), pp. 448–466, 2004.

700

[15] Pernot, P. and Cailliez, F., A critical review of statistical calibration/prediction models handling data inconsistency and model inadequacy, *AIChE Journal*, **63**(10), pp. 4642–4665, 2017.

[16] Zhang, P., Liu, J., Dong, J., Holovati, J.L., Letcher, B., and McGann, L.E., A bayesian adjustment for multiplicative measurement errors for a calibration problem with application to a stem cell study, *Biometrics*, **68**(1), pp. 268–274, 2012.

705

[17] Rasmussen, C.E. and Williams, C.K.I., *Gaussian processes for machine learning*, 2005.

[18] *Advanced Topics in MCMC*, chapter 8, pp. 237–283. John Wiley & Sons, Ltd, 2012.

710

[19] Roy, V. Convergence diagnostics for markov chain monte carlo, 2019.

[20] Higdon, D., Gattiker, J., Williams, B., and Rightley, M., Computer model calibration using high-dimensional output, *Journal of the American Statistical Association*, **103**(482), pp. 570–583, 2008.

715

[21] 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**, pp. 417–431, 2018.

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

720

[23] Maupin, K.A. and Swiler, L.P., Model discrepancy calibration across experimental settings, *Reliability Engineering and System Safety*, **200**, pp. 106818, 2020.

725

[24] Gardner, P., Rogers, T.J., Lord, C.E., and Barthorpe, R.J., Learning model discrepancy: A gaussian process and sampling-based approach, *Mechanical Systems and Signal Processing*, 2021.

[25] Bernstein, D.S., *Matrix Mathematics: Theory, Facts, and Formulas (Second Edition)*, Princeton University Press, 2009.

730
[26] Lapple, C.E. and Shepherd, C.B., Calculation of particle trajectories, *Industrial & Engineering Chemistry, 32*(5), pp. 605–617, 1940.

[27] Review of the empirical correlations for the drag coefficient of rigid spheres, *Powder Technology*, **352**, pp. 350–359, 2019.

735

[28] Leoni, N., Bayesian inference of model error for the calibration of two-phase cfd codes, *Thesis*, Institut Polytechnique de Paris, 2022.