# Uncertainty Quantification and Propagation in Surrogate-based Bayesian
Inference

Philipp Reiser1,∗, Javier Enrique Aguilar1,2, Anneli Guthke1, Paul-Christian
Bürkner1,2

( 1Cluster of Excellence SimTech, University of Stuttgart, Germany
2Department of Statistics, TU Dortmund University, Germany
∗Corresponding author, Email<EMAIL_ADDRESS>)

###### Abstract

Surrogate models are statistical or conceptual approximations for more complex
simulation models. In this context, it is crucial to propagate the uncertainty
induced by limited simulation budget and surrogate approximation error to
predictions, inference, and subsequent decision-relevant quantities. However,
quantifying and then propagating the uncertainty of surrogates is usually
limited to special analytic cases or is otherwise computationally very
expensive. In this paper, we propose a framework enabling a scalable, Bayesian
approach to surrogate modeling with thorough uncertainty quantification,
propagation, and validation. Specifically, we present three methods for
Bayesian inference with surrogate models given measurement data. This is a
task where the propagation of surrogate uncertainty is especially relevant,
because failing to account for it may lead to biased and/or overconfident
estimates of the parameters of interest. We showcase our approach in two
detailed case studies for both linear and nonlinear modeling scenarios.
Uncertainty propagation in surrogate models enables more reliable and safe
approximation of expensive simulators and will therefore be useful in various
fields of applications.

Keywords: Surrogate Modeling, Uncertainty Quantification, Uncertainty
Propagation, Bayesian Inference, Machine Learning

##  1 Introduction

Figure 1: Overview of two-step procedure. Left: In the surrogate training step
(T-Step), training data is generated using a simulator and a surrogate model
is fitted which allows to estimate the T-posterior. Right: In the surrogate-
based inference step (I-Step), measurements along with the T-posterior are
used to infer the I-posterior.

Simulations of complex phenomena are crucial in the natural sciences and
engineering for different scenarios, e.g., for gaining system understanding,
prediction of future scenarios, risk assessment, or system design. However,
often they are based on complex ordinary differential equations or partial
differential equations which may not have closed-form solutions and may have
to be solved using expensive numerical methods. To overcome computational
overhead, the field of surrogate models (Zhu and Zabaras,, 2018; Gramacy,,
2020; Lavin et al.,, 2021) has emerged which provide fast approximations of
computationally expensive simulation. Examples are polynomial chaos expansion
(Wiener,, 1938; Sudret,, 2008; Oladyshkin and Nowak,, 2012; Bürkner et al.,,
2023), Gaussian processes (Rasmussen and Williams,, 2005) or neural networks
(Goodfellow et al.,, 2016). Recently there has been a great interest in
applying surrogate models in relevant areas, for example in hydrology
(Mohammadi et al.,, 2018; Tarakanov and Elsheikh,, 2019; Zhang et al.,, 2020),
in fluid dynamics (Meyer et al.,, 2021), in climate prediction (Kuehnert et
al.,, 2022), or in systems biology (Renardy et al.,, 2018; Alden et al.,,
2020). Furthermore, great methodological advances have been made in the field
of surrogate modeling, for example, the incorporation of physical knowledge
(Raissi et al.,, 2019; Li et al.,, 2021; Brandstetter et al.,, 2023) or the
combination of surrogate models and simulation-based inference (Radev et al.,,
2023).

Despite these advances, a major remaining challenge is the trustworthiness and
reliability of the surrogate. Consequently, it is crucial to quantify
uncertainties associated with surrogate modeling, e.g., caused by limited
training data or inflexibility of the surrogate. To estimate the uncertainty
in surrogate model parameters, several methods for uncertainty quantification
(UQ) have been developed (e.g., Shao et al.,, 2017; Bürkner et al.,, 2023).
Often however, UQ is insufficient on its own: To account for surrogate
uncertainties in subsequent (surrogate-based) inference tasks, uncertainty
propagation (UP) is necessary as well (Smith,, 2013; Lavin et al.,, 2021;
Psaros et al.,, 2023).

Using probability theory as its fundamental basis, Bayesian statistics
provides a rigorous way for both UQ and UP (Gelman et al.,, 2013; McElreath,,
2020; Bürkner et al.,, 2023). One important use-case for UP in surrogate
modeling arises from the ”forward” problem, in which input parameters are
uncertain and the goal is to propagate them through the surrogate while
accounting for its uncertainty to compute a reliable output. For example,
Ranftl and von der Linden, (2021) proposed a method for UP in the forward
problem, for restricted cases where conjugate analytic models are available.
Similarily, Zhu and Zabaras, (2018) quantified uncertainty in neural network
surrogate outputs using approximate Bayesian inference. In the ”inverse”
problem, when surrogates are used for Bayesian inference of parameters given
observed data (e.g., Marzouk et al.,, 2007; Marzouk and Xiu,, 2009; Zeng et
al.,, 2012; Laloy et al.,, 2013; Li and Marzouk,, 2014), rigorous propagation
of the surrogate uncertainty becomes even more relevant and challenging. For
example, Zhang et al., (2020) applied surrogate-based Bayesian inference to
hydrological systems and propagated parts of the surrogate uncertainty
assuming normal surrogate posteriors. We review these methods in more detail
in Section 2.2.3.

In this paper, our primary focus lies on UP when solving probabilistic inverse
problems via surrogate models. Existing methods proposed for the same
challenge are scarce and propagate the surrogate-based uncertainty only in a
selective and simplified manner, which leaves a lot of room for both improved
theory and improved practical methods. This not only concerns UP itself, but
also diagnostic methods to assess the validity of the resulting inference.
This paper aims at addressing these challenges from a fully probabilistic
(Bayesian) perspective. Concretely, we make the following contributions: (i)
We introduce a formal framework for surrogate-based Bayesian inference that
specifies and explains how to quantify all relevant uncertainties. (ii) We
study three distinct methods for UP within a two-step inference procedure,
which consists of a surrogate training step (T-Step) and a surrogate-based
inference step (I-Step). For a high-level overview, see Fig. 1. (iii) We adapt
existing simulation-based procedures to validate the uncertainty calibration
achieved via the different UP approaches. Finally, we evaluate our methods in
two detailed case studies.

##  2 Method

In the following, we propose a framework for uncertainty propagation (UP) in
surrogate-based Bayesian inference. This framework is applicable to general
surrogate models, e.g., linear regression, polynomial chaos expansion,
Gaussian processes or neural networks. The framework consists of a two-step
procedure with a surrogate training step (T-Step) and an inference step
(I-Step), an overview of all uncertainties occurring in surrogate-based
inference, several UP schemes of selected uncertainties, and their evaluation.

###  2.1 Two-Step Procedure

We propose a two-step procedure for uncertainty propagation of surrogate
models which consists of (i) the Surrogate-Training Step (T-Step), where
training data is generated using a simulator and (ii) the Inference-Step
(I-Step) where the trained surrogate is used to infer a quantify of interest,
as illustrated in the graphical model shown in Fig. 2.1. Even though the
general setup may look relatively simple, it becomes challenging to quantify
and propagate all occurring uncertainties in a statistically rigorous manner.
In the following, we differentiate between aleatoric (irreducible) and
epistemic (reducible) uncertainty; for a review see Hüllermeier and Waegeman,
(2021); Gruber et al., (2023). In Table 1 we summarize all relevant
uncertainties occurring in the two-step procedure to be discussed in detail
below.

Figure 2: Graphical model of the T-Step and I-Step. Left: In the T-Step, the
$N_{T}$ observed quantities are simulation parameters $\omega_{T}$ and
simulation output $y_{T}$, and the unknowns are the surrogate parameters
$\theta$. Right: In the I-Step, measurement data $y_{I}$ is observed $N_{I}$
times and $S$ posterior samples of $\theta$ are propagated from the T-Step.
The unknowns to be inferred are the simulation parameters $\omega_{I}$.

####  2.1.1 First Step: Training the Surrogate (T-Step)

In the first step, we focus on training a surrogate model using artificial
data generated from the complex simulator we seek to approximate. This
involves computing a posterior over the surrogate parameters to account for
the induced uncertainties.

##### Simulator

We consider a given (potentially stochastic) simulator $\mathcal{M}$, which is
an arbitrarily complex model, for example describing a physical or biological
process, and which is hard to evaluate. Given an input simulation parameter
$\omega_{T}$, we obtain the output response from the simulator as follows:

$\displaystyle y_{T}=\mathcal{M}(\omega_{T};e_{S})\quad\text{with}\quad
e_{S}\sim p(e_{S}\mid\sigma_{S}),$ (1)

where $e_{S}$ is the noise of the simulation drawn from a simulator noise
distribution $p(e_{S}\mid\sigma_{S})$ with hyperparameters $\sigma_{S}$ that
describe the aleatoric (irreducible) uncertainty of the simulation. Note, that
the simulation model itself does not necessarily have to be stochastic; it
could be a deterministic physics-based or conceptual model that is
complemented with a stochastic representation of measurement noise. We
explicitly choose the subscript $S$ to highlight that the noise stems from the
simulation and not from the training process. This setup induces the true
generating distribution $p(y_{T}\mid\omega_{T},\sigma_{S})$ of simulation
responses $y_{T}$ given input simulation parameters $\omega_{T}$ and noise
hyperparamters $\sigma_{S}$.

##### Surrogate Model

A surrogate model $\widetilde{\mathcal{M}}$ is a statistical model that aims
to approximate the simulator $\mathcal{M}\approx\widetilde{\mathcal{M}}$ while
being computationally more efficient to evaluate. The parametric form of the
surrogate is defined by a set of surrogate approximation parameters $c$. Given
input parameters $\omega_{T}$, surrogate approximation parameters $c$, and the
simulator noise distribution $p(e_{S}\mid\sigma_{S})$ we can calculate the
surrogate response:

$\displaystyle\tilde{y}_{T}=\widetilde{\mathcal{M}}(\omega_{T},c,e_{S})\quad\text{with}\quad
e_{S}\sim p(e_{S}\mid\sigma_{S}).$ (2)

##### Surrogate Approximation Error

A surrogate model with a fixed architecture will always be misspecified, if
the true simulator is not included in the class of surrogate models that we
define. This is the standard use case, since we are specifically tailoring the
surrogate model to be much simpler than the simulator. The uncertainty of the
surrogate approximation parameters $c$ only captures the (epistemic)
uncertainty due to limited training data. To additionally account for the
approximation error of the surrogate with respect to the simulator, caused by
limited expressibility of the surrogate, we introduce an additional error term
$e_{A}$. We assume $e_{A}$ to follow a distribution $p(e_{A}\mid\sigma_{A})$
with hyperparameters $\sigma_{A}$. This distribution captures aleatoric
uncertainty which is not reducible by more training data. We can then model
the true response $y_{T}$ via a (potentially unknown) function $\tilde{f}$
that takes the output of the surrogate $\tilde{y}_{T}$ and the approximation
error $e_{A}$:

$\displaystyle y_{T}=\tilde{f}(\tilde{y}_{T},e_{A})\quad\text{with}\quad
e_{A}\sim p(e_{A}\mid\sigma_{A}).$ (3)

In practice, we might assume a simple additive error:
$y_{T}=\tilde{y}_{T}+e_{A}$ and a normal distribution for $e_{A}$, but our
framework is agnostic to these choices. For the sake of readability, we use
$\theta=\\{c,\sigma_{A}\\}$ to combine all trainable surrogate parameters in a
single vector. Together, this implies a surrogate likelihood of the simulator
responses $y_{T}$:

$y_{T}\sim p(y_{T}\mid\omega_{T},\sigma_{S},\theta).$ (4)

##### Surrogate Training

To train the surrogate, we use the simulator $\mathcal{M}$ to generate
training data
$\mathcal{D}_{T}=\\{\omega_{T}^{(i)},\sigma_{S}^{(i)},y_{T}^{(i)}\\}_{i=1}^{N_{T}}$
consisting of $N_{T}$ inputs $(\omega_{T}^{(i)},\sigma_{S}^{(i)})$ and
corresponding outputs $y_{T}^{(i)}\sim
p(y_{T}^{(i)}\mid\omega_{T}^{(i)},\sigma_{S}^{(i)})$. The goal of the T-Step
is to fit the parameters $\theta$ of the surrogate $\widetilde{\mathcal{M}}$
to approximate the data distribution implied by $\mathcal{M}$.

To train our surrogate parameters $\theta$, we perform Bayesian inference
using the fast-to-evaluate surrogate likelihood
$p(y_{T}\mid\omega_{T},\sigma_{S},\theta)$ and a (potentially non-informative)
prior $p(\theta)$. We obtain the joint posterior distribution over all
surrogate parameters given the simulation training data $\mathcal{D}_{T}$ as:

$\displaystyle
p(\theta\mid\mathcal{D}_{T})\propto\prod_{i=1}^{N_{T}}p(y_{T}^{(i)}\mid\omega_{T}^{(i)},\sigma_{S}^{(i)},\theta)\,p(\theta).$
(5)

This joint posterior describes the epistemic (reducible) uncertainty in the
surrogate model parameters. In the case of infinite training data, i.e.
$N_{T}\to\infty$ the posterior $p(\theta\mid\mathcal{D}_{T})$ converges to a
point mass under regularity conditions (van der Vaart,, 2000). To approximate
this posterior we can use sampling-based algorithms, such as Markov chain
Monte Carlo (MCMC) (Robert and Casella,, 2005), and represent the posterior in
the form of $S$ posterior samples $\\{\theta^{(1)},\dots,\theta^{(S)}\\}\sim
p(\theta\mid\mathcal{D}_{T})$.

####  2.1.2 Second Step: Inference on Real Data (I-Step)

In the second step, real-world measurement data is given and our goal is to
infer the unknown quantities of interest using the previously trained
surrogate model as an efficient replacement of the complex simulator.

##### Measurement Model

We assume that the (implicit) real-world data generator is well described by
the simulator $\mathcal{M}$, but with a potentially different measurement
noise $e_{I}$. Measurement data $y_{I}$ is then generated from an unknown
underlying parameter $\omega_{I}$:

$\displaystyle y_{I}=\mathcal{M}(\omega_{I};e_{I})\quad\text{with}\quad
e_{I}\sim p(e_{I}\mid\sigma_{I}),$ (6)

where the measurement noise $e_{I}$ is drawn from a distribution
$p(e_{I}\mid\sigma_{I})$ with hyperparameters $\sigma_{I}$ that describe the
aleatoric uncertainty of the measurement. This induces a generating
distribution $p(y_{I}\mid\omega_{I},\sigma_{I})$ of measurements $y_{I}$ given
inputs $\omega_{I}$ and noise hyperparameters $\sigma_{I}$. In contrast to the
simulator setup, we only have access to the measurements $y_{I}$ but the
underlying true input $\omega_{I}$ and $\sigma_{I}$ are unknown.

##### Inference

Given a set of real-world measurement data
$y_{I}=\\{y_{I}^{(i)}\\}^{N_{I}}_{i=1}$ with $y_{I}^{(i)}\sim
p(y_{I}^{(i)}\mid\omega_{I},\sigma_{I})$ for $i=1,\ldots,N_{I}$, our goal is
to infer the posterior of the unknown parameters $\omega_{I}$, which
constitute our primary quantity of interest. Additionally, we can also infer
$\sigma_{I}$ although it is only of secondary interest. To simplify the
presentation, we drop $\sigma_{I}$ from the notation in the following.

If the simulator $\mathcal{M}$ had a tractable and easy-to-evaluate likelihood
function $p(y_{I}\mid\omega_{I},\mathcal{M})$, we could simply calculate the
posterior of the unknown simulation parameters $\omega_{I}$ given the
measurement data $y_{I}$ and a prior $p(\omega_{I})$:

$\displaystyle p(\omega_{I}\mid y_{I},\mathcal{M})\propto
p(y_{I}\mid\omega_{I},\mathcal{M})\,p(\omega_{I}).$ (7)

However, for many real-world simulators, this likelihood is unavailable or
highly cumbersome to evaluate (see Section 1), which is why we seek to replace
it with simpler likelihood of the previously trained surrogate
$\widetilde{\mathcal{M}}$. To incorporate the uncertainties of the surrogate
parameters $\theta=(c,\sigma_{A})$ from the T-Step into our real-world
inference, we need to somehow propagate their T-Step posterior
$p(\theta\mid\mathcal{D}_{T})$ to the I-step. From the perspective of the
I-Step, the uncertainty in $p(\theta\mid\mathcal{D}_{T})$ becomes aleatoric
since it is no longer reducible. The main question how to propagate
$p(\theta\mid\mathcal{D}_{T})$ turns out to have multiple answers – even
multiple ones fully justified by probability theory. Each of these answers
corresponds to an uncertainty propagation method $u\in\mathcal{U}$, for which
we can obtain a surrogate-based posterior $p(\omega_{I}\mid
y_{I},\widetilde{\mathcal{M}},u)$ as detailed below.

Table 1: Uncertainties in the two-step procedure.  Uncertainty | T-Step | I-Step | Short name
---|---|---|---
Simulator noise | aleatoric | - | -
hyperparameter $\sigma_{S}$ |  |  |
Surrogate approx. | aleatoric | aleatoric | T-aleatoric uncertainty
error hyperparameter $\sigma_{A}$ |  |  |
Surrogate parameter | epistemic | aleatoric | T-epistemic uncertainty /
posterior $p(c,\sigma_{A}\mid\mathcal{D}_{T})$ |  |  | T-posterior
Measurement noise | - | aleatoric | -
hyperparamter $\sigma_{I}$ |  |  |
Simulator-based | - | epistemic | -
posterior $p(\omega_{I}\mid y_{I};\mathcal{M})$ |  |  |
Surrogate-based | - | epistemic & | I-posterior
posterior $p(\omega_{I}\mid y_{I};\widetilde{\mathcal{M}},u)$ |  | aleatoric* |
* Depends on uncertainty propagation method $u$.

###  2.2 Uncertainty Propagation in the Two-Step Procedure

In the following, we present four different uncertainty propagation methods to
perform the surrogate-based inference. These methods are namely (i) a Point
Estimate, (ii) the Expected-Posterior (E-Post), (iii) the Expected-Likelihood
(E-Lik), and (iv) the Expected-Log-Likelihood (E-Log-Lik). Accordingly, the
set of considered propagation methods is given by
$\mathcal{U}=\\{\text{Point},\text{E-Post},\text{E-Lik},\text{E-Log-Lik}\\}$.
The three latter approaches propagate the uncertainty from the T-Step using
the full posterior or samples from the posterior. In E-Post, we incorporate
the T-Step posterior $p(\theta\mid\mathcal{D}_{T})$ directly into the
posterior of the I-Step, while in the E-Lik and E-Log-Lik, the T-Step
posterior is incorporated into the likelihood
$p(y_{I}\mid\omega_{I},\widetilde{\mathcal{M}},u)$ of the I-step. In the
following, all posteriors are calculated using the surrogate model
$\widetilde{\mathcal{M}}$ and we will omit the dependency of
$\widetilde{\mathcal{M}}$ in the posterior: $p(\omega_{I}\mid
y_{I},u)=p(\omega_{I}\mid y_{I},\widetilde{\mathcal{M}},u)$ to improve
readability.

####  2.2.1 Point Estimate

First, we describe the I-Step using only a point estimator $\hat{\theta}$ of
the T-posterior $p(\theta\mid\mathcal{D}_{T})$, e.g., the posterior mean,
median, or mode. In this case, we condition the posterior of $\omega_{I}$ on
$\hat{\theta}$, i.e., we reduce the full posterior to a point estimate without
epistemic uncertainty. This is not a method we advocate for, but rather serves
as a baseline to compare more sophisticated methods against. The posterior
probability of $\omega_{I}$ given measurement data $y_{I}$ and $\hat{\theta}$
is then (up to normalizing constants):

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{Point})$ $\displaystyle\propto
p(y_{I}\mid\omega_{I},\hat{\theta})p(\omega_{I}).$ (8)

In the following, for all uncertainty propagation methods, we will assume that
$y_{I}$ consists of conditional i.i.d. measurements $y_{I}^{(i)}$, such that
the likelihood factorizes easily. This assumption is not necessary for our
framework but simplifies notation and allows for more efficient computation.
With this assumption, the (unnormaliezd) log I-posterior is

$\displaystyle\log p(\omega_{I}\mid
y_{I},u=\text{Point})=\sum_{i=1}^{N_{I}}\log
p(y_{I}^{(i)}\mid\omega_{I},\hat{\theta})+\log p(\omega_{I}).$ (9)

This log-posterior can easily be specified in a probabilistic programming
language such as Stan (Carpenter et al.,, 2017) and we can then sample from
the posterior with MCMC, leading to $K$ posterior samples
$\\{\omega_{I}^{(1)},\dots,\omega_{I}^{(K)}\\}\sim p(\omega_{I}\mid
y_{I},u=\text{Point})$. The benefit of the Point method is its comparably fast
evaluation, since we only use a point estimate of the T-posterior and
therefore can quickly evaluate the I-likelihood
$p(y_{I}\mid\omega_{I},\hat{\theta})$. While we are completely neglecting the
epistemic uncertainty contained in the T-posterior, we propagate (a point
estimate of) the surrogate approximation error parameter $\sigma_{A}$ through
$\hat{\theta}$. We expect the Point I-posterior $p(\omega_{I}\mid
y_{I},u=\text{Point})$ to be overconfident (too narrow); unless we have a
sufficiently large amount of simulation training data $\mathcal{D}_{T}$ such
that $p(\theta\mid\mathcal{D}_{T})$ converges to a point mass and then
corresponds exactly to the point estimator $\hat{\theta}$.

####  2.2.2 Expected-Posterior

Next, we present the Expected-Posterior (E-Post), a method that propagates
both aleatoric and epistemic uncertainty in the T-posterior by marginalizing
over the surrogate parameter posterior $p(\theta\mid\mathcal{D}_{T})$ in the
posterior of $\omega_{I}$:

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Post})=\int p(\omega_{I}\mid
y_{I},\theta)p(\theta\mid\mathcal{D}_{T})d\theta.$ (10)

Since we only have samples $\theta^{(s)}\sim p(\theta\mid\mathcal{D}_{T})$
from the T-step, we compute the Monte Carlo (MC) approximation of E-Post as:

$\displaystyle p(\omega_{I}\mid
y_{I},u=\text{E-Post})\stackrel{{\scriptstyle\text{MC}}}{{\approx}}\frac{1}{S}\sum_{i=1}^{S}p(\omega_{I}\mid
y_{I},\theta^{(s)}),$ (11)

which now is a finite mixture model with equal weights. This can easily be
implemented in a probabilistic program by fitting a separate model for each
T-posterior draw $\theta^{(s)}\sim p(\theta\mid\mathcal{D}_{T})$ using the
point I-posterior approach in Eq. (9), which results in $K$ draws of
$\omega_{I}$, i.e., $\\{\omega_{I}^{(s,1)},\dots,\omega_{I}^{(s,K)}\\}$. The
combination of these draws across all $s=1,\ldots,S$ then represents an MC-
estimate of the E-Post I-posterior as per Eq. (11), with a total of $S\cdot
K^{*}$ posterior draws.

##### Related work.

The idea of constructing a posterior via the aggregation of multiple posterior
distributions, each approximated via samples, has been explored in multiple
places in the literature. In the BayesBag method (Waddell et al.,, 2002;
Douady et al.,, 2003; Bühlmann,, 2014; Huggins and Miller,, 2020), posteriors
are obtained from bootstrapped copies (Efron,, 1979; Breiman,, 2004) of the
original dataset and subsequently averaging the resulting bootstrapped
posteriors. Similarily, in the context of missing value imputation (Little and
Rubin,, 2019), this approach has been used to combine models fitted on
multiple imputed datasets (Bürkner,, 2017, 2018). In both cases, models are
fitted to different datasets whereas in our E-Post method, the data is the
same but the model itself changes.

####  2.2.3 Expected-Likelihood

Next, we present the Expected-Likelihood (E-Lik) approach, where we
marginalize over the T-posterior $p(\theta\mid\mathcal{D}_{T})$ in the
I-likelihood:

$\displaystyle p(y_{I}\mid\omega_{I},u=\text{E-Lik})=\int
p(y_{I}\mid\omega_{I},\theta)p(\theta\mid\mathcal{D}_{T})d\theta.$ (12)

The full I-posterior of E-Lik is then given by

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Lik})$ $\displaystyle\propto
p(y_{I}\mid\omega_{I},u=\text{E-Lik})p(\omega_{I}).$ (13)

Assuming a factorizable likelihood, we compute the log I-posterior over the
whole dataset $y_{I}$ (without normalizing constant) as:

$\displaystyle\log p(\omega_{I}\mid
y_{I},u=\text{E-Lik})=\log\int\prod_{i=1}^{N_{I}}p(y_{I}^{(i)}\mid\omega_{I},\theta)p(\theta\mid\mathcal{D}_{T})d\theta+\log
p(\omega_{I}),$ (14)

for which we can obtain an MC approximation using draws $\theta^{(s)}\sim
p(\theta\mid\mathcal{D}_{T})$:

$\displaystyle\log p(\omega_{I}\mid
y_{I},u=\text{E-Lik})\stackrel{{\scriptstyle\text{MC}}}{{\approx}}\log\left(\frac{1}{S}\sum_{s=1}^{S}\prod_{i=1}^{N_{I}}p(y_{I}^{(i)}\mid\omega_{I},\theta^{(s)})\right)+\log
p(\omega_{I}).$ (15)

This representation has the problem that the product over likelihood
components becomes numerically unstable as $N_{I}$ grows larger, since the log
operator cannot simply be pulled into sum over draws. To circumvent this, we
use the log-sum-exp trick, i.e. calculate $\log\sum_{s=1}^{S}\exp(\log
p(y_{I}\mid\omega_{I},\theta^{(s)}))$, which has a numerically stable
implementation (Carpenter et al.,, 2017). For given $\theta^{(s)}$, the joint
log likelihood is then again a simple sum: $\log
p(y_{I}\mid\omega_{I},\theta^{(s)})=\sum_{i=1}^{N_{I}}\log
p(y_{I}^{(i)}\mid\omega_{I},\theta^{(s)})$. In contrast to E-Post, E-Lik is
expressed as a single probabilistic program and we can use MCMC to approximate
its I-posterior with $K$ samples
$\\{\omega_{I}^{(1)},\dots,\omega_{I}^{(K)}\\}$.

More general, in contrast to E-Post, which marginalizes over $\theta$ in the
I-posterior $p(\omega_{I}\mid y_{I},\theta)$, the E-Lik method marginalizes
over $\theta$ in the I-likelihood $p(y_{I}\mid\omega,\theta)$. E-Lik and
E-Post are not identical (see Section B.1 for a counterexample), but we
demonstrate in Section 3 that they usually yield very similar results.

##### Related work.

To our knowledge, E-Lik constitutes a novel method for full uncertainty
propagation. The approach used in Zhang et al., (2020) appears related,
although the details of their method are insufficiently described in the paper
for a definitive assessment. Based on our understanding, their method can be
seen as a special case of E-Lik, where the T-posterior is assumed to be normal
and surrogate approximation error is ignored. Further, E-Lik is related to
important quantities outside the area of UP: In particular, the log-likelihood
constructed in E-Lik resembles the expected log predictive density (ELPD), a
popular measure for predictive performance, which integrates the likelihood
over the posterior of the same model before taking the logarithm outside the
expectation (Vehtari and Ojanen,, 2012; Vehtari et al.,, 2017; Bürkner et
al.,, 2023). What is more, in the context of meta-analysis, Blomstedt et al.,
(2019) proposed a method to combine posteriors resulting from different
studies. While their sources of uncertainty are different than in our case,
they also integrate them with an expected likelihood approach, rendering at
least the core idea related to E-Lik.

####  2.2.4 Expected-Log-Likelihood

In the Expected-Log-Likelihood (E-Log-Lik) approach, instead of marginalizing
over the likelihood as in E-Lik, we marginalize over the T-posterior in the
log-likelihood:

$\displaystyle p(y_{I}\mid\omega_{I},u=\text{E-Log-Lik})=\exp\left(\int\log
p(y_{I}\mid\omega_{I},\theta)\,p(\theta\mid\mathcal{D}_{T})d\theta\right).$
(16)

The posterior of the E-Log-Lik for $\omega_{I}$ is then simply

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Log-Lik})$
$\displaystyle\propto p(y_{I}\mid\omega_{I},u=\text{E-Log-Lik})p(\omega_{I}).$
(17)

When assuming a factorizable likelihood and ignoring the normalizing constant,
the log posterior becomes

$\displaystyle\log p(\omega_{I}\mid y_{I},u=\text{E-Log-Lik})$
$\displaystyle=\int\log\left(\prod_{i=1}^{N_{I}}p(y_{I}^{(i)}\mid\omega_{I},\theta)\right)p(\theta\mid\mathcal{D}_{T})d\theta+\log
p(\omega_{I})$ (18) $\displaystyle=\int\sum_{i=1}^{N_{I}}\log
p(y_{I}^{(i)}\mid\omega_{I},\theta)\,p(\theta\mid\mathcal{D}_{T})d\theta+\log
p(\omega_{I}),$ (19)

The integral is readily approximated via draws from the T-posterior:

$\displaystyle\log p(\omega_{I}\mid y_{I},u=\text{E-Log-
Lik})\stackrel{{\scriptstyle\text{MC}}}{{\approx}}\frac{1}{S}\sum_{s=1}^{S}\sum_{i=1}^{N_{I}}\log
p(y_{I}^{(i)}\mid\omega_{I},\theta^{(s)})+\log p(\omega_{I}),$ (20)

which can be fitted with MCMC leading to $K$ draws
$\\{\omega_{I}^{(1)},\dots,\omega_{I}^{(K)}\\}$ in the I step. We can further
rewrite the MC approximation of the log-likelihood as

$\displaystyle\log p(y_{I}\mid\omega_{I},u=\text{E-Log-
Lik})=\sum_{s=1}^{S}\sum_{i=1}^{N_{I}}\log\left(p(y_{I}^{(i)}\mid\omega_{I},\theta^{(s)})^{\frac{1}{S}}\right).$
(21)

This shows that the E-Log-Lik can be interpreted as power-scaling the
likelihood components with equal weights $1/S$ (Geyer,, 1991; Kallioinen et
al.,, 2022), which readily generalizes to unequal weights as we illustrate in
Section 2.2.5. However, in contrast to E-Lik and E-Post, this does not
strictly follow rules of probability theory as we integrate over a probability
measure in the log space. Alternatively to E-Log-Lik, we could set up an
Expected Log Posterior (E-Log-Post) as

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Log-
Post})\propto\exp\left(\int\log p(\omega_{I}\mid
y_{I},\theta)\,p(\theta\mid\mathcal{D}_{T})d\theta\right),$ (22)

which however is equivalent to E-Log-Lik as we show in Appendix B.2.

##### Related work.

The computation of the E-Log-Lik is conceptually similar to a single E-Step in
an Expectation-Maximization (EM) algorithm (Dempster et al.,, 1977), where the
log likelihood is integrated over a discretized space of latent variable
values (instead of T-posterior draws as is done here). Furthermore, the Gibbs
loss, a computational convenient although uncommon measure of predictive
performance, also calculates an expectation over the log-likelihood (Vehtari
and Ojanen,, 2012; Bürkner et al.,, 2023).

####  2.2.5 Clustering of the T-Posterior Draws

To speed up computation of the I-Step while minimizing loss of information, we
can reduce the number of propagated T-posterior draws via clustering (see
Piironen et al., (2020) for a related use case in the context of variable
selection). For this purpose, any clustering algorithm can in principle be
used, for example, KMeans clustering (MacQueen,, 1967; Lloyd,, 1982). We apply
the clustering algorithm to the set of T-posterior draws of the surrogate
parameters $\\{\theta^{(s)}\\}_{s=1}^{S}$ to get cluster centroids
$\\{\mu^{(l)}\\}_{l=1}^{L}$. For each of the cluster centroids we additionally
store weights $\\{\alpha^{(l)}\\}_{l=1}^{L}$, where $\alpha^{(l)}$ is the
percentage of draws associated with the cluster. The cluster centroids
together with the corresponding weights allow for a reliable and
computationally efficient processing of the T-posterior draws as we show in
Section 3.

Clustering and re-weighting can be easily applied in all UP methods by
replacing $\theta^{(s)}$ with $\mu^{(l)}$ as well as replacing the equal
weight $1/S$ with $\alpha^{(l)}$ after moving it inside the sum. For example,
when using E-Log-Lik, the MC-approximated integral becomes:

$\displaystyle\log p(y_{I}\mid\omega_{I},u=\text{E-Log-Lik})$
$\displaystyle\stackrel{{\scriptstyle\text{MC}}}{{\approx}}\frac{1}{S}\sum_{s=1}^{S}\sum_{i=1}^{N_{I}}\log(p(y_{I}^{(i)}\mid\omega_{I},\theta^{(s)})$
(23) $\displaystyle\approx\sum_{l=1}^{L}\sum_{i=1}^{N_{I}}\alpha^{(l)}\log
p(y_{I}^{(i)}\mid\omega_{I},\mu^{(l)}).$ (24)

####  2.2.6 Parallelization

Inference based on all introduced UP methods can be parallized, but to a
different degree. To compute the I-posterior using E-Post, we fit a separate
model for each T-posterior draw (or each cluster of draws). This is
embarrassingly parallelizable, as the models are independent so can simply be
run on different cores. However, for each model, a separate MCMC warmup phase
is required, which induces computational overhead.

In contrast, for E-Lik, we fit only one model during the I-step, which loops
over the T-posterior draws when evaluating its likelihood. This however
defines a more complicated posterior, which is substantially slower to sample
from compared to the individual E-Post models. Fitting the single E-Lik model
can be speed up by between-chain parallization when running multiple MCMC
chains (say, one per core), but then we again create overhead due to separate
warmup phases per chain. Alternatively, even when running a single chain, the
E-Lik model can be parallelized via within-chain parallelization (aka
threading), where the likelihood contributions of the T-posterior draws are
evaluated in parallel. An overhead occurs due to variable passing and other
non-parallelized model parts (e.g., the prior density evaluation). This leads
to diminishing returns in terms of the number of cores used for threading. For
more details on threading in Stan, see Bürkner et al., (2022).

E-Log-Lik behaves as E-Lik in terms of parallelizability as they both fit only
a single model during the I-step. That said, in our experiments, E-Log-Lik
models sampled substantially faster than E-Lik models, presumably for two main
reasons. First, E-Log-Lik does not require the use of log-sum-exp in order to
obtain the joint log-likelihood, since we aggregate directly on the log-scale.
This reduces the number of required operations within each MCMC step. Second,
presumably, the geometry of the E-Log-Lik I-posterior is simpler than that of
the E-Lik I-posterior, thus implying a more efficient exploration with MCMC
for the former.

###  2.3 Evaluation of the Two-Step Procedure

Checking the calibration of uncertainty estimates is an important step to
improve the trustworthiness of any inference algorithm. As such, it is an
crucial aspect of the Bayesian workflow (Gelman et al.,, 2013). In our setup,
we are specifically interested in the uncertainty calibration of
$p(\omega_{I}\mid y_{I},\widetilde{\mathcal{M}},u)$, that is our I-posterior
implied by the surrogate. Simulation-based Calibration (SBC) checking (Cook et
al.,, 2006; Talts et al.,, 2018; Modrák et al.,, 2022) is a current gold-
standard approach to validate Bayesian computation, jointly testing the
trinity of the simulator, the probabilistic program, and the posterior
approximation algorithm, e.g., a sampling algorithm such as MCMC.

We extend SBC to the two-step procedure and use the notation of Bürkner et
al., (2023). For any quantile $q\in(0,1)$, let $U_{q}(\omega_{I}\mid
y_{I},\widetilde{\mathcal{M}},u)$ be any uncertainty region given by the
posterior $p(\omega_{I}\mid y_{I},\widetilde{\mathcal{M}},u)$ that depends on
the surrogate $\widetilde{\mathcal{M}}$ and the uncertainty propagation method
$u$ (see Section 2.1 and 2.2). If the generating distribution of the assumed
surrogate $\widetilde{\mathcal{M}}$ is equal to the true data-generating
distribution induced by the simulator the following property holds
simultaneously for all $q\in(0,1)$:

$\displaystyle q=\iint\mathbb{I}[\omega_{I}^{*}\in U_{q}(\omega_{I}\mid
y_{I},\widetilde{\mathcal{M}},u)]\,p(y_{I}\mid\omega_{I}^{*},\mathcal{M})\,p(\omega_{I}^{*})\,dy_{I}\,d\omega_{I}^{*},$
(25)

where we ignore $\sigma_{I}$ for simplicity and $\mathbb{I}[\cdot]$ denotes
the indicator function. This self-consistency property tests the correct
coverage of the uncertainty region conditional on the input for every
quantile. Practically, we check that the prior draws are uniformly distributed
in the surrogate-based samples of $\omega_{I}$. We calculate the ranks as
$r(\omega_{I}^{*},\\{\omega_{I}^{(1)},\dots,\omega_{I}^{(K)}\\})=\sum_{k=1}^{K}\mathbb{I}\left[\omega_{I}^{*}\leq\omega_{I}^{(k)}\right]$
and test them for uniformity using graphical tests as proposed in Säilynoja et
al., (2022). In addition to graphical tests, we calculate the
$\log(\gamma)$-statistic (Säilynoja et al.,, 2022) as a quantitative measure
of uniformity allowing for a faster comparison of calibration (or strength of
miscalibration) between different methods. Checking the calibration of a
specific uncertainty region, e.g., $U_{0.95}(\omega_{I}\mid
y_{I},\widetilde{\mathcal{M}},u)$, corresponds to evaluating Eq. (25) for
$q=0.95$ and is a special case of SBC.

First, as a point of comparison, we explain standard SBC for inference using
the simulator $\mathcal{M}$: (i) Sample a simulation input parameter
$\omega_{I}^{*}\sim p(\omega_{I})$, (ii) conditioned on $\omega_{I}^{*}$
generate measurement output data $y_{I}\sim p(y_{I}\mid\omega_{I}^{*})$, (iii)
using the simulator $\mathcal{M}$ and given the measurements $y_{I}$ draw
samples $\\{\omega_{I}^{(1)},\dots,\omega_{I}^{(K)}\\}$, and (iv) using the
posterior samples we calculate the rank statistics
$r(\omega_{I}^{*},\\{\omega_{I}^{(1)},\dots,\omega_{I}^{(K)}\\})$ from the
posterior $p(\omega_{I}\mid y_{I},\mathcal{M})$. This procedure is repeated
for a chosen number of I-Step trials and uniformity is tested on the stored
ranks, as described above.

However, in the two-step procedure (see Section 2.1), where we additionally
train a surrogate given simulation data (T-Step) and then propagate its
uncertainty to the I-Step, we need to extend SBC as shown in Algorithm 1: We
repeat the T-step multiple times (number of T-Step trials) and, in each
iteration, we simulate a new training dataset $\mathcal{D}_{T}$ used to train
the surrogate $\widetilde{\mathcal{M}}$. Within each such T-Step trial, we
repeat the I-Step for SBC: First, we draw a sample of the simulation
parameters $\omega_{I}^{*}$ and the simulation noise hyperparameter
$\sigma_{I}$ from their respective priors. Next, measurement data $y_{I}$ is
generated using the simulator $\mathcal{M}$. Using the surrogate
$\widetilde{\mathcal{M}}$ and the UP method $u$, we draw posterior samples
$\\{\omega_{I}^{(1)},\dots,\omega_{I}^{(K)}\\}$ from the I-posterior and store
the ranks $r(\omega_{I}^{*},\\{\omega_{I}^{(1)},\dots,\omega_{I}^{(K)}\\})$.
The repetition of the T-Step and I-Step, in addition to covering the whole
input space of $\omega_{I}$, helps to marginalize over both the noise of the
simulator and the noise in the (assumed) measurement process.

With this SBC variant of the two-step procedure, we simultaneously test six
different scenarios, where a failure can indicate one or more of the following
scenarios:

  * •

scenarios also tested in standard SBC

    1. (i)

incorrect implementation of simulator $\mathcal{M}$

    2. (ii)

incorrect implementation of probabilistic program of the surrogate

    3. (iii)

problems with the sampling algorithm

  * •

additional scenarios in proposed SBC for surrogate-based inference

    1. (iv)

inflexible surrogate

    2. (v)

insufficient training of surrogate because of too little simulation training
data

    3. (vi)

inappropriate uncertainty propagation in the surrogate-based inference.

Regarding points (i)-(iii), we assume the simulator to be correct and the
surrogate to be relatively simple (implementation-wise) as well as easy to fit
using MCMC. Concerning points (iv) and (v), we need a sufficiently flexible
surrogate in order to remove the approximation error (i.e., $\sigma_{A}\to 0$)
and an infinite amount of training samples ($N_{T}\to\infty$) in order to
remove the epistemic uncertainty in the posterior
$p(\theta\mid\mathcal{D}_{T})$. However, practically the latter will not be
the case and we expect the calibration to be imperfect. Nonetheless, we will
still be able to compare surrogate models $\widetilde{\mathcal{M}}$ and
different UP methods $u$ by comparing their SBC results, either graphically or
via test statistics.

Algorithm 1 SBC for the Two-Step Procedure

Choose simulator $\mathcal{M}$, surrogate $\widetilde{\mathcal{M}}$, and
uncertainty propagation method $u$

for $m$ in Number T-Step trials do

Draw a training dataset $\mathcal{D}_{T}$ using the simulator $\mathcal{M}$

Fit the surrogate $\widetilde{\mathcal{M}}$ and calculate the posterior
$p(\theta\mid\mathcal{D}_{T})$ (T-Step)

for $n$ in Number I-Step trials do

Draw a prior sample $\omega_{I}^{*}\sim p(\omega_{I})$, $\sigma_{I}\sim
p(\sigma_{I})$

Draw measurements $y_{I}\sim
p(y_{I}\mid\omega_{I}^{*},\sigma_{I},\mathcal{M})$ using simulator
$\mathcal{M}$

Draw posterior samples $\\{\omega_{I}^{(1)},\dots,\omega_{I}^{(K)}\\}\sim
p(\omega_{I}\mid y_{I},\widetilde{\mathcal{M}},u)$ (I-Step)

Store the rank of $\omega_{I}^{*}$ within the set of posterior samples
$\\{\omega_{I}^{(1)},\dots,\omega_{I}^{(K)}\\}$

Perform uniformity test on the stored ranks

##  3 Experiments

We evaluate our surrogate-based Bayesian inference framework in two case
studies: (1) A linear setup, where we propagate only epistemic uncertainty and
(2) a nonlinear setup, in which we propagate both epistemic and aleatoric
uncertainty. All code and material can be found on
GitHub111https://github.com/philippreiser/bayesian-surrogate-uncertainty-
paper.

###  3.1 Case Study 1: Uncertainty Propagation in a Linear Model

In the first case study, we use a linear setup leading to partly analytic
posteriors that allow us to study our framework in a simple, well-understood
scenario.

##### Setup

We consider a simple linear model as simulator:

$\displaystyle y_{T}=\mathcal{M}(\omega_{T},\sigma_{S}=0)=a+b\omega_{T},$ (26)

with simulation input parameter $\omega_{T}$ and two simulation control
parameters $a,b$ set to $a=0.5$, $b=2$, and a simulation noise parameter
$\sigma_{S}$ set to zero, i.e., we use a deterministic simulator. For the
T-Step, we consider $N_{T}=2$ training points
$\mathcal{D}_{T}=\\{\omega_{T}^{(i)},y_{T}^{(i)}\\}_{i=1}^{N_{T}}$ (chosen
according to Table 2), where $\omega_{T}^{(i)}$ denotes the simulation input
and $y_{T}^{(i)}=\mathcal{M}(\omega_{T}^{(i)})$ the corresponding simulation
output. We denote by
$\Omega_{T}=\bigl{[}\begin{smallmatrix}1&\omega_{T}^{(1)}\\\
1&\omega_{T}^{(2)}\end{smallmatrix}\bigr{]}$ the design matrix of all input
parameters and by $y_{T}$ the vector of all simulation outputs.

We set the surrogate model to a linear model as well, that is, we use the same
model class as for the simulator:

$\displaystyle\widetilde{\mathcal{M}}(\omega_{T},c)=c_{1}+c_{2}\omega_{T},$
(27)

where the intercept $c_{1}$ and slope $c_{2}$ form the surrogate approximation
parameters $c=[c_{1},c_{2}]^{\top}$. For the T-Step, we only consider the
surrogate approximation parameters $c$ as trainable surrogate parameters. Even
though the surrogate can match the simulator perfectly, we fix the surrogate
approximation error hyperparameter $\sigma_{A}$ to values greater than zero.
This allows us to control the width of the T-posterior and thus the epistemic
uncertainty, as shown below. As prior on the surrogate approximation
parameters we choose a bivariate normal distribution:

$\displaystyle p(c)=\mathcal{N}(c\mid\mu_{T0},\Sigma_{T0}),$ (28)

with mean $\mu_{T0}=0$ and a covariance matrix $\Sigma_{T0}=\sigma_{T0}^{2}I$
where $I$ is the identity matrix. We set the T-likelihood to a normal
distribution as well:

$\displaystyle p(y_{T}\mid c)=\prod_{i=1}^{N_{T}}\mathcal{N}(y_{T}^{(i)}\mid
c_{1}+c_{2}\omega_{T}^{(i)},\sigma_{A}^{2}).$ (29)

To generate data for the I-Step, a single measurement
$y_{I}=\mathcal{M}(\omega_{I}^{*})$ is obtained by inputting a true input
parameter $\omega_{I}^{*}$ to the deterministic simulator $\mathcal{M}$. In
the following, we will also use the augmented vector
$\hat{\omega}_{I}=[1,\omega_{I}]^{T}$ to simplify notation. Within the I-step,
we set a normal prior with mean $\mu_{I0}$ and variance $\sigma_{I0}^{2}$ on
$\omega_{I}$:

$\displaystyle
p(\omega_{I})=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2}).$ (30)

The I-likelihood $p(y_{I}\mid\omega_{I})$ is assumed to be normal with fixed
variance $\sigma_{I}^{2}$:

$\displaystyle p(y_{I}\mid\omega_{I},c)=\mathcal{N}(y_{I}\mid
c_{1}+c_{2}\omega_{I},\sigma_{I}^{2})$ (31)

##### Deriving the posteriors

In the following, we perform the T-Step (see Section 2.1.1) and I-Step (see
Section 2.2) for the linear setup. In the T-Step, we calculate the T-posterior
of the surrogate approximation parameters $c$, which contains epistemic
uncertainty. Then we derive the four different I-step methods when propagating
only the epistemic uncertainty from the T-posterior. The aleatoric uncertainty
is zero because, in our setup, the simulator is within the approximation space
of the surrogate model.

##### T-Step.

We calculate the T-posterior for the surrogate approximation parameters $c$
given the training data $\mathcal{D}_{T}$ from the simulator by using the
conjugate prior relation for a normal-normal model (Murphy,, 2007, 2012)
leading to a normal posterior:

$\displaystyle p(c\mid\mathcal{D}_{T})$
$\displaystyle=\mathcal{N}(c\mid\mu_{T1},\Sigma_{T1}),$ (32)

with

$\displaystyle\Sigma_{T1}$
$\displaystyle=(\Sigma_{T0}^{-1}+\sigma_{A}^{-2}\Omega_{T}^{\top}\Omega_{T})^{-1},$
(33) $\displaystyle\mu_{T1}$
$\displaystyle=\Sigma_{T1}(\Sigma_{T0}^{-1}\mu_{T0}+\sigma_{A}^{-2}\Omega_{T}y_{T}).$
(34)

We see that by fixing $\sigma_{A}$ we can control the width of the T-posterior
and hence the epistemic uncertainty of the surrogate parameters, even if we do
not propagate $\sigma_{A}$ itself.

##### I-Step.

Below, we derive the I-posteriors of all four uncertainty propagation
procedures. We compute the mean of the T-posterior: $\bar{c}=\mu_{T1}$ and use
it to calculate the Point I-posterior as follows:

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{Point})$ $\displaystyle\propto
p(\omega_{I})p(y_{I}\mid\omega_{I},\bar{c})$ (35)
$\displaystyle=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\mathcal{N}(y_{I}\mid\mu_{T1}^{(1)}+\mu_{T1}^{(2)}\omega_{I},\sigma_{I}^{2})$
(36)
$\displaystyle\propto\mathcal{N}(\omega_{I}\mid\mu_{I1},\sigma_{I1}^{2}),$
(37)

with

$\displaystyle\sigma_{I1}^{2}$
$\displaystyle=(\sigma_{I0}^{-2}+\sigma_{I}^{-2}\mu_{T1}^{(2)}\mu_{T1}^{(2)})^{-1},$
(38) $\displaystyle\mu_{I1}$
$\displaystyle=\sigma_{I1}^{2}(\sigma_{I0}^{-2}\mu_{I0}+\sigma_{I}^{-2}\mu_{T1}^{(2)}(y_{I}-\mu_{T1}^{(1)})).$
(39)

Again, this works due to the normal-normal conjugacy. We see that the
I-posterior variance $\sigma_{I1}^{2}$ does not depend on the T-posterior
covariance matrix $\Sigma_{T1}$, i.e. the uncertainty from the surrogate
training step is neglected.

The E-Log-Lik I-posterior is given by

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Log-Lik})\propto
p(\omega_{I})\exp\left\\{\int\log(p(y_{I}\mid\omega_{I},c))p(c\mid\mathcal{D}_{T})dc\right\\}$
(40)
$\displaystyle=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\exp\left\\{\int\log(\mathcal{N}(y_{I}\mid\hat{\omega}_{I}^{\top}c,\sigma_{I}^{2}))\mathcal{N}(c\mid\mu_{T1},\Sigma_{T1})dc\right\\}$
(41)
$\displaystyle\propto\mathcal{N}(\omega_{I}\mid\mu_{I1},\sigma_{I1}^{2}),$
(42)

with

$\displaystyle\sigma_{I1}^{2}$
$\displaystyle=(\sigma_{I0}^{-2}+\sigma_{I}^{-2}(\mu_{T1}^{(2)}\mu_{T1}^{(2)}+\Sigma_{T1}^{(2,2)}))^{-1},$
(43) $\displaystyle\mu_{I1}$
$\displaystyle=\sigma_{I1}^{2}(\sigma_{I0}^{-2}\mu_{I0}+\sigma_{I}^{-2}(\mu_{T1}^{(2)}y_{I}-\Sigma_{T1}^{(1,2)}-\mu_{T1}^{(1)}\mu_{T1}^{(2)}))$
(44)

Accordingly, it is also a normal distribution, but a different one from the
Point I-posterior. The detailed derivation of the E-Log-Lik is given in
Appendix B.3.1. If we look at the variance $\sigma_{I1}^{2}$ we see that
E-Log-Lik produces counter-intuitive results, as $\Sigma_{T1}^{(2,2)}$ and
$\sigma_{I1}^{2}$ are reciprocally related. That is, as the surrogate model
gets more uncertain in the T-step, the I-posterior gets more certain.

The E-Lik I-posterior is computed as

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Lik})$ $\displaystyle\propto
p(\omega_{I})\int p(y_{I}\mid\omega_{I},c)p(c\mid\mathcal{D}_{T})dc$ (45)
$\displaystyle=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\int\mathcal{N}(y_{I}\mid\hat{\omega}_{I}^{\top}c,\sigma_{I}^{2})\mathcal{N}(c\mid\mu_{T1},\Sigma_{T1})dc$
(46)
$\displaystyle=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\mathcal{N}(y_{I}\mid\hat{\omega}_{I}^{\top}\mu_{T1},\hat{\omega}_{I}^{\top}\Sigma_{T1}\hat{\omega}_{I}+\sigma_{I}^{2}).$
(47)

Here, we cannot apply the normal-normal conjugacy, because $\hat{\omega}_{I}$
is present in both the mean and the variance in the second term of the product
of the two normals. Instead, the resulting distribution is non-analytic and we
have to use numerical integration to calculate the normalization constant.
Looking at the derived quantity, we see that the variance of the term
resulting from the marginalization over the surrogate approximation parameters
$c$ increases with the variance of the T-posterior $\Sigma_{T1}$.

Finally, we calculate the I-posterior using the E-Post approach:

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Post})$
$\displaystyle=p(\omega_{I})\int\frac{p(y_{I}\mid\omega_{I},c)p(c\mid\mathcal{D}_{T})}{\int
p(y_{I}\mid\omega_{I},c)p(\omega_{I})d\omega_{I}}dc$ (48)
$\displaystyle=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\int\frac{\mathcal{N}(y_{I}\mid\hat{\omega}_{I}^{\top}c,\sigma_{I}^{2})\mathcal{N}(c\mid\mu_{T1},\Sigma_{T1})}{\mathcal{N}(y_{I}\mid
c_{1}+c_{2}\mu_{I0},c^{\top}c\sigma_{I0}^{2}+\sigma_{I}^{2})}dc.$ (49)

Again, this integral is non-analytic and numerical integration is required to
calculate the E-Post I-posterior. Nevertheless, we can expect similar behavior
to E-Lik, since we also marginalize over the T-posterior and only normalize
differently.

##### Results

Table 2: Parameters with realized values for Case study 1. Parameters | T-Step | I-Step
---|---|---
$\begin{bmatrix}\omega_{T}^{(1)}\\\ \omega_{T}^{(2)}\end{bmatrix}$ | $\begin{bmatrix}a\\\ b\end{bmatrix}$ | $\mu_{T0}$ | $\sigma_{T0}$ | $\sigma_{A}$ | $\omega_{I}^{*}$ | $\mu_{I0}$ | $\sigma_{I0}$ | $\sigma_{I}$
Values | $\begin{bmatrix}-0.9\\\ -0.3\end{bmatrix}$ | $\begin{bmatrix}0.5\\\ 2\end{bmatrix}$ | $\begin{bmatrix}0\\\ 0\end{bmatrix}$ | $10$ | $\\{0.1,0.5,1\\}$ | -0.5 | 0 | 1 | 0.1

In Table 2, we provide an overview over the chosen training data, measurement
data, and hyperparameters for the T and I-Step. We specifically vary the
surrogate approximation error parameter $\sigma_{A}=\\{0.1,0.5,1\\}$ to
control the epistemic uncertainty during the surrogate training, as explained
above. To compare the four I-Steps, we calculate their I-posteriors under the
given scenarios. The results are illustrated in Fig. 3. We see that Point
produces results that are independent of the standard deviation of the
T-likelihood. This is expected since it does not propagate the T-epistemic
uncertainty at all. Intuitively, for the other three methods, the I-posteriors
should become more uncertain as we increase the uncertainty in the T-step,
since our surrogate model gets less trustworthy. This is indeed what happens
for E-Lik and E-Post, which also produce very similar but not identical
results. In contrast, E-Log-Lik behaves counter-intuitively since its
I-posterior becomes more certain as the T-posterior becomes more uncertain, a
behavior that we also see clearly from Eq. (43), as noted above.

Figure 3: I-posterior densities for the linear surrogate with normal
priors/likelihoods in Case study 1. We use the data and parameters as
specified in Table 2. We use four different UP methods to compute the
I-posterior while the surrogate approximation error
$\sigma_{A}=\\{0.1,0.5,1\\}$ is varied.

###  3.2 Case study 2: Uncertainty Propagation in a Logistic Model

The second case study examines a nonlinear problem where we have to rely on
MCMC for the posterior approximations, because analytic posteriors are
unavailable. We choose the one-dimensional logistic function

$\displaystyle y=\mathcal{M}(\omega)=\frac{2}{1+\exp(-10\omega)}-1$ (50)

as the (true) simulator since it is invertible and smooth everywhere. We
consider two different surrogate models: The first is a parameterized
generalization of the simulator and the second is a polynomial chaos expansion
(PCE) surrogate. We will discuss these two cases separately below.

####  3.2.1 Logistic Surrogate Model

##### Setup.

As surrogate, we consider

$\displaystyle\widetilde{\mathcal{M}}(\omega;c)=\frac{\alpha}{1+\exp(-\beta(\omega-\gamma))}-\delta,$
(51)

with surrogate approximation parameters $c=(\alpha,\beta,\gamma,\delta)$.
Here, the true simulator is contained in the set of surrogate models (for
$\alpha=2,\beta=10,\gamma=0,\delta=-1$).

For the T-Step, we generate the training set by setting the input points
$\omega_{T}^{(i)}\in[-1,1]$ to the first $N_{T}$ points of a slightly modified
one-dimensional Halton sequence (Halton,, 1960), which starts with the
boundary points ($-1$ and $1$) and then progresses as the standard Halton
sequence with center point 0. The simulated output responses $y_{T}^{(i)}$ for
$i\in\\{1,\dots,N_{T}\\}$ are sampled from a normal distribution with mean
equal to the evaluated logistic simulator $\mathcal{M}$ at the input points
and standard deviation $\sigma_{S}$, i.e.
$y_{T}^{(i)}\sim\mathcal{N}(\mathcal{M}(\omega_{T}^{(i)}),\sigma_{S}^{2})$. To
avoid sampling issues during the training step, we induce a small simulation
noise $\sigma_{S}=0.01$, which is not explicitly modelled.

For the surrogate parameters $c$, we specify normal priors with means around
the true values and a standard deviation of 1 (except for $\beta$, where we
set the standard deviation to 10). For the T-likelihood, we consider a normal
distribution:

$\displaystyle p(y_{T}\mid
c)=\mathcal{N}(y_{T}\mid\widetilde{\mathcal{M}}(\omega_{T},c),\sigma_{A}^{2})$
(52)

with the surrogate approximation hyperparameter $\sigma_{A}$. On each training
data set, we fit the parameters of our surrogate using Markov chain Monte
Carlo (MCMC). Specifically, we use the no-U-turn-sampler (NUTS) (Hoffman and
Gelman,, 2014), an adaptive form of Hamiltonian Monte Carlo (HMC) (Neal,,
2011) which is a gradient based MCMC sampler via the probabilistic programming
language Stan (Stan Development Team,, 2022). We use four chains, each running
for 1250 iterations (1000 warmup and 250 post-warmup iterations), resulting in
a total of 1000 T-posterior draws. We checked convergence using standard
convergence checks (Vehtari et al.,, 2021).

For the I-Step, we generate $N_{I}=5$ random measurements
$y_{I}\sim\mathcal{N}(\mathcal{M}(\omega_{I}^{*}),\sigma_{I}^{2})$, based on
the simulator, true input parameters $\omega_{I}^{*}$, and the measurement
error $\sigma_{I}=0.01$. As priors we set
$\omega_{I}\sim\mathcal{N}_{[-1,1]}(0,0.5^{2})$ (with truncation bounds
$[-1,1]$) and $\sigma_{I}\sim\text{uniform}(0,0.05)$.

In contrast to Case Study 1, we propagate the T-posterior through samples and
hence use the MC approximation (see Section 2.2) for each of the four methods
(Point, E-Lik, E-Log-Lik, E-Post). We sample from the I-posterior of
$\omega_{I}$ and $\sigma_{I}$ using MCMC with four chains, each running for
5000 iterations (1000 warmup and 4000 post-warmup iterations), resulting in a
total of 16000 I-posterior draws.

Similar to Case Study 1, we propagate only the epistemic uncertainty that is
encoded in the T-posterior of the surrogate parameters $c$. For this purpose,
we either utilize all T-posterior draws or employ KMeans clustering (see
Section 2.2.5) with $L=25$ clusters. As the simulator is contained in the
class of surrogates, no approximation error is present ($\sigma_{A}=0$) and
hence, there is no aleatoric uncertainty to propagate.

##### Posterior Distributions.

We set $\sigma_{S}=0.01$ and vary the number of training points
$N_{T}=\\{5,7,10\\}$ to perform the T-Step. In Fig. 4, the left column shows
the mean of the T-posterior predictive with the 95%-credible interval (CI) to
display its epistemic uncertainty. As expected, increasing $N_{T}$ leads to
smaller T-epistemic uncertainty.

We compare the I-posteriors resulting from the four different methods using
the measurements resulting from exemplary underlying true inputs
$\omega_{I}^{*}\in\\{-0.05,0.1,0.3\\}$ in the three right columns in Fig. 4.
The Point method yields I-posteriors with constant width regardless of
variation in $N_{T}$. The I-posteriors of E-Lik and E-Post show similar
behavior and become more uncertain as T-epistemic uncertainty increases. The
E-Log-Lik follows a similar trend, but its I-posteriors have qualitatively
different shapes and are narrower. For $N_{T}=5$, only E-Lik and E-Post have
substantial I-posterior probability mass on the true inputs $\omega_{I}^{*}$.
Notably, all uncertainty propagation methods converge to the same results as
the epistemic uncertainty from the T-Step decreases.

Figure 4: Selected results for two-step procedure with the logistic surrogate
in Case study 2. Left: For $N_{T}=\\{5,7,10\\}$ the training data set
$\mathcal{D}_{T}$ (black dots) and the mean of the T-posterior predictive
distribution (red lines) is shown. Right: For each underlying true input
$\omega_{I}^{*}\in\\{-0.05,0.1,0.3\\}$ (black vertical lines) we depict the
I-posterior distributions for each Point, E-Lik, E-Post, and E-Log-Lik
(colored lines).

##### Calibration.

To check if the methods for estimating the I-posteriors are calibrated, we use
SBC checking (Talts et al.,, 2018; Modrák et al.,, 2022) with the SBC R
package (Kim et al.,, 2023). Concretely, we perform the adapted SBC procedure
for surrogate-based inference as detailed in Section 2.3. For the I-Step
trials we simulated the true values of the inputs $\omega_{I}^{*}$ and
measurement error $\sigma_{I}$ from the above chosen priors. We perform 20
I-step trials within each of 10 T-Step trials resulting in a total of 200 SBC-
trials.

To evaluate the calibration of the four methods graphically, we show the ECDF-
difference plots (Säilynoja et al.,, 2022) in the upper part of Fig. 5. We
choose two scenarios: high and low T-epistemic uncertainty, represented by the
number of training points $N_{T}=\\{5,10\\}$. In this graphical test,
calibration is achieved if the black line lies within the blue region, i.e.
the 95%-confidence envelops. We observe that, while high T-epistemic
uncertainty is present, E-Post and E-Lik are almost calibrated, whereas Point
and E-Log-Lik are overconfident. For $N_{T}=8$ all methods show good
calibration.

As an additional continuous calibration metric, we calculate the
$\log(\gamma)$-statistic (Modrák et al.,, 2022) from our SBC results. We vary
$N_{T}=\\{5,6,7,8,9,10\\}$ to control for the amount of to-be-propagated
T-epistemic uncertainty. The center part of Fig. 5 shows the corresponding
results. As a general trend, we observe that E-Lik and E-Post are similarly
calibrated. While first ($N_{T}=5$) slightly miscalibrated, for $N_{T}>5$
proper calibration is achieved. In contrast, E-Log-Lik and Point are both
miscalibrated for $N_{T}<7$, but with more training data, as the T-epistemic
uncertainty diminishes, they also become well calibrated.

In the bottom of Fig. 5 we depict the sharpness (Gneiting et al.,, 2007;
Bürkner et al.,, 2023) of the I-posteriors, here measured by the width of the
90 % CI of $\omega_{I}$. For similarly calibrated methods, we say that the
method with the smaller CI is sharper. We observe that Point and E-Log-Lik
produce overall sharper results, but given their bad calibration, its clear
that these two methods are just overconfident.

Figure 5: Calibration and sharpness of the I-posteriors using logistic
surrogate in Case study 2. Top: ECDF difference plots for the I-posterior
distributions of $\omega_{I}$ resulting from the four different methods. The
blue areas in the ECDF difference plots indicate 95%-confidence envelopes and
the black lines indicate the empirical cumulative distribution function (ECDF)
for two different number of simulation points $N_{T}=\\{5,10\\}$. Center: log-
gamma-statistics of SBC with calibration threshold depicted as black
horizontal line. Bottom: sharpness (90% CI) of I-posterior for four different
I-Steps (colored dots/lines) for $N_{T}=\\{5,6,7,8,9,10\\}$.

####  3.2.2 Polynomial Surrogate Model

##### Setup.

We now make the task substantially harder by not including the true simulator
in the set of considered surrogate models. For this purpose, we use a
polynomial chaos expansion (PCE) model as surrogate (Wiener,, 1938; Sudret,,
2008; Oladyshkin and Nowak,, 2012; Bürkner et al.,, 2023):

$\displaystyle\widetilde{\mathcal{M}}(\omega;c)=\sum_{i=0}^{d}c_{i}\psi_{i}(\omega),$
(53)

with the vector of surrogate coefficients $c=(c_{0},...,c_{d})$, the maximum
degree of polynomials $d$, and Legendre polynomials $\psi_{i}(\omega)$ (see
Sudret, (2008) for a detailed definition). We consider wide, independent
normal priors for all surrogate coefficients: $c_{i}\sim\mathcal{N}(0,5)$. In
the following, we fix the maximum polynomial degree to $d=5$. By doing so, we
create a scenario in which our surrogate is unable to fit the underlying true
model appropriately such that an approximation error $e_{A}$ is induced. This
creates a scenario in which it is important to propagate both T-epistemic and
T-aleatoric uncertainty.

Figure 6: Selected results for two-step procedure with PCE surrogate in Case
study 2. Left: For $N_{T}=\\{10,20\\}$ the training data set $\mathcal{D}_{T}$
(black dots), the T-posterior predictive distribution and the mean of the
T-posterior predictive distribution (dark red and red lines) is shown. Right:
For the true input $\omega_{I}^{*}=-0.05$ (black vertical line) we depict the
I-posterior distributions for each Point, E-Lik, E-Post, and E-Log-Lik
(colored lines). In the center column we only propagate T-epistemic
uncertainty via $\theta=c$ and in the right column we propagate both
T-epistemic and T-aleatoric uncertainty via $\theta=(c,\sigma_{A})$.

##### Posterior Distributions.

We perform the T- and I-Step, as previously described in Section 3.2.1. During
the T-Step, we vary the number of training points $N_{T}=\\{10,20\\}$ to
control for the T-epistemic uncertainty. For the I-Step, we set the true input
parameters to $\omega_{I}^{*}=-0.05$ and set $\sigma_{S}=0.01$ to avoid
sampling issues. We show the result of the T-Step in the left column of Fig.
6, where we present both the T-posterior predictive distribution (T-aleatoric
and T-epistemic uncertainty) and the mean of the T-posterior predictive
distribution (T-epistemic uncertainty only). The T-epistemic uncertainty
becomes smaller with increasing number of training data points $N_{T}$, but
the T-aleatoric uncertainty stays approximately constant.

In the middle column, we show the results of the I-posterior for the four
methods when propagating only the T-epistemic uncertainty by considering only
the T-posterior draws of the surrogate approximation parameters, i.e.
$\theta=c$. For $N_{T}=10$ the Point I-posterior and the E-Log-Lik I-posterior
are narrow despite high T-epistemic uncertainty. In contrast, E-Lik and E-Post
produce wider I-posteriors that are similar to each other. As the T-epistemic
uncertainty reduces, all methods produce I-posterior distributions that
converge towards a similar distribution.

In the right column, we show the I-posteriors when propagating both
T-epistemic and T-aleatoric uncertainty by also propagating posterior draws of
the surrogate approximation error: $\theta=(c,\sigma_{A})$ (see Section
2.1.2). In general, all I-posteriors now tend to be wider than for $\theta=c$
and are more similar to each other. However, in the presence of substantial
T-epistemic uncertainty (as for $N_{T}=10$), E-Post produces the widest
I-posterior, followed by E-Lik, Point, and E-Log-Lik. In Appendix C.3 we
provide further results for different true parameter values.

##### Calibration.

We performed SBC for our PCE surrogate setup and show the results in Fig. 7.
The top part shows the ECDF-difference plots for $N_{T}=10$ for two cases:
$\theta=c$ and $\theta=(c,\sigma_{A})$. When only T-epistemic uncertainty is
propagated, we see that only E-Lik and E-Post produce calibrated results,
while E-Log-Lik and Point are overconfident. When T-aleatoric uncertainty is
propagated as well, all calibrations improve, while E-Post still produces the
best calibration results closely followed by E-Lik. In the bottom of Fig. 7 we
compare the calibration via the $\log(\gamma)$-statistic of SBC under varying
$N_{T}=\\{10,20,30,40,50\\}$ and confirm the observed pattern described above.
As we use a surrogate model which was chosen on purpose to be inflexible and
$\sigma_{A}$ is modelled as constant over $\omega_{I}$, we cannot achieve
proper calibration with either method as epistemic uncertainty reduces.

Figure 7: Calibration and sharpness of the I-posteriors using PCE surrogate in
Case study 2. Top: ECDF difference plots for the I-posterior distributions of
$\omega_{I}$ resulting from the four different methods. We set the number
training points $N_{T}=10$ and show the results for T-epistemic uncertainty
propagation ($\theta=c$) and T-epistemic and T-aleatoric uncertainty
propagation ($\theta=(c,\sigma_{A})$). Center: log-gamma-statistics of SBC.
Bottom: sharpness (90 % CI) of I-posterior for four different I-Steps (colored
dots/lines) for $N_{T}=\\{10,20,30,40,50\\}$

##  4 Conclusion

We introduced a general two-step procedure for surrogate-based Bayesian
inference. Within our approach, we propagate all relevant surrogate
uncertainties (both aleatoric and epistemic) from the surrogate training to
the real-data inference step, thereby producing fully uncertainty aware
inference on the real data. The uncertainty propagation methods developed
within our framework are in principle agnostic to the chosen surrogate, but
require the ability to sample from its training-step posterior given data
generated from the simulator. To evaluate the uncertainty calibration of the
resulting inference, we proposed an extension of simulation-based calibration
suitable for our two-step surrogate approach.

As we demonstrate in our case studies, even in seemingly simple setups,
complex behavior occurs in terms of posterior shape and calibration when
propagating surrogate uncertainty. In particular two uncertainty propagation
methods (E-Lik & E-Post) showed substantial improvements in uncertainty
calibration compared to traditional surrogate-based inference that is
uncertainty-unaware. What is more, our results show the importance of
propagating the complete surrogate uncertainty (aleatoric and epistemic)
instead of propagating only parts of it. Intuitively, one might expect that
there is only one ”correct” uncertainty propagation method within the bounds
of probability theory. However, as we demonstrate, the two advocated methods
for uncertainty propagation (E-Lik & E-Post) produce non-equivalent inference
even in simple cases despite being equally justified by probability theory.
That said, at least in our case studies, the produced inference was very
similar. They differ in the computational requirements and ease of
parallelization though (see Section 2.2.6), such that either or the other may
be preferable depending the context and available resources.

##### Future Work.

Our surrogate-based Bayesian inference approach is agnostic to the input-
parameter dimensionality and its scaling is unaffected by said dimensionality
(but only by the number of draws propagated from the training to the inference
step). However, in our case studies, we focused on simulators with single-
dimensional input parameters in order to simplify the presentation and
establish a better intuition about the overall approach. To study the
applicability of our methods to higher-dimensional, real-world challenges,
biological systems requiring accurate, uncertainty-aware inference (Mitra and
Hlavacek,, 2019), will offer an interesting class of problems for future
research.

In our case studies, we modelled the surrogate approximation error to be
constant, thus assuming the aleatoric uncertainty of the surrogate to be
independent of the input parameters. However, this assumption is likely
unjustified in practice, if there is substantial approximation error due to
the surrogate inflexibility. Hence, one could improve the modeling of the
surrogate approximation error by conditioning it on the input parameters, for
example, as suggested in Kohlhaas et al., (2023). This will induce additional
parameters to model the approximation error, which are then seamlessly
propagated through our two-step procedure.

Our current implementation of uncertainty quantification relies on MCMC in
both training and inference step. This may become computationally expensive,
sometimes prohibitively so, as the number of model parameters grows (Izmailov
et al.,, 2021). Hence, for a surrogate model with hundreds of thousands of
parameters, full Bayesian UQ will become intractable and approximations are
needed. Approaches like variational inference (Hinton and van Camp,, 1993;
Graves,, 2011) and partial UQ on the last neural network layer (Kristiadi et
al.,, 2020; Fiedler and Lucia,, 2023) could be sensible alternatives. Our
uncertainty propagation methods only require the ability to draw samples from
an (approximate) posterior in the training step and thus can be combined with
any inference algorithm that produces samples. This flexibility ensures that
our approach is still applicable in scenarios where a full Bayesian UQ is
infeasible, but the implication on inference validity, in particular
uncertainty calibration, need to be further studied.

## Acknowledgments

Partially funded by Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy EXC 2075 – 390740016 and DFG
Project 500663361. We acknowledge the support by the Stuttgart Center for
Simulation Science.

## References

  * Alden et al., (2020) Alden, K., Cosgrove, J., Coles, M., and Timmis, J. (2020).  Using emulation to engineer and understand simulations of biological systems.  IEEE/ACM Transactions on Computational Biology and Bioinformatics, 17(1):302–315.
  * Blomstedt et al., (2019) Blomstedt, P., Mesquita, D., Lintusaari, J., Sivula, T., Corander, J., and Kaski, S. (2019).  Meta-analysis of bayesian analyses.
  * Brandstetter et al., (2023) Brandstetter, J., van den Berg, R., Welling, M., and Gupta, J. K. (2023).  Clifford neural layers for PDE modeling.  In The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.net.
  * Breiman, (2004) Breiman, L. (2004).  Bagging predictors.  Machine Learning, 24:123–140.
  * Bühlmann, (2014) Bühlmann, P. (2014).  Discussion of Big Bayes Stories and BayesBag.  Statistical Science, 29(1):91 – 94.
  * Bürkner et al., (2023) Bürkner, P.-C., Kröker, I., Oladyshkin, S., and Nowak, W. (2023).  A fully bayesian sparse polynomial chaos expansion approach with joint priors on the coefficients and global selection of terms.  J. Comput. Phys., 488:112210.
  * Bürkner et al., (2023) Bürkner, P.-C., Scholz, M., and Radev, S. T. (2023).  Some models are useful, but how do we know which ones? Towards a unified Bayesian model taxonomy.  Statistics Surveys, 17:216 – 310.
  * Bürkner, (2017) Bürkner, P.-C. (2017).  brms: An R package for Bayesian multilevel models using Stan.  Journal of Statistical Software, 80(1):1–28.
  * Bürkner, (2018) Bürkner, P.-C. (2018).  Advanced Bayesian Multilevel Modeling with the R Package brms.  The R Journal, 10(1):395–411.
  * Bürkner et al., (2022) Bürkner, P.-C., Gabry, J., Weber, S., Johnson, A., Modrak, M., Badr, H. S., Weber, F., Ben-Shachar, M. S., and Rabel, H. (2022).  Running brms models with within-chain parallelization.  URL https://CRAN.R-Project.org/package=brms. Vignette included in R package brms, version 2.17.0.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of Statistical Software, 76(1):1–32.
  * Cook et al., (2006) Cook, S. R., Gelman, A., and Rubin, D. B. (2006).  Validation of software for bayesian models using posterior quantiles.  Journal of Computational and Graphical Statistics, 15(3):675–692.
  * Dempster et al., (1977) Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977).  Maximum likelihood from incomplete data via the em algorithm.  Journal of the Royal Statistical Society. Series B (Methodological), 39(1):1–38.
  * Douady et al., (2003) Douady, C. J., Delsuc, F., Boucher, Y., Doolittle, W. F., and Douzery, E. J. P. (2003).  Comparison of Bayesian and Maximum Likelihood Bootstrap Measures of Phylogenetic Reliability.  Molecular Biology and Evolution, 20(2):248–254.
  * Efron, (1979) Efron, B. (1979).  Bootstrap Methods: Another Look at the Jackknife.  The Annals of Statistics, 7(1):1 – 26.
  * Fiedler and Lucia, (2023) Fiedler, F. and Lucia, S. (2023).  Improved uncertainty quantification for neural networks with bayesian last layer.  IEEE Access, 11:123149–123160.
  * Gelman et al., (2013) Gelman, A., Carlin, J., Stern, H., Dunson, D., Vehtari, A., and Rubin, D. (2013).  Bayesian Data Analysis, Third Edition.  Chapman & Hall/CRC Texts in Statistical Science. Taylor & Francis.
  * Geyer, (1991) Geyer, C. J. (1991).  Markov chain monte carlo maximum likelihood.
  * Gneiting et al., (2007) Gneiting, T., Balabdaoui, F., and Raftery, A. E. (2007).  Probabilistic Forecasts, Calibration and Sharpness.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 69(2):243–268.
  * Goodfellow et al., (2016) Goodfellow, I., Bengio, Y., and Courville, A. (2016).  Deep Learning.  MIT Press.  http://www.deeplearningbook.org.
  * Gramacy, (2020) Gramacy, R. B. (2020).  Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences.  Chapman Hall/CRC, Boca Raton, Florida.  http://bobby.gramacy.com/surrogates/.
  * Graves, (2011) Graves, A. (2011).  Practical variational inference for neural networks.  In Shawe-Taylor, J., Zemel, R., Bartlett, P., Pereira, F., and Weinberger, K., editors, Advances in Neural Information Processing Systems, volume 24. Curran Associates, Inc.
  * Gruber et al., (2023) Gruber, C., Schenk, P. O., Schierholz, M., Kreuter, F., and Kauermann, G. (2023).  Sources of uncertainty in machine learning – a statisticians’ view.
  * Halton, (1960) Halton, J. H. (1960).  On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals.  Numer. Math., 2(1):84–90.
  * Hinton and van Camp, (1993) Hinton, G. E. and van Camp, D. (1993).  Keeping the neural networks simple by minimizing the description length of the weights.  In Proceedings of the Sixth Annual Conference on Computational Learning Theory, COLT ’93, page 5–13, New York, NY, USA. Association for Computing Machinery.
  * Hoffman and Gelman, (2014) Hoffman, M. D. and Gelman, A. (2014).  The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo.  Journal of Machine Learning Research, 15(47):1593–1623.
  * Huggins and Miller, (2020) Huggins, J. H. and Miller, J. W. (2020).  Robust inference and model criticism using bagged posteriors.
  * Hüllermeier and Waegeman, (2021) Hüllermeier, E. and Waegeman, W. (2021).  Aleatoric and epistemic uncertainty in machine learning: an introduction to concepts and methods.  Machine Learning, 110(3):457–506.
  * Izmailov et al., (2021) Izmailov, P., Vikram, S., Hoffman, M. D., and Wilson, A. G. (2021).  What are bayesian neural network posteriors really like?  In Meila, M. and Zhang, T., editors, Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event, volume 139 of Proceedings of Machine Learning Research, pages 4629–4640. PMLR.
  * Kallioinen et al., (2022) Kallioinen, N., Paananen, T., Bürkner, P.-C., and Vehtari, A. (2022).  Detecting and diagnosing prior and likelihood sensitivity with power-scaling.
  * Kim et al., (2023) Kim, S., Moon, H., Modrák, M., and Säilynoja, T. (2023).  SBC: Simulation Based Calibration for rstan/cmdstanr models.  https://hyunjimoon.github.io/SBC/, https://github.com/hyunjimoon/SBC/.
  * Kohlhaas et al., (2023) Kohlhaas, R., Kröker, I., Oladyshkin, S., and Nowak, W. (2023).  Gaussian active learning on multi-resolution arbitrary polynomial chaos emulator: concept for bias correction, assessment of surrogate reliability and its application to the carbon dioxide benchmark.  Computational Geosciences, 27:1–21.
  * Kristiadi et al., (2020) Kristiadi, A., Hein, M., and Hennig, P. (2020).  Being bayesian, even just a bit, fixes overconfidence in relu networks.  In Proceedings of the 37th International Conference on Machine Learning, ICML’20. JMLR.org.
  * Kuehnert et al., (2022) Kuehnert, J., McGlynn, D., Remy, S. L., Walcott-Bryant, A., and Jones, A. (2022).  Surrogate Ensemble Forecasting for Dynamic Climate Impact Models.  arXiv:2204.05795 [physics].
  * Laloy et al., (2013) Laloy, E., Rogiers, B., Vrugt, J., Mallants, D., and Jacques, D. (2013).  Efficient posterior exploration of a high-dimensional groundwater model from two-stage mcmc simulation and polynomial chaos expansion.  Water Resources Research, 49.
  * Lavin et al., (2021) Lavin, A., Zenil, H., Paige, B., Krakauer, D., Gottschlich, J., Mattson, T., Anandkumar, A., Choudry, S., Rocki, K., Baydin, A. G., Prunkl, C., Isayev, O., Peterson, E., McMahon, P. L., Macke, J. H., Cranmer, K., Zhang, J., Wainwright, H. M., Hanuka, A., Veloso, M., Assefa, S., Zheng, S., and Pfeffer, A. (2021).  Simulation intelligence: Towards a new generation of scientific methods.  CoRR, abs/2112.03235.
  * Li and Marzouk, (2014) Li, J. and Marzouk, Y. M. (2014).  Adaptive construction of surrogates for the bayesian solution of inverse problems.  SIAM Journal on Scientific Computing, 36(3):A1163–A1186.
  * Li et al., (2021) Li, Z., Kovachki, N. B., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A. M., and Anandkumar, A. (2021).  Fourier neural operator for parametric partial differential equations.  In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net.
  * Little and Rubin, (2019) Little, R. J. and Rubin, D. B. (2019).  Statistical analysis with missing data, volume 793.  John Wiley & Sons.
  * Lloyd, (1982) Lloyd, S. (1982).  Least squares quantization in pcm.  IEEE Transactions on Information Theory, 28(2):129–137.
  * MacQueen, (1967) MacQueen, J. B. (1967).  Some methods for classification and analysis of multivariate observations.  In Cam, L. M. L. and Neyman, J., editors, Proc. of the fifth Berkeley Symposium on Mathematical Statistics and Probability, volume 1, pages 281–297. University of California Press.
  * Marzouk and Xiu, (2009) Marzouk, Y. and Xiu, D. (2009).  A stochastic collocation approach to bayesian inference in inverse problems.  PRISM: NNSA Center for Prediction of Reliability, Integrity and Survivability of Microsystems, 6.
  * Marzouk et al., (2007) Marzouk, Y. M., Najm, H. N., and Rahn, L. A. (2007).  Stochastic spectral methods for efficient bayesian solution of inverse problems.  Journal of Computational Physics, 224(2):560–586.
  * McElreath, (2020) McElreath, R. (2020).  Statistical Rethinking: A Bayesian Course with Examples in R and STAN.  Chapman & Hall/CRC Texts in Statistical Science. CRC Press.
  * Meyer et al., (2021) Meyer, L., Pottier, L., Ribés, A., and Raffin, B. (2021).  Deep surrogate for direct time fluid dynamics.  CoRR, abs/2112.10296.
  * Mitra and Hlavacek, (2019) Mitra, E. D. and Hlavacek, W. S. (2019).  Parameter estimation and uncertainty quantification for systems biology models.  Current Opinion in Systems Biology, 18:9–18.
  * Modrák et al., (2022) Modrák, M., Moon, A. H., Kim, S., Bürkner, P., Huurre, N., Faltejsková, K., Gelman, A., and Vehtari, A. (2022).  Simulation-based calibration checking for bayesian computation: The choice of test quantities shapes sensitivity.
  * Mohammadi et al., (2018) Mohammadi, F., Kopmann, R., Guthke, A., Oladyshkin, S., and Nowak, W. (2018).  Bayesian selection of hydro-morphodynamic models under computational time constraints.  Advances in Water Resources, 117:53–64.
  * Murphy, (2007) Murphy, K. P. (2007).  Conjugate bayesian analysis of the gaussian distribution.  -.
  * Murphy, (2012) Murphy, K. P. (2012).  Machine Learning: A Probabilistic Perspective.  The MIT Press.
  * Neal, (2011) Neal, R. (2011).  Handbook of Markov Chain Monte Carlo.  Chapman and Hall/CRC.
  * Oladyshkin and Nowak, (2012) Oladyshkin, S. and Nowak, W. (2012).  Data-driven uncertainty quantification using the arbitrary polynomial chaos expansion.  Reliability Engineering & System Safety, 106:179–190.
  * Petersen and Pedersen, (2012) Petersen, K. B. and Pedersen, M. S. (2012).  The matrix cookbook.  Version 20121115.
  * Piironen et al., (2020) Piironen, J., Paasiniemi, M., and Vehtari, A. (2020).  Projective inference in high-dimensional problems: Prediction and feature selection.  Electronic Journal of Statistics, 14(1):2155 – 2197.
  * Psaros et al., (2023) Psaros, A. F., Meng, X., Zou, Z., Guo, L., and Karniadakis, G. E. (2023).  Uncertainty quantification in scientific machine learning: Methods, metrics, and comparisons.  Journal of Computational Physics, 477:111902.
  * Radev et al., (2023) Radev, S. T., Schmitt, M., Pratz, V., Picchini, U., Köthe, U., and Bürkner, P. (2023).  Jana: Jointly amortized neural approximation of complex bayesian models.  In Evans, R. J. and Shpitser, I., editors, Uncertainty in Artificial Intelligence, UAI 2023, July 31 - 4 August 2023, Pittsburgh, PA, USA, volume 216 of Proceedings of Machine Learning Research, pages 1695–1706. PMLR.
  * Raissi et al., (2019) Raissi, M., Perdikaris, P., and Karniadakis, G. (2019).  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707.
  * Ranftl and von der Linden, (2021) Ranftl, S. and von der Linden, W. (2021).  Bayesian surrogate analysis and uncertainty propagation.  Physical Sciences Forum, 3:6.
  * Rasmussen and Williams, (2005) Rasmussen, C. E. and Williams, C. K. I. (2005).  Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning).  The MIT Press.
  * Renardy et al., (2018) Renardy, M., Yi, T.-M., Xiu, D., and Chou, C.-S. (2018).  Parameter uncertainty quantification using surrogate models applied to a spatial model of yeast mating polarization.  PLOS Computational Biology, 14(5):1–26.
  * Robert and Casella, (2005) Robert, C. and Casella, G. (2005).  Monte Carlo Statistical Methods.  Springer Texts in Statistics. Springer New York.
  * Shao et al., (2017) Shao, Q., Younes, A., Fahs, M., and Mara, T. A. (2017).  Bayesian sparse polynomial chaos expansion for global sensitivity analysis.  Computer Methods in Applied Mechanics and Engineering, 318:474–496.
  * Smith, (2013) Smith, R. C. (2013).  Uncertainty Quantification: Theory, Implementation, and Applications.  Society for Industrial and Applied Mathematics, USA.
  * Stan Development Team, (2022) Stan Development Team (2022).  Stan Modeling Language Users Guide and Reference Manual.  Version 2.30.
  * Sudret, (2008) Sudret, B. (2008).  Global sensitivity analysis using polynomial chaos expansions.  Reliability Engineering & System Safety, 93(7):964–979.  Bayesian Networks in Dependability.
  * Säilynoja et al., (2022) Säilynoja, T., Bürkner, P.-C., and Vehtari, A. (2022).  Graphical test for discrete uniformity and its applications in goodness-of-fit evaluation and multiple sample comparison.  Statistics and Computing, 32.
  * Talts et al., (2018) Talts, S., Betancourt, M., Simpson, D., Vehtari, A., and Gelman, A. (2018).  Validating Bayesian Inference Algorithms with Simulation-Based Calibration.  arXiv e-prints, page arXiv:1804.06788.
  * Tarakanov and Elsheikh, (2019) Tarakanov, A. and Elsheikh, A. H. (2019).  Regression-based sparse polynomial chaos for uncertainty quantification of subsurface flow models.  Journal of Computational Physics, 399:108909.
  * van der Vaart, (2000) van der Vaart, A. (2000).  Asymptotic Statistics.  Asymptotic Statistics. Cambridge University Press.
  * Vehtari et al., (2017) Vehtari, A., Gelman, A., and Gabry, J. (2017).  Practical bayesian model evaluation using leave-one-out cross-validation and waic.  Statistics and Computing, 27.
  * Vehtari et al., (2021) Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C. (2021).  Rank-normalization, folding, and localization: An improved $\widehat{R}$ for assessing convergence of MCMC.  Bayesian Analysis, 16(2).  arXiv: 1903.08008.
  * Vehtari and Ojanen, (2012) Vehtari, A. and Ojanen, J. (2012).  A survey of Bayesian predictive methods for model assessment, selection and comparison.  Statistics Surveys, 6(none):142 – 228.
  * Waddell et al., (2002) Waddell, P., Kishino, H., and Ota, R. (2002).  Very fast algorithms for evaluating the stability of ml and bayesian phylogenetic trees from sequence data.  Genome informatics. International Conference on Genome Informatics, 13:82–92.
  * Wiener, (1938) Wiener, N. (1938).  The Homogeneous Chaos.  American Journal of Mathematics, 60(4):897–936.
  * Zeng et al., (2012) Zeng, L., Shi, L., Zhang, D., and Wu, L. (2012).  A sparse grid bayesian method for contaminant source identification.  Advances in Water Resources, 37:1–9.
  * Zhang et al., (2020) Zhang, J., Zheng, Q., Chen, D., Wu, L., and Zeng, L. (2020).  Surrogate-based bayesian inverse modeling of the hydrological system: An adaptive approach considering surrogate approximation error.  Water Resources Research, 56(1):e2019WR025721.  e2019WR025721 2019WR025721.
  * Zhu and Zabaras, (2018) Zhu, Y. and Zabaras, N. (2018).  Bayesian deep convolutional encoder–decoder networks for surrogate modeling and uncertainty quantification.  Journal of Computational Physics, 366:415–447.

##  Appendix A Notation

  * •

$\omega$: input

  * •

$y$: output

  * •

$\mathcal{M}$: simulator

  * •

$\widetilde{\mathcal{M}}$: surrogate model

  * •

$c$: surrogate approximation parameters

  * •

$d$: number of surrogate parameters

  * •

T-Step: surrogate training step

    * –

$N_{T}$: number of simulation training pairs

    * –

$\omega_{T}$: simulation input

    * –

$y_{T}$: simulation output

    * –

$e_{S}$: simulation noise

    * –

$\sigma_{S}$: simulation noise hyperparamters

    * –

$e_{A}$: approximation error

    * –

$\sigma_{A}$: surrogate approximation error hyperparameters

    * –

$\theta=(c,\sigma_{A})$: trainable surrogate parameters

    * –

$S$: number of T-posterior samples

  * •

I-Step: surrogate inference step

    * –

$N_{I}$: number of measurements

    * –

$\omega_{I}$: quantity of interest (QoI)

    * –

$y_{I}$: measurement data

    * –

$e_{I}$: measurement noise

    * –

$\sigma_{I}$: measurement noise hyperparameters

    * –

$K$: number of I-posterior samples

##  Appendix B Derivations

###  B.1 Inequality of E-Lik and E-Post via counterexample

In the following, we show the inequality of E-Post (Section 2.2.2) and E-Lik
(Section 2.2.3) using a counterexample with discrete random variables. We
denote the T-posterior as $p(\theta)$ and omit the dependence on the training
data $\mathcal{D}_{T}$ to simplify notation.

Let $\omega\in\\{0,1\\}$, $y\in\\{0,1\\}$, $\theta\in\\{0,1\\}$. We set the
probability values

$\begin{array}[]{c|c}&\\\ \hline\cr p(\omega=0)&1/2\\\
p(\omega=1)&1/2\end{array}\qquad\begin{array}[]{c|c}&\\\ \hline\cr
p(\theta=0)&1/2\\\
p(\theta=1)&1/2\end{array}\quad\begin{array}[]{c|c|c|c|c}&\omega=0,\theta=0&\omega=1,\theta=0&\omega=0,\theta=1&\omega=1,\theta=1\\\
\hline\cr p(y=0\mid\omega,\theta)&1/4&1/2&1/2&1/2\\\
p(y=1\mid\omega,\theta)&3/4&1/2&1/2&1/2\end{array}$

First, we state the general formulation of E-Lik and E-Post for discrete
pmf’s:

$\displaystyle p(\omega\mid y,u=\text{E-Lik})$
$\displaystyle=\frac{p(\omega)\sum_{i}p(y\mid\omega,\theta=i)p(\theta=i)}{\sum_{i}\sum_{j}p(y\mid\omega=j,\theta=i)p(\omega=j)p(\theta=i)}$
$\displaystyle p(\omega\mid y,u=\text{E-Post})$
$\displaystyle=\sum_{i}\frac{p(\omega)p(y\mid\omega,\theta=i)p(\theta=i)}{\sum_{j}p(y\mid\omega=j,\theta=i)p(\omega=j)}\frac{p(\omega=0)p(y=0\mid\omega=0,\theta=1)p(\theta=1)}{\sum_{j}p(y\mid\omega=j,\theta=1)p(\omega=j)}$

We calculate E-Post for $\omega=0$ and $y=0$ by plugging in the probability
values:

$\displaystyle p(\omega=0\mid y=0,u=\text{E-Post})$
$\displaystyle=\sum_{i}\frac{p(\omega=0)p(y=0\mid\omega=0,\theta=i)p(\theta=i)}{\sum_{j}p(y=0\mid\omega=j,\theta=i)p(\omega=j)}$
$\displaystyle=\frac{1/2\cdot 1/4\cdot 1/2}{3/8}+\frac{1/2\cdot 1/2\cdot
1/2}{1/2}=\frac{1/16}{3/8}+\frac{1/8}{1/2}$ $\displaystyle=\frac{5}{12}\approx
0.417$

Next, we calculate E-Lik for $\omega=0$ and $y=0$:

$\displaystyle p(\omega=0\mid y=0,u=\text{E-Lik})$
$\displaystyle=\frac{p(\omega=0)\sum_{i}p(y=0\mid\omega=0,\theta=i)p(\theta=i)}{\sum_{i}\sum_{j}p(y=0\mid\omega=j,\theta=i)p(\omega=j)p(\theta=i)}$
$\displaystyle=\frac{1/2\cdot(1/4\cdot 1/2+1/2\cdot 1/2)}{7/16}$
$\displaystyle=\frac{3}{7}\approx 0.429.$

We see that the results produced by E-Post and E-Lik for $\omega=0$ and $y=0$
are similar, but not equal.

The normalization constants for E-Post were given by:

$\displaystyle\sum_{j}p(y\mid\omega=j,\theta=0)p(\omega=j)$
$\displaystyle=p(y=0\mid\omega=0,\theta=0)p(\omega=0)+p(y=0\mid\omega=1,\theta=0)p(\omega=1)$
$\displaystyle=1/4\cdot 1/2+1/2\cdot 1/2=3/8$

and

$\displaystyle\sum_{j}p(y\mid\omega=j,\theta=1)p(\omega=j)$
$\displaystyle=p(y=0\mid\omega=0,\theta=1)p(\omega=0)+p(y=0\mid\omega=1,\theta=1)p(\omega=1)$
$\displaystyle=1/2\cdot 1/2+1/2\cdot 1/2=1/2.$

The normalization constant for E-Lik was given by:

$\displaystyle p_{\text{norm}}(y=0)=$
$\displaystyle\sum_{i}\sum_{j}p(y=0\mid\omega=j,\theta=i)p(\omega=j)p(\theta=i)$
$\displaystyle=$ $\displaystyle
p(y=0\mid\omega=0,\theta=0)p(\omega=0)p(\theta=0)+p(y=0\mid\omega=1,\theta=0)p(\omega=1)p(\theta=0)+$
$\displaystyle
p(y=0\mid\omega=0,\theta=1)p(\omega=0)p(\theta=1)+p(y=0\mid\omega=1,\theta=1)p(\omega=1)p(\theta=1)$
$\displaystyle=$ $\displaystyle 1/4\cdot 1/2\cdot 1/2+1/2\cdot 1/2\cdot
1/2+1/2\cdot 1/2\cdot 1/2+1/2\cdot 1/2\cdot 1/2=7/16$

###  B.2 Equivalence of E-Log-Lik and E-Log-Post

Here we show that the two formulations E-Log-Post and E-Log-Lik (Section
2.2.4) are equivalent. The E-Log-Lik is defined as:

$\displaystyle\log(p(\omega_{I}\mid y_{I},u=\text{E-Log-Lik}))$
$\displaystyle\propto\log(p(\omega_{I}))+\int\log(p(y_{I}\mid\omega_{I},\theta))p(\theta\mid\mathcal{D}_{T})d\theta$
(54)

Similarly, to E-Post and E-Lik we can define the E-Log-Post:

$\displaystyle\log(p(\omega_{I}\mid y_{I},u=\text{E-Log-
Post}))=\int\log(p(\omega_{I}\mid y_{I}))p(\theta\mid y_{T})d\theta$ (55)
$\displaystyle=\int\left[\log(p(\omega_{I}))+\log(p(y_{I}\mid\omega_{I},\theta)-\log(C(\theta))\right]p(\theta\mid
y_{T})d\theta$ (56) $\displaystyle=\log(p(\omega_{I}))\int p(\theta\mid
y_{T})d\theta+\int\log(p(y_{I}\mid\omega_{I},\theta))p(\theta\mid
y_{T})d\theta-\int\log(C(\theta))p(\theta\mid y_{T})d\theta$ (57)
$\displaystyle\propto\log(p(\omega_{I}))+\int\log(p(y_{I}\mid\omega_{I},\theta))p(\theta\mid
y_{T})d\theta$ (58)

In the last equation we used that the integral over a probability distribution
is one, i.e. $\int p(\theta\mid y_{T})d\theta=1$ and since the posterior is a
function of $\omega_{I}$ we can define a new constant
$C_{1}(\theta):=-\int\log(C(\theta))p(\theta\mid y_{T})d\theta$.

###  B.3 Case Study 1: Slope Intercept Model

In this section we derive the E-Log-Lik for the slope intercept model stated
in Section 3.1. Further we show the E-Lik with full normalization constant.
Let $c=[c_{1},c_{2}]^{T}$, Ω_T = [1 ωT(1)1 ωT(2)⋮⋮1 ωT(NT)],
$\hat{\omega}_{I}=[1,\omega_{I}]^{T}$,
$\mu_{T1}=[\mu_{T1}^{(1)},\mu_{T2}^{(2)}]^{T}$, Σ_T1 = [ΣT1(1, 1)ΣT1(1,
2)ΣT1(2, 1)ΣT1(2, 2)],

####  B.3.1 E-Log-Lik

Now we derive the stated Expected-Log-Likelihood result:

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Log-Lik})\propto{}$
$\displaystyle
p(\omega_{I})\exp\left\\{\int\log(p(y_{I}\mid\omega_{I},c))p(c\mid\mathcal{D}_{T})dc\right\\}$
(59) $\displaystyle={}$ $\displaystyle
p(\omega_{I})\exp\left\\{\int\log(\mathcal{N}(y_{I}\mid\hat{\omega}_{I}^{\top}c,\sigma_{I}^{2}))\mathcal{N}(c\mid\mu_{T1},\Sigma_{T1})dc\right\\}$
(60) $\displaystyle={}$ $\displaystyle
p(\omega_{I})\exp\left\\{\int-\frac{1}{2\sigma_{I}^{2}}(\hat{\omega}_{I}^{\top}c-y_{I})^{2}\mathcal{N}(c\mid\mu_{T1},\Sigma_{T1})dc\right\\}$
(61) $\displaystyle={}$ $\displaystyle
p(\omega_{I})\exp\left\\{-\frac{1}{2\sigma_{I}^{2}}\int(y_{I}^{2}-2y_{I}\hat{\omega}_{I}^{\top}c+c^{\top}\hat{\omega}_{I}\hat{\omega}_{I}^{\top}c)\mathcal{N}(c\mid\mu_{T1},\Sigma_{T1})dc\right\\}$
(62) $\displaystyle={}$ $\displaystyle
p(\omega_{I})\exp\left\\{-\frac{1}{2\sigma_{I}^{2}}(y_{I}^{2}-2y_{I}\hat{\omega}_{I}^{\top}\mathbb{E}[c]+\mathbb{E}[c^{\top}\hat{\omega}_{I}\hat{\omega}_{I}^{\top}c])\right\\}$
(63) $\displaystyle={}$ $\displaystyle
p(\omega_{I})\exp\left\\{-\frac{1}{2\sigma_{I}^{2}}(y_{I}^{2}-2y_{I}\hat{\omega}_{I}^{\top}\mu_{T1}+\text{Tr}(\hat{\omega}_{I}\hat{\omega}_{I}^{\top}\Sigma_{T1})+\mu_{T1}^{\top}\hat{\omega}_{I}\hat{\omega}_{I}^{\top}\mu_{T1})\right\\}$
(64)
$\displaystyle\begin{split}\propto{}&p(\omega_{I})\exp\left\\{-\frac{1}{2\sigma_{I}^{2}}(y_{I}^{2}-2y_{I}\hat{\omega}_{I}^{\top}\mu_{T1}+\omega_{I}^{2}\Sigma_{T1}^{(2,2)}+\omega_{I}(\Sigma_{T1}^{(2,1)}+\Sigma_{T1}^{(1,2)})+\right.\\\
&\left.\mu_{T1}^{(2)}\mu_{T1}^{(2)}\omega_{I}^{2}+2\mu_{T1}^{(1)}\mu_{T1}^{(2)}\omega_{I})\right\\}\end{split}$
(65) $\displaystyle\propto{}$ $\displaystyle
p(\omega_{I})\exp\left\\{-\frac{1}{2\sigma_{I}^{2}}(\mu_{T1}^{(2)}\mu_{T1}^{(2)}+\Sigma_{T1}^{(2,2)})(\omega_{I}-\frac{y\mu_{T1}^{(2)}-\Sigma_{T1}^{(1,2)}-\mu_{T1}^{(1)}\mu_{T1}^{(2)}}{\mu_{T1}^{(2)}\mu_{T1}^{(2)}+\Sigma_{T1}^{(2,2)}})^{2}\right\\}$
(66) $\displaystyle\propto{}$
$\displaystyle\mathcal{N}(\omega_{I}\mid\mu_{I1},\sigma_{I1}^{2}),$ (67)

with

$\displaystyle\sigma_{I1}^{2}={}$
$\displaystyle(\sigma_{I0}^{-2}+\sigma_{I}^{-2}(\mu_{T1}^{(2)}\mu_{T1}^{(2)}+\Sigma_{T1}^{(2,2)}))^{-1},$
(68) $\displaystyle\mu_{I1}={}$
$\displaystyle\sigma_{I1}^{2}(\sigma_{I0}^{-2}\mu_{I0}+\sigma_{I}^{-2}(\mu_{T1}^{(2)}y_{I}-\Sigma_{T1}^{(1,2)}-\mu_{T1}^{(1)}\mu_{T1}^{(2)}))$
(69)

where we used (Petersen and Pedersen,, 2012, p. 43):

$\displaystyle\text{Tr}(\hat{\omega}_{I}\hat{\omega}_{I}^{\top}\Sigma_{T1})$
$\displaystyle=\omega_{I}^{2}\Sigma_{T1}^{(2,2)}+\omega_{I}(\Sigma_{T1}^{(2,1)}+\Sigma_{T1}^{(1,2)})+\Sigma_{T1}^{(1,1)}$
(70)

and

$\displaystyle\mu_{T1}^{\top}\hat{\omega}_{I}\hat{\omega}_{I}^{\top}\mu_{T1}$
$\displaystyle=\mu_{T1}^{(2)}\mu_{T1}^{(2)}\omega_{I}^{2}+2\mu_{T1}^{(1)}\mu_{T1}^{(2)}\omega_{I}+\mu_{T1}^{(1)}\mu_{T1}^{(1)}$
(71)

####  B.3.2 E-Lik

E-Lik with normalization constant reads as:

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Lik})$
$\displaystyle=\frac{p(\omega_{I})\int
p(y_{I}\mid\omega_{I},c)p(c\mid\mathcal{D}_{T})dc}{\int p(\omega_{I})\int
p(y_{I}\mid\omega_{I},c)p(c\mid\mathcal{D}_{T})dcd\omega_{I}}$ (72)
$\displaystyle=\frac{\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\int\mathcal{N}(y_{I}\mid\hat{\omega}_{I}^{\top}c,\sigma_{I}^{2})\mathcal{N}(c\mid\mu_{T1},\Sigma_{T1})dc}{\int\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\int\mathcal{N}(y_{I}\mid\hat{\omega}_{I}^{\top}c,\sigma_{I}^{2})\mathcal{N}(c\mid\mu_{T1},\Sigma_{T1})dcd\omega_{I}}$
(73)
$\displaystyle=\frac{\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\mathcal{N}(y_{I}\mid\hat{\omega}_{I}^{\top}\mu_{T1},\hat{\omega}_{I}^{\top}\Sigma_{T1}\hat{\omega}_{I}+\sigma_{I}^{2})}{\int\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\mathcal{N}(y_{I}\mid\hat{\omega}_{I}^{\top}\mu_{T1},\hat{\omega}_{I}^{\top}\Sigma_{T1}\hat{\omega}_{I}+\sigma_{I}^{2})d\omega_{I}}$
(74)

##  Appendix C Further results

###  C.1 Case Study 1: Slope only Model

In Section 3.1 we showed the I-posteriors for the linear model with a slope
and an intercept. In this section, we consider an even simpler setup, where
the simulator and surrogate model are both a slope only model:

$\displaystyle
y=\mathcal{M}(\omega)=\widetilde{\mathcal{M}}(\omega,c)=c\cdot\omega,$ (75)

where $\omega\in\mathbb{R}$ is the input, $y\in\mathbb{R}$ the output and
$c\in\mathbb{R}$ is the unknown parameter. We consider to propagate only the
surrogate approximation parameter, i.e. we set $\theta=c$. We set normal
priors and normal likelihoods.

##### T-Step

As our simulation data set we use one input-output pair
$\mathcal{D}_{T}=\\{\omega_{T},y_{T}\\}$ and $N_{T}=1$. We set a normal prior
on $c$ with mean $\mu_{T0}$ and variance $\sigma_{T0}^{2}$. We calculate the
posterior using the conjugate prior relation for a normal-normal model
(Murphy,, 2007, 2012):

$\displaystyle p(c\mid y_{T})$ $\displaystyle=\frac{p(c)p(y_{T}\mid
c)}{p(y_{T})}$ (76)
$\displaystyle\begin{split}&=\mathcal{N}(c\mid\mu_{T1},\sigma_{T1}^{2}),\\\
\sigma_{T1}^{2}&=(\sigma_{T0}^{-2}+\sigma_{A}^{-2}\omega_{T}^{2})^{-1},\\\
\mu_{T1}&=\sigma_{T1}^{2}(\sigma_{T0}^{-2}\mu_{T0}+\sigma_{A}^{-2}\omega_{T}y_{T})\end{split}$
(77)

##### I-Step: Point

Figure 8: Two-Step-Procedure for linear model with normal likelihoods and
normal priors.

We set a normal prior on $\omega_{I}$:

$\displaystyle
p(\omega_{I})=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2}).$ (78)

We use the following normal likelihood:

$\displaystyle p(y_{I}\mid\omega_{I},c)=\mathcal{N}(y_{I}\mid
c\cdot\omega_{I},\sigma_{I}^{2})$ (79)

We compute the mean of the T-Step posterior: $\bar{c}=\mu_{T1}$ Now we can
compute the Point I-posterior using the conjugate prior relation (Murphy,,
2007, 2012):

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{Point})$ $\displaystyle\propto
p(\omega_{I})p(y_{I}\mid\omega_{I},\bar{c})$ (80)
$\displaystyle=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\mathcal{N}(y_{I}\mid\mu_{T1}\cdot\omega_{I},\sigma_{I}^{2})$
(81)
$\displaystyle\begin{split}&=\mathcal{N}(\omega_{I}\mid\mu_{I1},\sigma_{I1}^{2})\\\
&\text{, with:}\\\
\sigma_{I1}^{2}&=(\sigma_{I0}^{-2}+\sigma_{I}^{-2}\mu_{T1}^{2})^{-1},\\\
\mu_{I1}&=\sigma_{I1}^{2}(\sigma_{I0}^{-2}\mu_{I0}+\sigma_{I}^{-2}\mu_{T1}y_{I})\end{split}$
(82)

##### I-Step: E-Lik

Derivation

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Lik})$ $\displaystyle\propto
p(\omega_{I})\int p(y_{I}\mid\omega_{I},c)p(c\mid\mathcal{D}_{T})dc$ (83)
$\displaystyle=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\int\mathcal{N}(y_{I}\mid\omega_{I}\cdot
c,\sigma_{I}^{2})\mathcal{N}(c\mid\mu_{T1},\sigma_{T1}^{2})dc$ (84)
$\displaystyle=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\mathcal{N}(y_{I}\mid\omega_{I}\cdot\mu_{T1},\omega_{I}^{2}\cdot\sigma_{T1}^{2}+\sigma_{I}^{2})$
(85)

##### I-Step: E-Log-Lik

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Log-Lik})$
$\displaystyle\propto
p(\omega_{I})\exp\left\\{\int\log(p(y_{I}\mid\omega_{I},c))p(c\mid\mathcal{D}_{T})dc\right\\}$
(86)
$\displaystyle=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\exp\left\\{\int\log(\mathcal{N}(y_{I}\mid\omega_{I}\cdot
c,\sigma_{I}^{2}))\mathcal{N}(c\mid\mu_{T1},\sigma_{T1}^{2})dc\right\\}$ (87)
$\displaystyle\propto\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\exp\left\\{\int\left(\frac{1}{2\sigma_{I}^{2}}(y_{I}-\omega_{I}c)^{2}\right)\mathcal{N}(c\mid\mu_{T1},\sigma_{T1}^{2})dc\right\\}$
(88)
$\displaystyle\begin{split}&=\mathcal{N}(\omega_{I}\mid\mu_{I1},\sigma_{I1}^{2})\\\
&\text{, with:}\\\
\sigma_{I1}^{2}&=(\sigma_{I0}^{-2}+\sigma_{I}^{-2}(\mu_{T1}^{2}+\sigma_{T1}^{2}))^{-1},\\\
\mu_{I1}&=\sigma_{I1}^{2}(\sigma_{I0}^{-2}\mu_{I0}+\sigma_{I}^{-2}\mu_{T1}y_{I})\end{split}$
(89)

##### I-Step: E-Post

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Post})$
$\displaystyle=p(\omega_{I})\int\frac{p(y_{I}\mid\omega_{I},c)p(c\mid\mathcal{D}_{T})}{\int
p(y_{I}\mid\omega_{I},c)p(\omega_{I})d\omega_{I}}dc$ (90)
$\displaystyle=\mathcal{N}(\omega_{I}\mid\mu_{I0},\sigma_{I0}^{2})\int\frac{\mathcal{N}(y_{I}\mid\omega_{I}\cdot
c,\sigma_{I}^{2})\mathcal{N}(c\mid\mu_{T1},\sigma_{T1}^{2})}{\mathcal{N}(y_{I}\mid
c\cdot\mu_{I0},c^{2}\cdot\sigma_{I0}^{2}+\sigma_{I}^{2})}dc$ (91)

Figure 9: I-posterior distributions for the slope-only surrogate with normal
priors/likelihoods. We used exemplary hyperparameters and data. We use the
four UPs to compuete the I-posterior while varying $\sigma_{A}$.

In Fig. 9 we compare the four uncertainty propagation methods for exemplary
parameters. To control for T-epistemic uncertainty we vary $\sigma_{A}$. In
general, we notice similar I-posteriors to the slope-intercept model (see
3.1). Additionally, for large $\sigma_{A}$ we observe that both E-Lik and
E-Post produce bimodal I-posteriors.

####  C.1.1 Derviation E-Log-Lik

We derive the Expected-Log-Likelihood result for the slope-only model:

$\displaystyle p(\omega_{I}\mid y_{I},u=\text{E-Log-Lik})\propto{}$
$\displaystyle
p(\omega_{I})\exp\left\\{\int\log(p(y_{I}\mid\omega_{I},c))p(c\mid\mathcal{D}_{T})dc\right\\}$
(92) $\displaystyle={}$ $\displaystyle
p(\omega_{I})\exp\left\\{\int\log(\mathcal{N}(y_{I}\mid\omega_{I}c,\sigma_{I}^{2}))\mathcal{N}(c\mid\mu_{T1},\sigma_{T1})dc\right\\}$
(93) $\displaystyle={}$ $\displaystyle
p(\omega_{I})\exp\left\\{\int-\frac{1}{2\sigma_{I}^{2}}(\omega_{I}c-y_{I})^{2}\mathcal{N}(c\mid\mu_{T1},\sigma_{T1})dc\right\\}$
(94) $\displaystyle={}$ $\displaystyle
p(\omega_{I})\exp\left\\{-\frac{1}{2\sigma_{I}^{2}}\int(y_{I}^{2}-2y_{I}\omega_{I}c+c^{2}\omega_{I}^{2})\mathcal{N}(c\mid\mu_{T1},\sigma_{T1})dc\right\\}$
(95) $\displaystyle={}$ $\displaystyle
p(\omega_{I})\exp\left\\{-\frac{1}{2\sigma_{I}^{2}}(y_{I}^{2}-2y_{I}\omega_{I}\mathbb{E}[c]+\omega_{I}\mathbb{E}[c^{2}])\right\\}$
(96) $\displaystyle={}$ $\displaystyle
p(\omega_{I})\exp\left\\{-\frac{1}{2\sigma_{I}^{2}}(y_{I}^{2}-2y_{I}\omega_{I}\mu_{T1}+\mu_{T1}^{2}+\sigma_{T1}^{2})\right\\}$
(97) $\displaystyle\propto{}$ $\displaystyle
p(\omega_{I})\exp\left\\{-\frac{1}{2\sigma_{I}^{2}}(\mu_{T1}^{2}+\sigma_{T1}^{2})(\omega_{I}-\frac{\mu_{T1}y_{I}}{\mu_{T1}^{2}+\sigma_{T1}^{2}})^{2}\right\\}$
(98)
$\displaystyle\begin{split}\propto{}&\mathcal{N}(\omega_{I}\mid\mu_{I1},\sigma_{I1}^{2})\\\
&\text{, with:}\\\
\sigma_{I1}^{2}&=(\sigma_{I0}^{-2}+\sigma_{I}^{-2}(\mu_{T1}^{2}+\sigma_{T1}^{2}))^{-1},\\\
\mu_{I1}&=\sigma_{I1}^{2}(\sigma_{I0}^{-2}\mu_{I0}+\sigma_{I}^{-2}\mu_{T1}y_{I})\end{split}$
(99)

###  C.2 Case Study 2: Logistic Model

###  C.3 PCE posterior distributions

We show I-posterior density plots using the PCE surrogate to approximate the
logistic function (same setup as in Section 3.2.2) for additional true input
parameters $\omega_{I}^{*}=\\{-0.05,0.1,0.4\\}$. In Fig. 10 we propagate only
epistemic uncertainty and in Fig. 11 we propagate both epistemic and aleatoric
uncertainty.

Figure 10: Additional results for T-epistemic UP ($\theta=c$) with PCE
surrogate in Case study 2. Left: For $N_{T}=\\{10,30,50\\}$ the training data
set $\mathcal{D}_{T}$ (black dots), the T-posterior predictive distribution
and the mean of the T-posterior predictive distribution (dark red and red
lines) is shown. Right: For the true inputs
$\omega_{I}^{*}=\\{-0.05,0.1,0.4\\}$ (black vertical lines) we depict the
I-posterior distributions for each Point, E-Lik, E-Post, and E-Log-Lik
(colored lines). Figure 11: Additional results for T-epistemic and T-aleatoric
UP ($\theta=(c,\sigma_{A})$) with PCE surrogate in Case study 2. Left: For
$N_{T}=\\{10,30,50\\}$ the training data set $\mathcal{D}_{T}$ (black dots),
the T-posterior predictive distribution and the mean of the T-posterior
predictive distribution (dark red and red lines) is shown. Right: For the true
inputs $\omega_{I}^{*}=\\{-0.05,0.1,0.4\\}$ (black vertical lines) we depict
the I-posterior distributions for each Point, E-Lik, E-Post, and E-Log-Lik
(colored lines).

####  C.3.1 Computational Complexity

Figure 12: Sampling times estimated for Point, E-Lik, E-Log-Lik and E-Post
shown on a logarithmic scale. For each method, we sample 4 chains with 1000
warmup iterations and 2000 total iterations. We vary the maximum degree of
polynomials $d\in\\{2,5,9,25\\}$ and the number of clusters
$K\in\\{2,5,25,100,250,500,750,1000\\}$.

To evaluate the computational complexity of the different UP methods, we
estimate their sampling times using a PCE surrogate approximating a logistic
function (see Section 3.2.2). For each UP method, we sample 4 chains with 1000
warmup iterations and 2000 total iterations. However, for a fair comparison,
for E-Post we run for each model a full warmup phase with 1000 iterations, but
sample only $2000/K$-times. The computations were not performed in parallel.
We vary the maximum degree of polynomials $d\in\\{2,5,9,25\\}$ and the number
of clusters $K\in\\{2,5,25,100,250,500,750,1000\\}$. We plot the times on a
logarithmic scale in Fig. 12.