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

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. 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, Pietro Marco Congedo, Olivier Le Maitre, Enora Denimal Goy 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� Efficient treatment of the model error in the calibration of computer codes: the Complete Maximum a Posteriori method Omar Kahola,b,, Pietro Marco Congedoa, Olivier Le Maîtrec, Enora Denimal Goya aInria, Centre de Mathématiques Appliquées, Ecole polytechnique, IPP, Route de Saclay, 91128 Palaiseau Cedex, France bDipartimento di Scienze e Tecnologie Aerospaziali, Politecnico di Milano, Via La Masa 34, 20156 Milano, Italy cCNRS, 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 phenomena. 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 introduced 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 Maximum 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 estimate 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 Maximum a Posteriori (CMP) method, it provides a closed-form expression for the marginalization integral over the model error hyperparameters, significantly Email address: [email protected] (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, θ. These param-5 eters can be observable quantities, such as dimensional constants, or nonobservable 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 prob-10 lem. According to this perspective, the parameters, θ, 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 experi-15 mental data is observed. The two are related by Bayes’ formula; see [1–3] or §2for 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 predic-20 tion 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]. 2 The seminal work of Kennedy and O’Hagan [6] was the first to consider25 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 infinitedimensional form (such as Gaussian Processes) [1]. The model error term30 usually requires an additional parametrization; we call its parameters hyperparameters and denote them with ψto distinguish them from the model parameters, θ. 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] and35 climate [10,11]. In a full Bayesian approach, the hyperparameters, ψ, are random variables and must be learned jointly with the parameters, θ. The introduction of the model error term therefore increases the dimensionality and complexity of the calibration problem, calling for simplified solutions [12,13]. These40 solutions usually target sampling from a lower-dimensional distribution. First, sequential techniques propose to find a deterministic estimate of the hyperparameters, ψ, 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 the45 Kennedy O’Hagan (KOH) method. 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, ψand the parameters, θ, and conditioning the joint posterior distribution on that relationship.50 The marginal posterior of the parameters is the integral of the joint distribution over all possible values of the hyperparameters, ψ. 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 out-55 performs 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 method60 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 3 hyperparameters’ conditional distribution. Moreover, the FMP method does65 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 expression70 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 uncertainty75 in the model error term, from different values of the hyperparameters ψto potentially different shapes of their distribution. The paper is organized in the following way. Section 2presents a selfcontained description of the calibration problem using the Bayesian framework. Section 3is dedicated to the presentation of the CMP method. In80 §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, §6concludes the discussion and outlines possible future investigations.85 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 dis-90 cuss 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∈ U that represents the system’s state95 given a value of some measurable control variables, x∈ X. The mathematical model, M, is represented as an operator that operates on the response and is parameterized by some parameters, θ∈Θ⊂Rdθ. If uis the solution of the model, the following equation holds: M(u|x,θ) = 0 .(1) 4 The solution, u, is usually not directly measurable as it is an abstract repre-100 sentation of the response of the physical system. One usually has access to some observables,{y(i)}i=1...Nh:y(i)∈ 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 straightfor-105 ward and will be discussed later. We model the observable, y, as the sum between an observation operator, G, which operates on the solution of the model, u, and a residual term, r, y|x=G(u|x,θ) + r(x|θ).(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 is110 usually unknown but depends on the model parameters, θ, since a different choice of θwill lead to a different residual term [13]. We call a model wellspecified if there exists a value of the parameters, θ∗, that makes the residual term, r, zero: r(x|θ∗)=0. We call a model miss-specified if otherwise. In this work, we consider the case of miss-specified models.115 The term model calibration refers to the process of learning the value of the parameters θgiven some experimental data, D={(xi obs, yi obs)}i=1...N , where xi obs represents the value of the control variables and yi obs is the experimental measurement of the observable, y. In this work, we consider the case in which the experimental observations, yobs, are affected by measurement error.120 To proceed, we need to define a statistical model (and some assumptions) that explains the observations. We model the observations as a vector of Nrandom variables, Y, which is the sum of the outputs of three random functions: Y=Gθ+z+ϵ.(3) The random vector Gθcontains the output of the observation operator at125 each control variable, Gθ= (G(u|xi obs,θ))i=1...N . Following the Bayesian framework, θis a random variable and it is distributed according to the prior distribution, π(θ). The random vector z= (z(xi obs |ψz))i=1...N is the model of the residual, r(x|θ), and is called model error [6]. In this work, the model error term130 is considered to be a Gaussian Process (GP) with mean µψz(x)and covariance cψz(x,x′). The GP is assumed to have zero mean, leading to better identifiability [14]. The parameters ψzare the hyperparameters that parameterize the model error term. Similarly to the model parameters, the 5 hyperparameters are random variables and are distributed according to the135 prior distribution π(ψz). The random vector ϵcontains the measurement error at each control point, ϵ= (ϵi)i=1...N . In this work, we assume that the ϵiare i.i.d. centered Gaussian random variables with equal variance σ2 eand call π(σe)its prior distribution. The hyperparameters ψ= (ψz, σe)∈Ψ⊂Rdψare the vector of all the140 hyperparameters that parameterize the model and experimental errors. The prior distributions are updated using the available experimental data to obtain the posterior distribution, p (θ,ψ|D). The posterior distribution is related to the prior distributions by Bayes’ formula [1]: p(θ,ψ|D) = L(D|θ,ψ)π(θ)π(ψ) p(D).(4) The term L(D|θ,ψ)is the likelihood that measures the probability of ob-145 serving 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: L(D|θ,ψ) = 1 √(2π)Ndet (Kψz+σ2 eIN) exp (−1 2rT θ(Kψz+σ2 eIN)−1rθ). (5) The vector rθcontains the residuals, rθ= (yi obs −G(u|θ,xi obs))i=1...N , and Kψzis the covariance matrix of the model error term evaluated at the ob-150 servation points (Kψz)ij =cψz(xi obs,xj obs). The notation det (.)denotes the determinant of a matrix and INis the identity matrix of size N. In the literature, one can find different expressions of Eq. 5that can account also for data inconsistency [15], multiplicative experimental error [16], and model errors with non-zero mean [1].155 When multiple observables are available, the likelihood is the product of the likelihoods of each observable, L(D(1), . . . D(Nh)|θ,ψ(1), . . . ψ(Nh))= Nh ∏ i=1 L(D(i)|θ,ψ(i)),(6) where ψ(i)is the hyperparameters vector that parameterize the model and experimental error terms for the i-th observable and D(i)is the experimental data concerning the i-th observable.160 Finally, the denominator in Eq. 4is the model evidence given by: 6 p(D) = ∫Θ∫ΨL(D|θ,ψ)π(θ)π(ψ)dψdθ.(7) Integrating the joint posterior in Eq. 4over the hyperparameters yields the parameters’ marginal posterior, p(θ|D) = ∫Ψ p(θ,ψ|D)dψ.(8) The corrected model prediction of an observable at a new coordinate, x∗, denoted with y∗ corr |x∗=G(u|x∗,θ) + z(x∗|ψz), is the sum of the model165 and model error prediction. Its marginal distribution can be computed using p(y∗ corr |D) = ∫Θ∫Ψz p(y∗ corr |θ,ψz,D)p(θ,ψz|D)dθdψz,(9) where the distribution p (y∗ corr |θ,ψz,D)is normal with mean and variance specified by the Gaussian process predictive equations [17]: p(y∗ corr |θ,ψz,D) =N(y∗ corr |µpred, σ2 pred), µpred =G(u|x∗,θ) + k∗TK−1 ψzrθ, σ2 pred =cψz(x∗,x∗)−k∗TK−1 ψzk∗, (10) where (k∗)i=cψz(x∗,xi obs). This framework, together with the above-stated assumptions, constitutes the170 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 Problem175 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, one180 can use the so-called full Bayesian approach. In this case, one samples from the joint posterior distribution, p (θ,ψ|D), and then marginalizes over the hyperparameters to obtain the parameters’ posterior distribution, p (θ|D). Examples of this technique can be found in [13,14,20]. The full Bayesian 7 approach is the most general but it is also the most computationally expen-185 sive 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, ¯ ψ, and sample from the190 parameters’ posterior distribution conditioned on that value, p(θ|D,¯ ψ)∝ L(D|θ,¯ ψ)π(θ).(11) Different authors argue for different choices for the deterministic estimate of the hyperparameters, ¯ ψ, 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,195 ¯ ψ=ψKOH = arg max ψ∈Ψ∫Θ p(θ,ψ|D)dθ, = arg max ψ∈Ψ p(ψ|D). (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,200 ψMAP (θ) = arg max ψ∈Ψ p(ψ,θ|D).(13) The method then proceeds by conditioning the likelihood on the functional relationship, ψMAP (θ), and sampling from the following distribution: pFMP (θ|D)∝ L(D|θ,ψMAP (θ))π(θ).(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.205 Note that ψKOH and ψMAP (θ)can be related. This is done, initially, by using the marginalization formula, p(θ,ψ|D) = p(ψ|θ,D)p(θ|D).(15) Substituting Eq. 15 in Eq. 12 yields 8 σMAP e(θ) = √rT θrθ N+a, |det (Sθ)|=√2a+N rT θrθ∝√1 rT θrθ . (32) Calling t=σe √rT θrθ and by using Eq. 32, we can rewrite Eq. 31 as p(σe|θ,D)∝ | det (Sθ)|e−1 2t2(1 t)N+a .(33) We can recover the affine transformation by noting that ϕ=Sθ(σe−σMAP e) = t−√1 N+a=t−ϕ0.(34) Finally, substituting Eq. 34 in Eq. 33 it yields p(σe|θ,D)∝ | det (Sθ)|e −1 2(ϕ+ϕ0)2(1 ϕ+ϕ0)N+a .(35) As it is possible to see, Eq. 35 respects exactly the hypothesis of the CMP355 method. It can be rewritten as described by Eq. 18 with f(ϕ) : (−ϕ0,∞)→R+, ∝e −1 2(ϕ+ϕ0)2(1 ϕ+ϕ0)N+a . (36) 0.2 0.2 0.4 0.6 0.8 1.0 1 2 3 4 5 Figure 2: Shape function. 15 Fig. 2depicts 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. Omitting360 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.365 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 approxima-370 tions. 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.375 We consider the case in which the joint posterior, p (γ= (θ,ψ)|D), is a mixture of mGaussian distributions with weights wi: m ∑ i=1 wi= 1, means µi= (µθ,i µψ,i), and covariance matrices Ki=[Kθ,i Kθ,ψ,i KT θ,ψ,i Kψ,i ]. The posterior distribution is given by p(θ,ψ|D)∝ m ∑ i=1 wi 1 √(2π)dθ+dψdet (Ki) exp (−1 2(γ−µi)TK−1 i(γ−µi)). (37) We assume that the modes are well-separated. This translates into requir-380 ing that “the intervals: Θi={θ: (µi,θ−θ)TK−1 θ,i (µi,θ−θ)≤tdθ 95%}i=1...m (with tdθ 95% equal to the 95% quantile of the χ2law with dθ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 by385 16 p(ψ|θ,D)≈ m ∑ i=1 I(θ∈Θi)1 √(2π)dψdet (Kψ,i|θ) exp (−1 2(ψ−µψ,i|θ)TK−1 ψ,i|θ(ψ−µψ,i|θ)). (38) The conditional means can be expressed as µψ,i|θ=µψ,i +KT θ,ψ,iK−1 θ,i (θ− µθ,i)and the conditional covariances with Kψ,i|θ=Kψ,i −KT θ,ψ,iK−1 θ,i Kθ,ψ,i. The function I(θ∈Θi)is the indicator function that is equal to 1 if θ∈Θ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 Gaussian390 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, ψMAP (θ), and the correction factor, |det (Sθ)|, are well approximated by the following piecewise constant functions395 ψMAP (θ)≈ m ∑ i=1 µψ,i|θI(θ∈Θi), |det (Sθ)| ≈ m ∑ i=1 1 √det (Kψ,i|θ)I(θ∈Θi). (39) 2 1 0 1 2 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Figure 3: Example of a mixture of Gaussian distribution. The red line represents the MAP value of the hyperparameters, ψMAP (θ), and the blue line represents the KOH value. Fig. 3shows 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 17 represents the MAP value of the hyperparameters, ψMAP (θ), and the blue line represents the KOH value whose expression is reported in [13]. Using400 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, ˆ pCMP (θ|D)∝ m ∑ i=1 wi√det (Kψ,i|θ) √(2π)dθ+dψdet (Ki) exp (−1 2(θ−µθ,i µψ,i|θ−µψ,i)T K−1 i(θ−µθ,i µψ,i|θ−µψ,i)). (40) We can exploit the block structure of the covariance matrix, Ki, to express its determinant as det (Ki) = det (Kθ,i)det (Kψ,i|θ)[25]. Also, the argument405 of the exponential simplifies to −1 2(θ−µθ,i)TK−1 θ,i (θ−µθ,i), see [13] for the details. Introducing these simplifications in Eq. 40 and normalizing the distribution yields ˆ pCMP (θ|D)≈ m ∑ i=1 wi 1 √(2π)dθdet (Kθ,i) exp (−1 2(θ−µθ,i)TK−1 θ,i (θ−µθ,i))=p(θ|D). (41) The expression in Eq. 41 is the exact representation of the parameters’ marginal posterior distribution [13].410 2 1 1 2 0.2 0.4 0.6 0.8 1.0 1.2 FMP KOH CMP Figure 4: Parameters’ marginal posterior distribution. 18 Fig. 4shows 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 their415 weights. The reader is referred to [13] for the derivation of the KOH and FMP approximations. 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 of420 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. 5. Examples The first example in §5.1 addresses the calibration of an inadequate model,425 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 novel430 method accurately captures the variability of this mode, leading to more robust predictions and correcting the false certitude effect characteristic of modular methods. 5.1. Example 1: Bimodal Posterior In this example, the true function is y(x) = x. The experimental data435 consists of 11 synthetically generated observations, equally spaced within the interval [0,1], based on the true model with the addition of white noise, ϵ∼ N(0,0.012). The computer model, which is an inadequate representation of the real process, depends on a single parameter called θ. The statistical model is defined440 in Eq. 3, where the measurement error is modeled as a white noise process ϵ|σe∼ N(0, σ2 e)and the model error is represented by a Gaussian process with zero mean and a squared exponential kernel as the covariance function: Computer model: f(x, θ) = (1 −θ)(x+ 0.15) + xsin(2θx) Model error covariance: cψz(x, x′) = σ2exp −1 2(x−x′ l)2(42) 19 The prior for θis considered uniform, while a set of inverse gamma functions, IG(x, α, β), is used for the hyperparameters. The scale and shape parameters445 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 to450 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 5: Slice of the un-normalized posterior distribution with σ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. Fig. 5shows a slice of the resulting un-normalized posterior distribution, obtained by fixing two hyperparameters: the measurement noise σe= 0.01455 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θ)|, 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 of460 the KOH method, or along a curve for the CMP and FMP methods. The 20 posterior has two modes corresponding to θ≈ −0.1,0.9with 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 all465 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θ)|, which increases at the tails.470 Figure 6: Posterior distribution computed using different approximation techniques. The histogram corresponds to samples from the full-Bayesian posterior. Fig. 6presents the results of the calibration. It shows that the posterior has two modes at θ≈ −0.1and θ≈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 first475 mode and also displays a similar false certitude effect. In contrast, the CMP approximation provides an excellent approximation of the true posterior. 21 (a) Hyperparameters (MAP in red and KOH in blue) (b) Correction factor Figure 7: Hyperparameters and correction factor used by the modular approaches. Fig. 7reports 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 posterior480 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 θ≈ −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 to485 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,490 Cd, 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 empirical495 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 on500 the modified Stokes law (Eq. 8 in [27]) and depends on three parameters: 22 Cd=A ReB+C . (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. This505 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 inversegamma functions, detailed in Tab. 1, is used for the hyperparameters. Variable α β mean mode σe4 0.15 0.05 0.03 σ3 0.2 0.1 0.05 l3 4 2 1 Table 1: Parameters of the inverse-gamma priors of the hyperparameters and the corresponding mean and mode of the distribution. Inverse gamma priors are a popular choice of hyperparameter prior as they510 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.515 We choose to perform the calibration of the log Cdversus 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.520 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 uniform525 distributions for the prior of the hyperparameters hence this problem does not appear. 23 (a) A (b) B (c) C No Model Error Full Bayes CMP KOH Figure 8: Comparison between different approximations of the marginal posterior. A B C Parameter µ σ µ σ µ σ BAYES 31.2 6.87 0.89 0.054 0.41 0.07 CMP 30.9 5.65 0.89 0.050 0.41 0.06 KOH 29.7 2.09 0.89 0.029 0.41 0.03 Table 2: Mean and variance of the parameters. Fig. 8reports 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 to530 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 of535 the variance, on the other hand, is much closer to the full Bayesian one. This 24 [3] Berger, J.O. Statistical Decision Theory and Bayesian Analysis. Springer665 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. [5] Ling, Y., Mullins, J., and Mahadevan, S., Selection of model discrepancy670 priors in bayesian calibration, Journal of Computational Physics,276, pp. 665–680, 2014. [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.675 [7] Oliver, T.A. and Moser, R.D., Bayesian uncertainty quantification applied to rans turbulence models, Vol. 318, 2011. [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,680 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. [10] Stainforth, D., Allen, M., Tredger, E., and Smith, L., Confidence, un-685 certainty 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. [11] Sansó, B. and Forest, C., Statistical calibration of climate system properties, Journal of the Royal Statistical Society: Series C (Applied Statis-690 tics),58(4), pp. 485–503, 2009. [12] Kennedy, M., Supplementary details on bayesian calibration of computer models, 2001. [13] Leoni, N., Maître, O.L., Rodio, M.G., and Congedo, P.M., Bayesian calibration with adaptive model discrepancy, International Journal for695 Uncertainty Quantification,14(1), pp. 19–41, 2024. 31 [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 for705 a calibration problem with application to a stem cell study, Biometrics, 68(1), pp. 268–274, 2012. [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,710 Ltd, 2012. [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 bayesian720 analysis, with emphasis on analysis of computer models, Bayesian Analysis,4(1), pp. 119–150, 2009. [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 32 [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. [28] Leoni, N., Bayesian inference of model error for the calibration of two-735 phase cfd codes, Theses, Institut Polytechnique de Paris, 2022. 33