Bayesian Surrogate Analysis and Uncertainty Propagation

The quantification of uncertainties of computer simulations due to input parameter uncertainties is paramount to assess a model's credibility. For computationally expensive simulations, this is often feasible only via surrogate models that are learned from a small set of simulation samples. The surrogate models are commonly chosen and deemed trustworthy based on heuristic measures, and substituted for the simulation in order to approximately propagate the simulation input uncertainties to the simulation output. In the process, the contribution of the uncertainties of the surrogate itself to the simulation output uncertainties are usually neglected. In this work, we specifically address the case of doubtful surrogate trustworthiness, i.e. non-negligible surrogate uncertainties. We find that Bayesian probability theory yields a natural measure of surrogate trustworthiness, and that surrogate uncertainties can easily be included in simulation output uncertainties. For a Gaussian likelihood for the simulation data, with unknown surrogate variance and given a generalized linear surrogate model, the resulting formulas reduce to simple matrix multiplications. The framework contains Polynomial Chaos Expansions as a special case, and is easily extended to Gaussian Process Regression. Additionally, we show a simple way to implicitly include spatio-temporal correlations. Lastly, we demonstrate a numerical example where surrogate uncertainties are in part negligible and in part non-negligible.


Introduction
Uncertainty quantification of simulations has gained increasing attention, e.g. in the field of Computational Engineering, in order to address doubtful parameter choices and assess the models' credibility.Surrogate models have become a popular tool to propagate simulation input uncertainties to the simulation output, particularly for modern day applications with high computational cost and many uncertain model parameters.For that, a parametrized surrogate model (synonyms: meta-model, emulator) is learned from a finite set of simulation samples.I.e. the surrogate is a function of the uncertain simulation input parameters that is 'fitted' to the simulation output data.The quality of this fit is then judged by heuristic diagnostics, and the surrogate deemed trustworthy respectively.A key aspect of this procedure is, that the surrogate can be evaluated much faster than the simulation, and still retains a reasonable approximation to the simulation.
The simulation is then substituted with the surrogate model in order to compute the marginal probability density function of the simulation output.The simulation uncertainties are thus inferred from the surrogate model instead of the original simulation model at a significantly reduced computational effort.While this practise allows to obtain estimates on uncertainties of expensive simulations in the first place, the contribution of the uncertainty of the surrogate itself to the total simulation uncertainty is commonly neglected.In other words, the estimation of the surrogate parameters based on the finite set of simulation samples entails an additional uncertainty in the sought-for uncertainty of the simulation output.The purpose of this paper is to investigate this surrogate uncertainty as the natural measure for the surrogate's trustworthiness, and how the surrogate uncertainty affects the simulation output uncertainty.
In many cases, the surrogate uncertainty is indeed small if the heuristic diagnostics naively imply so.If the heuristic diagnostics imply that the surrogate is not trustworthy, one may resort to two options: (i) Acquire more simulation data until the surrogate is trustworthy.This is limited by the computational budget and the surrogate's convergence properties.(ii) Shrink the parameter space, e.g.omit a number of uncertain simulation parameters by assuming definite parameter values.However, in some cases (i) is not feasible and (ii) is not desired.In this contribution, we demonstrate how to include surrogate uncertainties if the user deals with a surrogate model with doubtful trustworthiness.
Popular surrogate models are Polynomial Chaos Expansions [Xiu and Karniadakis, 2005, O'Hagan, 2013, Crestaux et al., 2009] and Gaussian Process Regression [O'Hagan, 1978, Rasmussen andWilliams, 2006], the latter of which has had its renaissance recently from within the machine learning community.In this work, we assume a Gaussian likelihood for the simulation data with unknown variance and given a generalized linear surrogate model (i.e.linear in the surrogate parameters) that includes Polynomial Chaos Expansions as a special case and is easily extended to Gaussian Process Regression.Other Bayesian perspectives on Uncertainty Quantification of computer simulations with these popular surrogate models are given in [Sraj et al., 2016, O'Hagan et al., 1999, Kennedy and O'Hagan, 2000, Arnst et al., 2010, Madankan et al., 2013, Karagiannis and Lin, 2014, Lu et al., 2015, Hwai and Tan, 2015].A comprehensive collection of reviews on Uncertainty Quantification, from the point of view of computational engineering and applied mathematics, can be found in [Ghanem et al., 2017].In [O'Hagan, 2006, O'Hagan et al., 1999], a statistician's perspective is discussed.Here, we will use Bayesian Probability Theory [von der Linden et al., 2014].

Bayesian Uncertainty Quantification
We start with the general structure of uncertainty propagation problems based on surrogate models in Sec.2.1.In Sec.
2.2, we analyze a generalized linear surrogate model with a Gaussian likelihood for the simulation data with unknown surrogate variance.In Sec.2.3 we proceed to use the surrogate model to propagate the input uncertainties to the output, and show how the surrogate's uncertainties too can be included.

General structure of the problem
The goal in this paper is to quantify the uncertainties of the simulation results for the observable z (x) at different measurement points x = 1, ..., N x in the simulation domain.E.g., z could be the mechanical stress resulting from a structural analysis with a finite element simulation, where x could denote the location of the measurement probe in or on the analysed structure.z (x) depends on unknown or uncertain model parameters a = {a i } Na i=1 , which are generally inferred from experimental data d exp .Based on these data, Bayes theorem allows to determine the posterior probability density function (pdf) for a, where all background information on the experiment is subsumed in I.The implications of the background information I will be discussed later.This pdf will be assumed to be (almost) arbitrary but given in the following considerations.It usually is the result of a statistical data analysis of the foregoing experiment.This experiment could be the measurement of some material property needed for the simulation, e.g.viscosity for a computational fluid dynamics simulation.
The uncertainty of the model parameters a entails an uncertainty in the simulated observable z (x) , and the latter is determined by the marginalization rule, In the first pdf we have striked out d exp because the knowledge of a suffices to perform the simulation to obtain z (x) .If a consists of only one or two parameters (N p = 1, 2), then the numerical evaluation of the integral over the model parameters a will typically require a few dozens to hundreds of simulations.The uncertainty propagation is then done and no surrogate is needed.However this is the trivial case, and usually a will consist of way more parameters.Let's assume that a consists of, e.g., four parameters.That would imply the need for performing simulations at least 10 5 times, which is way too CPU expensive for most real problems.This can be avoided if the simulations are replaced by a surrogate model that approximates the observable z by a suitable parametrized surrogate function z sur = g(a|c), where c are yet unknown parameters.The simulation may yield the observable z (x) at different sites x in the domain, however x could also denote the time-instance in non-static problems.Clearly, the parameters of the surrogate model will also depend on those positions.So we actually have z (x) ≈ z (x) sur = g(a|c (x) ) . (3) The unknown parameters will be inferred from a suitable training data.To this end, simulations are performed for a finite set of model parameters A s = {a (i) s } Ns i=1 and the corresponding observables Z s = {z (x),(i) s } Ns,Nx i,x=1 are computed and combined in D sim = {A s , Z s }.The surrogate parameters c (x) are then inferred from D sim , and the surrogate is so constructed.We now proceed to substitute the simulation for the surrogate, z (x) → z sur , in order to solve eq. ( 2) at a significantly reduced computational cost.This implies that the background information has changed.We will denote this as Ĩ, suggesting that we take the observable z entering the integral in eq. ( 2) from the surrogate model eq.( 3) rather than from the expensive simulation.More precisely, instead of eq. ( 2) we now need to consider As far as the (second) pdf for the model parameters is concerned, we can omit the information on the training set D sim , as it does not tell us anything about the model parameters.This pdf is actually the same as that in eq. ( 1), i.e. p a | d exp , Ĩ = p a | d exp , I , as it makes no difference for the model parameters how we solve the equations underlying the simulation.In the first pdf, we can omit the information on the experiment d exp , as we only need the simulation data D sim to fix the surrogate model, which in turn defines the observable z.The first pdf can be further specified by the marginalization rule upon introducing the surrogate parameters The first pdf is uniquely fixed by the knowledge of a and C, hence D sim is superfluous.Similarly in the second pdf, where C is inferred from the training data, additional model parameters without the corresponding observables' values z, are useless.In summary we have The first pdf is rather simple.According to the background information Ĩ we will determine the observable via the surrogate model.Since the necessary parameters c (x) ∈ C and a are part of the conditional complex, the surrogate model allows only one value for the observable.That means p z (x) | C, a, Ĩ is equivalent to the probability density function for z (x) given z Finally, we have Before we can evaluate this integral, we first need to determine the two ingredients, which have their own independent significance.The last term is the result of a data analysis of a specific foregoing experiment, and will therefore not be treated here.We will suppress the background information in the following, as ambiguities should no longer occur.

Bayesian Analysis and Selection of the Surrogate Model
We recall that eq. ( 7) allows to determine the pdf for the observable based on the pdf for the model parameters, and the pdf for the parameters of the surrogate model p C | D sim , that we will determine now.To this end we have to specify the form of the surrogate model.We use the expansion in terms of basis functions Φ ν (a) and expansion coefficients c ν .No further specification is needed at this point.Without loss of generality, we will use multi-variate Legendre polynomials for the numerical examples.This expansion is very similar to the frequently used generalized Polynomial Chaos Expansion [Xiu and Karniadakis, 2005], where the polynomials Φ ν (a) are orthogonal with respect to the L 2 inner product with the prior of a, p a , as integration measure.However, here we actually consider a posterior p a | d exp that generally has no standard form, for which no standard orthogonal polynomial basis is known, and for which conditional independence of the model parameters a does not hold.The procedure has been extented to arbitrary probability measures [Oladyshkin and Nowak, 2012] and dependent parameters [Jakeman et al., 2019], but in the present context, however, these polynomials are not of primary interest and would only complicate the numerical evaluation.As outlined in section 2.1, N s simulations are performed for a set of model parameters A s = {a s } Ns i=1 and the corresponding observables Z s are computed.The theory is so far agnostic to the experimental design of these simulations, and it is therefore not of concern here.Now, we want to determine the pdf for the surrogate parameters C, which are combined in a matrix with elements C ν,x , where ν enumerates the surrogate basis functions and x enumerates the measurement positions in the domain for which the observables are computed.We abbreviated the simulation data by the quantity D sim = {A s , Z s }, where the matrix Z s has the elements (Z s ) i,x , which represent the observable z (x) at position x corresponding to the model parameter vector a (i) s .The sought-for pdf follows from Bayes' theorem The proportionality constant is not required in the ensuing considerations and we have assumed an ignorant, uniform prior for the coefficients c, i.e. p C | I = const.We note that this is also the transformation invariant Riemann prior (see Appendix B).However, any prior that is conjugate to the likelihood will retain analytical tractability.For the likelihood we need the total misfit, which is given by We assume a Gaussian type of likelihood, i.e.
with normalization Z.We have mentioned the ∆-dependence of the normalization explicitly, while the rest of of the normalization is irrelevant in the present context.Usually, the misfit entering the likelihood comes from the noise of the data.In the present case, however, there is no noise (merely a tiny numerical error), but the surrogate model is presumably not an exact description of the simulation data and ∆ covers the corresponding uncertainty.However, the uncertainty level ∆ is not known and has to be marginalized over.Along with the appropriate Jeffreys' prior, with terms independent of C subsumed in the normalization Z .For computing the mean, variance and evidence we first complete the square in eq. ( 9) to get a quadratic form in C, which can then be integrated analytically (see Appendix A).The result is and we argue that the prefactor of H −1 s is the Bayesian estimate for ∆ 2 , the variance of the Gaussian in eq. ( 10).This reasoning is similar as in [von der Linden et al., 1996].Note that Z s is a matrix of size N s × N x , containing the data vectors of length N s for each measurement site x.As shown in Appendix A, the evidence for a particular set of surrogate models is computed as where Ω N bx is the solid angle in N bx dimensions.The evidence is the probability for a surrogate model given the data.
Note that this quantity does not depend on p a | d exp , I .This is reasonable because the analysis of the experimental data should be independent of the analysis of the simulation data.However p a | d exp , I will typically be used for the experimental design of the simulation data acquisition.By comparing the evidences for different models, the user can choose a particular surrogate model or, if the results do not overwhelmingly suggest one single model, average the results for the surrogate analysis and the following uncertainty propagation over several plausible models.Note that the evidence is the pillar of a Bayesian procedure to select a surrogate model, and is distinct to the procedure of incorporating the trustworthiness or uncertainty of the surrogate in the subsequent uncertainty propagation.

Bayesian Uncertainty Propagation with Surrogate Models
Now that we have selected the surrogate model and determined the ingredients of eq. ( 7), we can determine the pdf for the observables in the light of the experimental data and the simulation results of the training set.The form in eq. ( 5) allows an easy evaluation of the mean value Similarly we obtain The covariance then follows from If we neglected the uncertainty of the surrogate, i.e.
then we retain the widely known special case of 'perfectly trustworthy' surrogates Thus, the first part in the integral of eq. ( 15) is the uncertainty of the observable due to experimental uncertainties and given the surrogate model, while the second term adds the uncertainty of the surrogate itself.The term ∆C νx ∆C ν x a is commonly neglected, but easily computed.This result also suggests a natural measure for the trustworthiness of the surrogate model, which is directly linked to the specific experiment: If the surrogate uncertainties are, on average, smaller than the experimental uncertainties by a few orders of magnitude, e.g.= 10 −3 , then they may be neglected.However is the user's choice.Note that this result does not spare the user to solve the foregoing surrogate model selection problem by e.g.computing evidences.This work only demonstrates how surrogate uncertainties can be incorporated and a practical rule when they could be neglected, given the surrogate model has already been selected before.
Figure 1: Simulation data (black dots) and simulation uncertainty (1σ) according to our Bayesian approach (red, including surrogate uncertainty) as well as the naive simulation uncertainty (blue, neglecting surrogate uncertainty).The black line is the surrogate mean.

Numerical Example
Here, we demonstrate an application where surrogate uncertainties were in part negligible and in part non-negligible.We apply our method to a computational fluid dynamics simulation of aortic hemodynamics, i.e. blood flow in an aorta resembled by the simplified geometry of an upside down umbrella stick.The simulation depends on a non-Newtonian viscosity model with four parameters a = {a 1 , a 2 , a 3 , a 4 }.The model was accompanied by viscosity measurements of human blood samples, thus determining p a | d exp , I .This posterior turned out to have a complex landscape that cannot be reasonably approximated by standard distributions.Particularly, strong correlations and multi-modality was observed, i.e. p a | d exp , I = i p a i | d exp , I .This means, that vanilla Polynomial Chaos Expansions could not be applied without an undesirable transformation to conditionally independent variables.The posterior is described in detail in [Ranftl et al., 2021].Based on p a | d exp , I , N s = 100 parameter samples a s were chosen and the simulation evaluated accordingly.The output, Z s , were the absolute values of the wall shear stress, that the blood flow exerts on the aortic wall, for N x = 10 measurement probes at different locations, each for N t = 101 time-instances equidistantly spaced over one cardiac cycle (ca.1 sec).Further details on the simulation are not relevant here, but are documented [Ranftl et al., 2021].A simulation time on the order of 150 CPU hours per sample suggested to use a surrogate for the inference.For the surrogate's basis functions, Φ ν (a), we found multi-variate Legendre polynomials up to order two sufficient.The numerical integrals were computed with Riemannian quadrature and convergence checked with successive grid refinement, however stochastic integration would work just as well.A sketch example on how to implement this procedure computationally efficient via vectorisation in parameter space can be found at https://github.com/Sranf/BayesianSurrogate_sketch.git In fig. 1, we compare the simulation uncertainty (including surrogate uncertainty) as computed with our Bayesian approach (eq.( 15)) to the naive estimate for the simulation uncertainty (without surrogate uncertainty, i.e. neglecting ∆C νx ∆C ν x a in eq. ( 15)).The surrogate uncertainties in the first half (left hand side) are relatively small, comprising only a few percent of the total uncertainty, and could possibly be neglected.In the second half (right hand side) however, the surrogate uncertainties make up to ∼ 50% of the total uncertainty.This demonstrates that simulation uncertainties inferred via surrogate models can be severely underestimated if the surrogate uncertainties are neglected, and subsequently lead to overconfidence in the simulation model.In practise, one would acquire more data in order to reduce the surrogate uncertainties, e.g. more data at later time-instances in fig. 1 is particularly promising.This was here limited not only by the computational budget, but also unpractical in that dynamic simulations require the full evaluation of all previous time instances where the surrogate is already reasonably accurate.Thus, the procedure of instead explicitly including the surrogate uncertainties here also has proven to be practical.A similar situation is to be expected for most transient simulations, as uncertainties will usually increase as time progresses.
In this work, we have assumed a Gaussian likelihood for the simulation data, with unknown variance, for a surrogate that is linear in its parameters.Surrogates that are non-linear in its parameters (e.g.neural networks) may promise higher capacity, however at the expense of losing analytical tractability of the surrogate uncertainty entirely.Other likelihood functions might be useful if further information is available, such as bounds on the observable (Gamma-or Beta-likelihood).
The result is a simple formula to incorporate surrogate uncertainties in the simulation uncertainties.This formula will be particularly useful if 'convergence' in the sense of finding the coefficients of e.g. a Polynomial Chaos Expansion is doubtful or not achievable due to a limit to the computational budget.The formula immediately suggests an intrinsic measure for the trustworthiness of the surrogate, distinct from commonly used ad-hoc diagnostics.This measure is not to be confused with the evidence, and should not be used for model selection because it would not preclude over-fitting etc.It is merely a measure for the trustworthiness of the already selected surrogate.
Let us now explore the connections of this work to Polynomial Chaos Expansions (PCE) and Gaussian Process Regression (GPR).PCE is a special case of our generalized linear surrogate model, in that the basis functions of the surrogate are chosen such that The double sum in eq. ( 15) then contracts to a single sum, and the diagonal of the term for the surrogate uncertainty, ∆C νx ∆C ν x a , survives.This is expected, in that PCE is defined such that the basis functions are uncorrelated, but still the expansion coefficients must be uncertain to a finite degree and this must carry over to the simulation uncertainty.A severe limitation of PCE is, that it is rather difficult to find basis function sets {Φ ν } that fulfill eq. ( 18), depending on p a | d exp , I .For most practical purposes, one demands (i) conditional independence of the simulation parameters, i.e. p a | d exp , I = i p a i | d exp , I , as well as (ii) simple standard distributions for p a i | d exp , I , in order to find a solution (usually a tensor-product) to eq. ( 18).Known albeit tedious work-arounds are for (i) variable transformations and numerical orthonormalisation [Jakeman et al., 2019] and for (ii) PCE-constructions for arbitrary pdfs [Oladyshkin and Nowak, 2012].Note that also [Oladyshkin and Nowak, 2012] demands (i) conditional independence of the simulation parameters.In the numerical example above, neither (i) nor (ii) were applicable.Finding a variable transformation in order to fulfill (i) or numerical construction of orthonormal basis functions can be difficult, and particularly inconvenient if sophisticated priors p a | I are being used, e.g.Jeffreys' generalized prior.An interesting alternative would be to model the input dependencies with vine copulas [Torre et al., 2019] in order to overcome the limitations of PCE addressed here.Unfortunately, no obvious vine copula was found for the here presented example.
Gaussian Process Regression would correspond to a change in the prior for z (x) in eq. ( 6) as follows where N denotes a normal distribution and K is the prior's covariance matrix and defined by the parametrized covariance function 9).By again completing the square and following the same procedure, the corresponding results for mean eq. ( 12a), variance eq. ( 12b) and evidence eq. ( 13) are then retained by a simple substitution of s | θ) is the likelihood's covariance matrix evaluated at A s for the data set Z s at given θ.eq. ( 14) and eq.( 15) would preserve their form with the substitution where the subscript θ acknowledges that the right hand side now depends on θ and s | θ) is the covariance between the training set A s and the 'test set', i.e. the integration variable a.Note that the additionally introduced hyperparameters θ would require the choice of a prior for θ and marginalization wrt θ in eq. ( 7), and subsequently also (12)(13)(14)(15).
We now discuss the implications of Ĩ in contrast to the original background information I .I contains, most importantly, that the observable z is uniquely determined by the simulation for a given set of input parameters a.A prerequisite here was, that the simulation is converged.E.g. for finite element simulations this would be a given mesh-converged spatial discretization.The proposition Ĩ additionally assumes eq. ( 3) and eq.( 6), so that it can be used to replace eq. ( 2) by eq. ( 4).Formally this means, to get from eq. (2) to eq. ( 4 3) holds and p(C | D sim ) ≈ δ(C − Ĉ), i.e. if the surrogate is indeed a good approximation and the posterior for the surrogate parameters is sharply peaked at Ĉ. Very often this posterior is not sharply peaked, then we can just gather more data until it is, or, if that is not possible, we can at least avoid overconfidence induced by neglecting these uncertainties.A numerical example of where this is the case has been demonstrated above.
We have modelled spatial correlations by introducing a location index x, and assumed that the expansion coefficients at different sites, c (x) and c (x ) , are conditionally independent.This assumption is reasonable, in that the expansion coefficients are arbitrary mathematical constructs and no physically motivated model for their correlation is known.The spatial correlation however is retained in z (x) , as was originally intended.More general models for spatial correlations can easily be implemented by substitution of δ xx with a spatial covariance matrix in eq. ( 12b).Note, that this would require an additional marginalization wrt the (typically non-linear) hyperparameters of the spatial covariance matrix.By introducing a compound index x = (x, t) and substituting x → x, we find a simple generalization to spatio-temporal correlations.This is equivalent to re-ordering spatial and temporal indices into a single sequence.While this procedure is convenient and requires only minor changes in the numerical implementation, it implicitly assumes conditional independence of spatial and temporal correlations.Analogous to above, general temporal correlations can be modelled by a substitution of δ xx in eq. ( 12b) with a temporal covariance matrix, again requiring an additional marginalization wrt the latter's hyperparameters.

Conclusions
We presented a Bayesian analysis of surrogate models and its associated uncertainty propagation problem in the context of uncertainty quantification of computer simulations.The assumptions were a generalized linear surrogate model (linear in its parameters, not the variable) and a Gaussian likelihood with unknown variance.Additionally, spatial and temporal correlations have been discussed.The result suggests a measure of trustworthiness of the surrogate by quantifying the ratio of the surrogate uncertainty to the total uncertainty, in contrast to commonly used heuristic diagnostics.The main result however is a rather simple rule to include surrogate uncertainties in the sought-for uncertainties of the simulation output.This is useful particularly for problems where the surrogate's trustworthiness is doubtful and cannot be improved.The connections to Polynomial Chaos Expansions and Gaussian Process Regression have been discussed.A numerical example demonstrated that simulation uncertainties can be significantly underestimated if surrogate uncertainties are neglected.

A Mathematical proofs
Here we want to determine norm, mean and covariance of the marginalized Gaussian (Student-t distribution) in eq. ( 11), which is In order to perform the integration, we first complete the square to get a quadratic form in C, which can then be integrated analytically, i.e. we bring the misfit χ 2 into a form that elucidates the C dependence Now, the first moment is easily obtained.Along with the variable transformation under the integral we obtain where we have used the symmetry properties of the likelihood.Next, we transform the expression for normalization based on eq. ( 22) Now we combine and reorder the double indices (ν, x) into a single index l, which turns the matrix X of dimension In this representation we have where we substituted x → H − 1 2 y.Next we introduce hyper-spherical coordinates, which leads to where Ω N bx is the solid angle in N bx dimensions.Finally, based on the substitution ρ = t • χ 2 min , we recover an identity of the Beta-function and we obtain This result is valid only for N sx > N bx , which is fulfilled in the present application.For future use we rewrite this as Finally, we calculate the covariance, based also on the compound index l = (ν, x), and by using the variable transformation in eq. ( 22).
In the last step we have used that H is a symmetric matrix.This is a very reasonable result because if the variance ∆ 2 in the Gaussian in eq. ( 10) would be known, then the covariance is ∆ 2 H −1 .Consequently, the prefactor represents the Bayesian estimate for the variance ∆ 2 based on the data.Now we go back to the original meaning of the compound index eq.( 24), i.e.H −1 ll → (H −1 s ) νν δ xx , and obtain the final result.

B The transformation invariant prior for the surrogate coefficients
Bayesian probability theory allows to rigorously and consistently incorporate any prior knowledge we have about the experiment before taking a look at the data.This knowledge shall be elicited here.Our inference must not depend on the exact parametrization.E.g. if we re-parametrize the surrogate, re-label or re-order the surrogate parameters, the surrogate still should describe the same simulation.This is reasonable because the surrogate is a purely mathematical, auxiliary construct.This rescaling-invariance is ensured by Jeffreys' generalized prior and is given by the Riemann metric R (or the determinant of the Fisher information matrix) [von der Linden et al., 2014] with multi-indices i, j = (ν, x).With the likelihood and the generalized surrogate model defined in the manuscript, the result is ∂g(a ), we replace p(z | a, I) → p(z | a, c, Ĩ), where p(z | a, I) = δ(z − z(a)), p(z | a, c, Ĩ) = δ(z − z sur (a)) = δ(z − g(a | c)).I.e.Ĩ contains in comparison to I the additional assumption that we can use the value for z as predicted/approximated by the surrogate model.It also means that we introduce additional, artificial, and usually unknown regression parameters c that need to be marginalized over.The additional uncertainty introduced by this approximation (i.e. the surrogate assumption) is encoded in p(c|D sim , Ĩ), and is correctly incorporated into the simulation observable uncertainties in eq.(15).What is important here, is that p(z | d exp , D sim , Ĩ) = p(z | d exp , ¨D sim , I) in general (the latter is computationally infeasible) but p(z | d exp , D sim , Ĩ) ≈ p(z | d exp , D sim , I) if eq. (