Statistics and Computing (2025) 35:66
https://doi.org/10.1007/s11222-025-10597-8

ORIG INAL PAPER

Uncertainty quantification and propagation in surrogate-based
Bayesian inference

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

Received: 9 April 2024 / Accepted: 23 February 2025 / Published online: 13 March 2025
© The Author(s) 2025

Abstract
Surrogatemodels are statistical or conceptual approximations formore complex simulationmodels. 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 three detailed case studies for
linear and nonlinear real-world 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

1 Introduction

Simulations of complex phenomena are crucial in the natu-
ral sciences and engineering for different scenarios, e.g., for
gaining systemunderstanding, 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 solu-
tions 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 approxi-
mations of computationally expensive simulation. Examples
are polynomial chaos expansion (Wiener 1938; Sudret 2008;
Oladyshkin and Nowak 2012; Bürkner et al. 2023), Gaus-
sian processes (Kennedy and O’Hagan 2001; Rasmussen
and Williams 2005) or neural networks (Goodfellow et al.

B Philipp Reiser
philipp-luca.reiser@simtech.uni-stuttgart.de

1 Cluster of Excellence SimTech, University of Stuttgart,
Stuttgart, Germany

2 Department of Statistics, TU Dortmund University,
Dortmund, Germany

2016). Recently, there has been a great interest in applying
surrogate models in relevant areas, for example in hydrol-
ogy (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). Further-
more, great methodological advances have been made in the
field of surrogatemodeling, for example, the incorporation of
physical knowledge (Raissi et al. 2019; Li et al. 2021; Brand-
stetter 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. Conse-
quently, 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 uncer-
tainty quantification (UQ) have been developed (e.g., Shao
et al. 2017; Bürkner et al. 2023, ). Uncertainty propagation
(UP), a sub-field of UQ, is particularly important for address-
ing surrogate uncertainties in subsequent (surrogate-based)
inference tasks (Smith 2013; Lavin et al. 2021; Psaros et al.
2023).

123

http://crossmark.crossref.org/dialog/?doi=10.1007/s11222-025-10597-8&domain=pdf


66 Page 2 of 28 Statistics and Computing (2025) 35 :66

Fig. 1 Overview of two-step procedure. Left: In the surrogate training
step (T-Step), training data is generated using a simulator and a surro-
gate 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

Usingprobability theory as its fundamental basis,Bayesian
statistics provides a rigorous way for UQ generally and
specifically for UP (Gelman et al. 2013; McElreath 2020;
Bürkner et al. 2023). One important use-case for UP in sur-
rogatemodeling arises from the “forward” problem, inwhich
input parameters are uncertain and the goal is to propagate
them through the surrogate while accounting for its uncer-
tainty 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 ana-
lytic models are available. Similarly, 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.,
Kennedy and O’Hagan 2001; Marzouk et al. 2007; Marzouk
and Xiu 2009; Zeng et al. 2012; Laloy et al. 2013; Li and
Marzouk 2014; Cleary et al. 2021), rigorous propagation of
the surrogate uncertainty becomes even more relevant and
challenging.

The framework introduced by Kennedy and O’Hagan
(2001), for instance, attempts to jointly infer unknown input
and surrogate parameters using data from both simulation
models and observations. However, this approach leads to a
loss of control over which parameters are updated by which
data source, also implying complex posteriors that are hard
to sample from (Bayarri et al. 2009).

To address these issues, modularization (Bayarri et al.
2009) has been introduced to surrogate modeling, aiming to
update only specific parameters using selected data. Despite
these advances, the uncertainty of the surrogate parameters
is often neglected or simplified in the context of surrogate
modeling. For example, Bayarri et al. (2009) propagated only

point estimates of the surrogate parameters betweenmodules,
thus neglecting relevant uncertainty in subsequent calcula-
tions. Further, Zhang et al. (2020) applied surrogate-based
Bayesian inference to hydrological systems and propagated
parts of the surrogate uncertainty assuming normal surro-
gate posteriors. We review these methods in more detail in
Sect. 2.2.3.

In this paper, our primary focus lies on UP when solving
probabilistic inverse problems via surrogate models. Exist-
ing methods proposed for the same challenge are scarce and
propagate the surrogate-based uncertainty only in a selec-
tive 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) Within a formal framework for surrogate-
based Bayesian inference, we specify and categorize all
relevant uncertainties. (ii) We present three distinct meth-
ods for propagating surrogate uncertainty 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 cal-
ibration achieved via the different UP approaches. Finally,
we evaluate our methods in three detailed case studies.

2 Method

In the following, we propose a framework for uncertainty
propagation (UP) in surrogate-based Bayesian inference.

123



Statistics and Computing (2025) 35 :66 Page 3 of 28 66

Fig. 2 Graphical model of the T-Step and I-Step. Left: In the T-Step,
the observed quantities are simulation parametersωT , simulation output
yT , and the noise hyperparameters σS . The unknowns are the surrogate
parameters θ . Right: In the I-Step, measurement data yI is observed
NI times and S posterior samples of θ are propagated from the T-Step.
The dashed arrow indicates that uncertainty in θ is propagated to the
I-Step while θ is not updated using the data yI . The unknowns to be
inferred are the simulation parameters ωI and the measurement error
hyperparameters σI

This framework is applicable to general surrogate models,
e.g., linear regression, polynomial chaos expansion, Gaus-
sian 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), anoverviewof all uncertainties
occurring in surrogate-based inference, several UP schemes
of selected uncertainties, and their evaluation.

2.1 Two-step procedure

Wepropose 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 quantity of interest, as
illustrated in the graphical model (Jordan 1999; Ben-Gal
2008) shown in Fig. 2. 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 uncertain-
ties occurring in the two-step procedure to be discussed in
detail below.

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 induceduncertainties.

SimulatorWe consider a given (potentially stochastic) simu-
latorM, 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 ωT ,
we obtain the output response from the simulator as follows:

yT = M(ωT ; eS) with eS ∼ p(eS | σS), (1)

where eS is the noise of the simulation drawn from a sim-
ulator noise distribution p(eS | σS) with hyperparameters
σ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 determinis-
tic 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 pro-
cess. This setup induces the true generating distribution
p(yT | ωT , σS) of simulation responses yT given input sim-
ulation parameters ωT and noise hyperparameters σS .

Surrogate model A surrogate model ˜M is a statistical model
that aims to approximate the simulatorM ≈ ˜Mwhile being
computationally more efficient to evaluate. The parametric
form of the surrogate is defined by a set of surrogate approx-
imation parameters c. For example, if the surrogate model is
a polynomial, then c are the polynomial coefficients. Given
input parameters ωT , surrogate approximation parameters c,
and the simulator noise distribution p(eS | σS), we can cal-
culate the surrogate response:

ỹT = ˜M(ωT , c, eS) with eS ∼ p(eS | σS). (2)

Surrogate approximation error A surrogate model with a
fixed architecture will always bemisspecified, if the true sim-
ulator is not included in the class of surrogate models that
we define. This is the standard use case, since we are specif-
ically 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 approx-
imation error of the surrogate with respect to the simulator,
caused by limited expressibility of the surrogate, we intro-
duce an additional error term eA. We assume eA to follow a
distribution p(eA | σA) with hyperparameters σA. This dis-
tribution captures aleatoric uncertaintywhich is not reducible
by more training data. We can then model the true response
yT via a (potentially unknown) function f̃ that takes the out-
put of the surrogate ỹT and the approximation error eA:

yT = f̃ (ỹT , eA) with eA ∼ p(eA | σA). (3)

In practice, we might assume a simple additive error: yT =
ỹT +eA and a normal distribution for eA, but our framework is
agnostic to these choices. For the sake of readability, we use
θ = {c, σA} to combine all trainable surrogate parameters in
a single vector. Together, this implies a surrogate likelihood
of the simulator responses yT :

123



66 Page 4 of 28 Statistics and Computing (2025) 35 :66

Table 1 Uncertainties in the two-step procedure

Parameter/Posterior distribution T-Step I-Step Synonyms

Simulator noise hyperparameter σS Aleatoric – –

Surrogate approx error hyperparameter σA Aleatoric Aleatoric T-aleatoric uncertainty

Surrogate parameter Epistemic Aleatoric T-epistemic uncertainty /

Posterior p(c, σA | DT ) T-posterior

Measurement noise – Aleatoric –

Hyperparameter σI

Simulator-based – Epistemic –

Posterior p(ωI | yI ;M)

Surrogate-based – Epistemic & I-posterior

Posterior p(ωI | yI ; ˜M, u) aleatoric*

* Depends on uncertainty propagation method u.
For each parameter and posterior distribution, we list the type of uncertainty (epistemic or aleatoric) in the T-/I-Step. We also list the synonyms
used throughout the text

yT ∼ p(yT | ωT , σS, θ). (4)

Surrogate training To train the surrogate, we use the simu-
lator M to generate training data DT = {ω(i)

T , σ
(i)
S , y(i)

T }NT
i=1

consisting of NT inputs (ω
(i)
T , σ

(i)
S ) and corresponding out-

puts y(i)
T ∼ p(y(i)

T | ω
(i)
T , σ

(i)
S ). The goal of the T-Step is to

fit the parameters θ of the surrogate ˜M to approximate the
data distribution implied byM.

To train our surrogate parameters θ , we perform Bayesian
inference using the fast-to-evaluate surrogate likelihood
p(yT | ωT , σS, θ) and a (potentially non-informative) prior
p(θ). We obtain the joint posterior distribution over all sur-
rogate parameters given the simulation training data DT as:

p(θ | DT ) ∝
NT
∏

i=1

p(y(i)
T | ω

(i)
T , σ

(i)
S , θ) p(θ). (5)

This joint posterior describes the epistemic (reducible) uncer-
tainty in the surrogate model parameters. In the case of
infinite training data, i.e. NT→∞ the posterior p(θ |DT )

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 repre-
sent the posterior in the form of S posterior samples
{θ(1), . . . , θ (S)} ∼ p(θ | DT ). The framework is in prin-
ciple agnostic to the choice of estimation algorithm, as long
as the algorithm can be used to obtain posterior samples. This
includes MCMC but also other methods such as variational
inference (Kucukelbir et al. 2017) or integrated Laplace
approximation (Ruiz-Cárdenas et al. 2012; Martino and
Riebler 2019). For a surrogate model with non-identifiable

parameters θ , the algorithm would potentially have to deal
with multimodalities (Medina-Aguayo and Christen 2022).

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, i.e.
the input simulation parameters, 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
M, but with a potentially different measurement noise eI .
Measurement data yI is then generated from an unknown
underlying parameter ωI (which has the same dimension as
ωT and serves as input to the simulator):

yI = M(ωI ; eI ) with eI ∼ p(eI | σI ), (6)

where the measurement noise eI is drawn from a distribu-
tion p(eI | σI ) with hyperparameters σI that describe the
aleatoric uncertainty of themeasurement. This induces a gen-
erating distribution p(yI | ωI , σI ) ofmeasurements yI given
inputs ωI and noise hyperparameters σI . In contrast to the
simulator setup, we only have access to the measurements
yI , but the underlying true input ωI and σI are unknown.
Inference Given a set of real-world measurement data yI =
{y(i)

I }NI
i=1 with y(i)

I ∼ p(y(i)
I | ωI , σI ) for i = 1, . . . , NI , our

goal is to infer the posterior of the unknown parameters ωI ,
which constitute our primary quantity of interest. Addition-
ally, we can also infer σI although it is only of secondary
interest. To simplify the presentation, we drop σI from the
notation in the following.

If the simulator M had a tractable and easy-to-evaluate
likelihood function p(yI | ωI ,M), we could simply calcu-

123



Statistics and Computing (2025) 35 :66 Page 5 of 28 66

late the posterior of the unknown simulation parameters ωI

given the measurement data yI and a prior p(ωI ):

p(ωI | yI ,M) ∝ p(yI | ωI ,M) p(ωI ). (7)

However, for many real-world simulators, this likelihood is
unavailable or highly cumbersome to evaluate (see Sect. 1),
which is why we seek to replace it with simpler likelihood
of the previously trained surrogate ˜M. To incorporate the
uncertainties of the surrogate parameters θ = (c, σA) from
the T-Step into our real-world inference, we need to some-
how propagate their T-Step posterior p(θ | DT ) to the
I-step. From the perspective of the I-Step, the uncertainty in
p(θ | DT ) becomes aleatoric since it is no longer reducible.
The main question how to propagate p(θ | DT ) 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 ∈ U , for which we
can obtain a surrogate-based posterior p(ωI | yI , ˜M, u) as
detailed below.

2.2 Uncertainty propagation in the two-step
procedure

In the following, we present four different uncertainty prop-
agation 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 U = {Point,E-Post,E-Lik,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(θ | DT )

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(yI | ωI , ˜M, u) of the I-step. In the follow-
ing, all posteriors are calculated using the surrogate model
˜M and we will omit the dependency of ˜M in the posterior:
p(ωI | yI , u) = p(ωI | yI , ˜M, u) to improve readabil-
ity. Furthermore, for all uncertainty propagation methods,
we will assume that yI consists of conditional i.i.d. mea-
surements y(i)

I , such that the likelihood factorizes easily.
This assumption is not necessary for our framework, but
it simplifies the notation and allows for more efficient
computation.

2.2.1 Point estimate

First, we describe the I-Step using only a point estimator θ̂ of
the T-posterior p(θ | DT ), e.g., the posterior mean, median,
or mode. In this case, we condition the posterior of ωI on
θ̂ , i.e., we reduce the full posterior to a point estimate with-

out epistemic uncertainty. This is not a method we advocate
for, but rather use it as a baseline to compare against more
sophisticated methods. The posterior probability of ωI given
measurement data yI and θ̂ is then (up to normalizing con-
stants):

p(ωI | yI , u = Point) ∝ p(yI | ωI , θ̂ )p(ωI ), (8)

where the dependency of the posterior of ωI on θ̂ is rep-
resented via u = Point. With the i.i.d. assumption, the
(unnormalized) log I-posterior is

log p(ωI | yI , u = Point)

∝
NI
∑

i=1

log p(y(i)
I | ωI , θ̂ ) + log p(ωI ),

(9)

where we use the ∝ symbol for log-probability statements
to imply that an additive constant C (here, the log marginal
likelihood) is not shown in the equation. Such a constant is
independent of the parameters and thus irrelevant for pos-
terior inference via MCMC or related sampling methods.
The log-posterior can easily be specified in a probabilis-
tic programming language (Gorinova et al. 2019) such as
Stan (Carpenter et al. 2017) and we can then sample from
the posterior with MCMC, leading to K posterior samples
{ω(1)

I , . . . , ω
(K )
I } ∼ p(ωI | yI , u = 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(yI | ωI , θ̂ ). While we are
completely neglecting the epistemic uncertainty contained in
the T-posterior, we propagate (a point estimate of) the surro-
gate approximation error parameter σA through θ̂ . We expect
the Point I-posterior p(ωI | yI , u = Point) to be overconfi-
dent (too narrow); unless we have a sufficiently large amount
of simulation training dataDT such that p(θ | DT ) converges
to a point mass and then corresponds exactly to the point esti-
mator θ̂ .

Related work Training a surrogate model using simulation
data and subsequently employing its point estimate to infer
unknown input parameters from observed data is a stan-
dard and well known approach (e.g., Marzouk et al. 2007;
Laloy et al. 2013; Li and Marzouk 2014). This method has
been extensively studied in the context of Bayesian surro-
gate models, particularly through modularization introduced
by Bayarri et al. (2009).

2.2.2 Expected-posterior

Next, we present the Expected-Posterior (E-Post), a method
that propagates both aleatoric and epistemic uncertainty in
theT-posterior bymarginalizing over the surrogate parameter

123



66 Page 6 of 28 Statistics and Computing (2025) 35 :66

posterior p(θ | DT ) in the posterior of ωI :

p(ωI | yI , u = E-Post) =
∫

p(ωI | yI , θ)p(θ | DT )dθ.

(10)

Since we only have samples θ(s) ∼ p(θ | DT ) from the T-
step, we compute the Monte Carlo (MC) approximation of
E-Post as:

p(ωI | yI , u = E-Post)
MC≈ 1

S

S
∑

s=1

p(ωI | yI , θ(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 θ(s) ∼
p(θ | DT ) using the point I-posterior approach in Eq. (9),
which results in K draws of ωI , i.e., {ω(s,1)

I , . . . , ω
(s,K )
I }.

The combination of these draws across all s = 1, . . . , S then
represents an MC-estimate of the E-Post I-posterior as per
Eq. (11), with a total of S · K posterior draws. We note that
the number of propagated T-posterior draws S is a hyperpa-
rameter that needs to be tuned depending on the complexity
and dimensionality of the problem.
Related work The idea of constructing a posterior via the
aggregation of multiple posterior distributions, each approx-
imated 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 sub-
sequently averaging the resulting bootstrapped posteriors.
Similarly, in the context of missing value imputation (Lit-
tle and Rubin 2019), this approach has been used to combine
models fitted on multiple imputed datasets (Bürkner 2017,
2018). In terms of how many imputed data sets are needed,
Austin et al. (2021) report that between 20 and 100 imputa-
tions are typically used. In both multiple data imputation and
bagging, models are fitted to different datasets whereas in
our E-Post method, the data is the same but the model itself
changes. Furthermore, within the context of modularization
of Bayesian models (Plummer 2014; Jacob et al. 2017), this
approach is known as the cut distribution.

2.2.3 Expected-likelihood

Next, we present the Expected-Likelihood (E-Lik) approach,
where we marginalize over the T-posterior p(θ | DT ) in the
I-likelihood:

p(yI | ωI , u = E-Lik) =
∫

p(yI | ωI , θ)p(θ | DT )dθ.

(12)

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

p(ωI | yI , u = E-Lik) ∝ p(yI | ωI , u = E-Lik)p(ωI ).

(13)

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

log p(ωI | yI , u = E-Lik)

∝ log
∫ NI
∏

i=1

p(y(i)
I | ωI , θ)p(θ | DT )dθ + log p(ωI ),

(14)

for which we can obtain an MC approximation using draws
θ(s) ∼ p(θ | DT ):

log p(ωI | yI , u = E-Lik)

MC≈ log

(

1

S

S
∑

s=1

NI
∏

i=1

p(y(i)
I | ωI , θ

(s))

)

+ log p(ωI ).
(15)

This representation has the problem that the product over
likelihood components becomes numerically unstable as
NI 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

∑S
s=1 p(yI |ωI , θ

(s))

= log
∑S

s=1 exp(log p(yI | ωI , θ
(s))), which has a numer-

ically stable implementation (Carpenter et al. 2017). For
given θ(s), the joint log likelihood is then again a simple
sum: log p(yI | ωI , θ

(s)) = ∑NI
i=1 log p(y(i)

I | ωI , θ
(s)).

In contrast to E-Post, E-Lik is expressed as a single proba-
bilistic program and we can use MCMC to approximate its
I-posterior with K samples {ω(1)

I , . . . , ω
(K )
I }.

More general, in contrast to E-Post, which marginalizes
over θ in the I-posterior p(ωI | yI , θ), the E-Lik method
marginalizes over θ in the I-likelihood p(yI | ω, θ). E-Lik
and E-Post are not identical (see Appendix Section B.1 for
a counterexample), but we demonstrate in Sect. 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 approxi-
mation 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 predic-
tive performance, which integrates the likelihood over the

123



Statistics and Computing (2025) 35 :66 Page 7 of 28 66

posterior of the same model before taking the logarithm out-
side 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 theExpected-Log-Likelihood (E-Log-Lik) approach, instead
ofmarginalizingover the likelihood as inE-Lik,wemarginal-
ize over the T-posterior in the log-likelihood. We define the
E-Log-Lik I-likelihood as:

p(yI | ωI , u = E-Log-Lik)

:= exp

(∫

log(p(yI | ωI , θ)) p(θ | DT )dθ

)

.
(16)

The posterior of the E-Log-Lik for ωI is then defined as:

p(ωI | yI , u = E-Log-Lik)

∝ p(yI | ωI , u = E-Log-Lik)p(ωI ).
(17)

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

log p(ωI | yI , u = E-Log-Lik)

∝
∫

log

( NI
∏

i=1

p(y(i)
I | ωI , θ)

)

p(θ | DT )dθ + log p(ωI )

=
∫ NI
∑

i=1

log p(y(i)
I | ωI , θ) p(θ | DT )dθ + log p(ωI ),

(18)

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

log p(ωI | yI , u = E-Log-Lik)

MC≈ 1

S

S
∑

s=1

NI
∑

i=1

log p(y(i)
I | ωI , θ

(s)) + log p(ωI ),
(19)

which can be fitted with MCMC leading to K draws
{ω(1)

I , . . . , ω
(K )
I } in the I step. We can further rewrite the

MC approximation of the log-likelihood as

log p(yI | ωI , u = E-Log-Lik)

=
S
∑

s=1

NI
∑

i=1

log
(

p(y(i)
I | ωI , θ

(s))
1
S

)

.
(20)

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 gen-
eralizes to unequal weights as we illustrate in Sect. 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

p(ωI | yI , u = E-Log-Post)

∝ exp

(∫

log p(ωI | yI , θ) p(θ | DT )dθ

)

,
(21)

which however is equivalent to E-Log-Lik as we show in
Appendix B.1
Related work The computation of the E-Log-Lik is con-
ceptually 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 computa-
tional convenient although uncommonmeasure 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 {θ(s)}Ss=1
to get cluster centroids {μ(l)}Ll=1. The sufficient number of
clusters L for a trustworthy approximation of the I-posterior
highly depends on the complexity of the I-posterior and is
a hyperparameter that needs to be tuned. Here, we carried
out a visual convergence analysis, given that rigorous guide-
lines on the choice of the number of clusters are lacking to
date and an open scientific challenge. For each of the clus-
ter centroids we additionally store weights {α(l)}Ll=1, where
α(l) is the percentage of draws associated with the cluster.
While inducing an approximation error, 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 Sect. 3.

Clustering and re-weighting can be easily applied in all
UP methods by replacing θ(s) with μ(l) as well as replacing
the equal weight 1/S with α(l) after moving it inside the sum.
For example, when using E-Log-Lik, the MC-approximated

123



66 Page 8 of 28 Statistics and Computing (2025) 35 :66

integral becomes:

log p(yI | ωI , u = E-Log-Lik)

MC≈ 1

S

S
∑

s=1

NI
∑

i=1

log(p(y(i)
I | ωI , θ

(s))

≈
L
∑

l=1

NI
∑

i=1

α(l) log p(y(i)
I | ωI , μ

(l)).

(22)

2.2.6 Parallelization

Inference based on all introduced UP methods can be paral-
lelized, 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 par-
allelizable, 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 computa-
tional overhead.

In contrast, for E-Lik, we fit only one model during the
I-step, which loops over the T-posterior draws when evaluat-
ing its likelihood. This however defines a more complicated
posterior, which is substantially slower to sample from com-
pared to the individual E-Post models. Fitting the single
E-Lik model can be sped 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 contri-
butions 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 asE-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 impor-
tant step to improve the trustworthiness of any inference

algorithm.As such, it is a crucial aspect of theBayesianwork-
flow (Gelman et al. 2013). In our setup, we are specifically
interested in the uncertainty calibration of p(ωI | yI , ˜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 algo-
rithm such as MCMC.

Algorithm 1 SBC for the Two-Step Procedure
Choose simulator M, surrogate ˜M, and uncertainty propagation
method u
for m in Number T-Step trials do

Draw a training dataset DT using the simulator M
Fit the surrogate ˜M and calculate the posterior p(θ | DT ) (T-Step)
for n in Number I-Step trials do

Draw a prior sample ω∗
I ∼ p(ωI ), σI ∼ p(σI )

Draw measurements yI ∼ p(yI | ω∗
I , σI ,M) using simula-

tor M
Draw posterior samples {ω(1)

I , . . . , ω
(K )
I } ∼ p(ωI |

yI , ˜M, u) (I-Step)
Store the rank of ω∗

I within the set of posterior sam-

ples {ω(1)
I , . . . , ω

(K )
I }

end for
end for
Perform uniformity test on the stored ranks

We extend SBC to the two-step procedure and use the
notation of Bürkner et al. (2023). For any quantile q ∈ (0, 1),
let Uq(ωI | yI , ˜M, u) be any uncertainty region (e.g., the
quantile-based credible intervals or highest density intervals)
given by the posterior p(ωI | yI , ˜M, u) that depends on
the surrogate ˜M and the uncertainty propagation method u
(see Sect. 2.1 and 2.2). If the generating distribution of the
assumed surrogate ˜M is equal to the true data-generating
distribution induced by the simulator the following property
holds simultaneously for all q ∈ (0, 1):

q =
∫∫

I[ω∗
I ∈ Uq(ωI | yI , ˜M, u)] p(yI | ω∗

I ,M)

× p(ω∗
I ) dyI dω∗

I ,

(23)

where we ignore σI for simplicity and I[·] denotes the indica-
tor 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
ωI . We calculate the ranks as r(ω∗

I , {ω(1)
I , . . . , ω

(K )
I }) =

∑K
k=1 I

[

ω∗
I ≤ ω

(k)
I

]

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(γ )-statistic
(Säilynoja et al. 2022) as a quantitativemeasure of uniformity

123



Statistics and Computing (2025) 35 :66 Page 9 of 28 66

allowing for a faster comparison of calibration (or strength
of miscalibration) between different methods. Checking
the calibration of a specific uncertainty region, e.g., U0.95

(ωI | yI , ˜M, u), corresponds to evaluating Eq. (23) 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 simulatorM: (i) Sample a simulation
input parameterω∗

I ∼ p(ωI ), (ii) conditioned onω∗
I generate

measurement output data yI ∼ p(yI | ω∗
I ), (iii) using the

simulator M and given the measurements yI draw samples
{ω(1)

I , . . . , ω
(K )
I }, and (iv) using the posterior samples, we

calculate the rank statistics r(ω∗
I , {ω(1)

I , . . . , ω
(K )
I }) from the

posterior p(ωI | yI ,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 Sect. 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 DT used
to train the surrogate ˜M. Within each such T-Step trial, we
repeat the I-Step for SBC: First, we draw a sample of the
simulation parameters ω∗

I and the simulation noise hyperpa-
rameter σI from their respective priors. Next, measurement
data yI is generated using the simulatorM. Using the surro-
gate ˜M and the UP method u, we draw posterior samples
{ω(1)

I , . . . , ω
(K )
I } from the I-posterior and store the ranks

r(ω∗
I , {ω(1)

I , . . . , ω
(K )
I }). The repetition of the T-Step and I-

Step, in addition to covering the whole input space of ω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

(i) Incorrect implementation of simulator M
(ii) Incorrect implementation of probabilistic program of

the surrogate
(iii) Problems with the sampling algorithm

• Additional scenarios in proposed SBC for surrogate-
based inference

(iv) Inflexible surrogate
(v) Insufficient training of surrogate because of too little

simulation training data
(vi) Inappropriate uncertainty propagation in the surrogate-

based inference.

Regarding points (i)–(iii), we assume the simulator to be cor-
rect 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., σA → 0) and an
infinite amount of training samples (NT → ∞) in order to
remove the epistemic uncertainty in the posterior p(θ | DT ).
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 surrogatemodels ˜M and differentUP
methods u by comparing their SBC results, either graphically
or via test statistics.

3 Experiments

We evaluate our surrogate-based Bayesian inference frame-
work in three case studies: (1) A linear setup, where we
propagate only epistemic uncertainty, (2) a nonlinear setup,
in which we propagate both epistemic and aleatoric uncer-
tainty, and (3) a real-world model. All code and material can
be found on GitHub.1

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:

yT = M(ωT , σS = 0) = a + bωT , (24)

with simulation input parameter ωT and two simulation con-
trol parameters a, b set to a = 0.5, b = 2, and a simulation
noise parameter σS set to zero, i.e., we use a deterministic
simulator. For theT-Step,we consider NT = 2 trainingpoints
DT = {ω(i)

T , y(i)
T }NT

i=1 (chosen according to Table 2), where

ω
(i)
T denotes the simulation input and y(i)

T = M(ω
(i)
T ) the cor-

responding simulation output. We denote by �T = [ 1 ω
(1)
T

1 ω
(2)
T

]

the designmatrix of all input parameters and by yT 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:

˜M(ωT , c) = c1 + c2ωT , (25)

where the intercept c1 and slope c2 form the surrogate approx-
imation parameters c = [c1, c2]
. For the T-Step, we only
consider the surrogate approximation parameters c as train-
able surrogate parameters. Even though the surrogate can

1 https://github.com/philippreiser/bayesian-surrogate-uncertainty-
paper.

123

https://github.com/philippreiser/bayesian-surrogate-uncertainty-paper
https://github.com/philippreiser/bayesian-surrogate-uncertainty-paper


66 Page 10 of 28 Statistics and Computing (2025) 35 :66

match the simulator perfectly, we fix the surrogate approxi-
mation error hyperparameter σ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:

p(c) = N (c | μT 0, �T 0), (26)

with mean μT 0 = 0 and a covariance matrix �T0 = σ 2
T 0 I

where I is the identity matrix. We set the T-likelihood to a
normal distribution as well:

p(yT | c) =
NT
∏

i=1

N (y(i)
T | c1 + c2ω

(i)
T , σ 2

A). (27)

To generate data for the I-Step, a single measurement yI =
M(ω∗

I ) is obtained by inputting a true input parameter ω∗
I

to the deterministic simulator M. In the following, we will
also use the augmented vector ω̂I = [1, ωI ]T to simplify
notation. Within the I-step, we set a normal prior with mean
μI0 and variance σ 2

I0 on ωI :

p(ωI ) = N (ωI | μI0, σ
2
I0). (28)

The I-likelihood p(yI | ωI ) is assumed to be normal with
fixed variance σ 2

I :

p(yI | ωI , c) = N (yI | c1 + c2ωI , σ
2
I ) (29)

Deriving the posteriors In the following, we perform the T-
Step (see Sect. 2.1.1) and I-Step (see Sect. 2.2) for the linear
setup. In the T-Step, we calculate the T-posterior of the surro-
gate 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-StepWe calculate the T-posterior for the surrogate approx-
imation parameters c given the training data DT from the
simulator by using the conjugate prior relation for a normal-
normal model (Murphy 2007, 2012) leading to a normal
posterior:

p(c | DT ) = N (c | μT 1, �T 1), (30)

with

�T 1 = (�−1
T 0 + σ−2

A �

T �T )−1, (31)

μT 1 = �T 1(�
−1
T0μT 0 + σ−2

A �T yT ). (32)

We see that by fixing σ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 σA itself.
I-Step Below, we derive the I-posteriors of all four uncer-
tainty propagation procedures. We compute the mean of the
T-posterior: c̄ = μT 1 and use it to calculate the Point I-
posterior as follows:

p(ωI | yI , u = Point) ∝ p(ωI )p(yI | ωI , c̄) (33)

= N (ωI | μI0, σ
2
I0)N (yI | μ

(1)
T 1 + μ

(2)
T 1ωI , σ

2
I ) (34)

∝ N (ωI | μI1, σ
2
I1), (35)

with

σ 2
I1 = (σ−2

I0 + σ−2
I μ

(2)
T 1μ

(2)
T 1)

−1, (36)

μI1 = σ 2
I1(σ

−2
I0 μI0 + σ−2

I μ
(2)
T 1(yI − μ

(1)
T 1)). (37)

Again, this works due to the normal-normal conjugacy. We
see that the I-posterior variance σ 2

I1 does not depend on the
T-posterior covariance matrix �T 1, i.e. the uncertainty from
the surrogate training step is neglected.

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

p(ωI | yI , u = E-Log-Lik)

∝ p(ωI ) exp

{∫

log(p(yI | ωI , c))p(c | DT )dc

}

= N (ωI | μI0, σ
2
I0)

× exp

{∫

log(N (yI | ω̂

I c, σ

2
I ))N (c | μT 1, �T 1)dc

}

∝ N (ωI | μI1, σ
2
I1),

(38)

with

σ 2
I1 =(σ−2

I0 + σ−2
I (μ

(2)
T 1μ

(2)
T 1 + �

(2,2)
T 1 ))−1,

μI1 =σ 2
I1(σ

−2
I0 μI0

+ σ−2
I (μ

(2)
T 1yI − �

(1,2)
T1 − μ

(1)
T 1μ

(2)
T 1))

(39)

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 σ 2

I1 we see that E-Log-Lik produces counter-

intuitive results, as �
(2,2)
T 1 and σ 2

I1 are reciprocally related.
That is, as the surrogate model gets more uncertain in the
T-step, the I-posterior gets more certain.

123



Statistics and Computing (2025) 35 :66 Page 11 of 28 66

The E-Lik I-posterior is computed as

p(ωI | yI , u = E-Lik)

∝ p(ωI )

∫

p(yI | ωI , c)p(c | DT )dc

= N (ωI | μI0, σ
2
I0)

×
∫

N (yI | ω̂

I c, σ

2
I )N (c | μT 1, �T 1)dc

= N (ωI | μI0, σ
2
I0)

× N (yI | ω̂

I μT 1, ω̂



I �T1ω̂I + σ 2

I ).

(40)

Here,we cannot apply the normal-normal conjugacy, because
ω̂I is present in both the mean and the variance in the second
term of the product of the two normals. Instead, the result-
ing distribution is non-analytic and we have to use numerical
integration to calculate the normalization constant. Look-
ing 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 �T1.

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

p(ωI | yI , u = E-Post)

=p(ωI )

∫

p(yI | ωI , c)p(c | DT )
∫

p(yI | ωI , c)p(ωI )ddωI
dc

=N (ωI | μI0, σ
2
I0)

×
∫ N (yI | ω̂


I c, σ
2
I )N (c | μT 1, �T 1)

N (yI | c1 + c2μI0, c
cσ 2
I0 + σ 2

I )
dc.

(41)

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

In Table 2, we provide an overview over the chosen train-
ing data, measurement data, and hyperparameters for the T
and I-Step. We specifically vary the surrogate approximation
error parameter σ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 becomemore uncertain aswe increase the uncertainty
in the T-step, since our surrogate model gets less trustwor-
thy. 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. (39), as noted above.

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

y = M(ω) = 2

1 + exp(−10ω)
− 1 (42)

as the (true) simulator since it is invertible and smooth every-
where. 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

˜M(ω; c) = α

1 + exp(−β(ω − γ ))
+ δ, (43)

with surrogate approximation parameters c = (α, β, γ, δ).
Here, the true simulator is contained in the set of surrogate
models (for α = 2, β = 10, γ = 0, δ = −1).

For the T-Step, we generate the training set by setting the
input points ω

(i)
T ∈ [−1, 1] to the first NT 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(i)

T for i ∈ {1, . . . , NT }
are sampled from a normal distribution with mean equal to
the evaluated logistic simulator M at the input points and
standard deviation σS , i.e. y

(i)
T ∼ N (M(ω

(i)
T ), σ 2

S ). To avoid
sampling issues during the training step, we induce a small
simulation noiseσS = 0.01,which is not explicitlymodelled.

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

p(yT | c) = N (yT | ˜M(ωT , c), σ 2
A) (44)

with the surrogate approximation hyperparameter σ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

123



66 Page 12 of 28 Statistics and Computing (2025) 35 :66

Table 2 Parameters with
realized values for case study 1

Parameters T-Step I-Step
[

ω
(1)
T

ω
(2)
T

]

[

a
b

]

μT 0 σT 0 σA ω∗
I μI0 σI0 σI

Values

[−0.9
−0.3

] [

0.5
2

] [

0
0

]

10 {0.1, 0.5, 1} −0.5 0 1 0.1

Fig. 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 σA = {0.1, 0.5, 1} is varied

2014), an adaptive form of HamiltonianMonte Carlo (HMC)
(Neal 2011) which is a gradient based MCMC sampler via
the probabilistic programming language Stan (StanDevelop-
ment Team 2024).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 assessed
convergence using standard convergence checks (i.e., the R-
hat diagnostic of all parameters (Vehtari et al. 2021)).

For the I-Step,we generate NI = 5 randommeasurements
yI ∼ N (M(ω∗

I ), σ
2
I ), based on the simulator, true input

parameters ω∗
I , and the measurement error σI = 0.01. As

priors we set ωI ∼ N[−1,1](0, 0.52) (with truncation bounds
[−1, 1]) and σI ∼ uniform(0, 0.05). These hyperparame-
ters were chosen so that the true simulator could make valid
inference about the input parameters given the measurement
data.

In contrast to case study 1, we propagate the T-posterior
through samples and hence use the MC approximation (see
Sect. 2.2) for each of the four methods (Point, E-Lik, E-
Log-Lik, E-Post). We sample from the I-posterior of ωI and
σI using NUTS 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 epis-
temic 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
Sect. 2.2.5) with L = 25 clusters. As the simulator is con-

tained in the class of surrogates, no approximation error is
present (σA = 0) and hence, there is no aleatoric uncertainty
to propagate.

Posterior distributionsWe set σS = 0.01 and vary the num-
ber of training points NT = {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 NT leads to
smaller T-epistemic uncertainty. In Appendix Fig. 11 we
depict the pairs plot of the T-posterior draws of the logis-
tic surrogate.

We compare the I-posteriors resulting from the four
different methods using the measurements resulting from
exemplary underlying true inputs ω∗

I ∈ {−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 NT .
The I-posteriors of E-Lik and E-Post show similar behav-
ior 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 nar-
rower. For NT = 5, only E-Lik and E-Post have substantial
I-posterior probabilitymass on the true inputsω∗

I .Notably, all
uncertainty propagationmethods converge to the same results
as the epistemic uncertainty from the T-Step decreases.

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

123



Statistics and Computing (2025) 35 :66 Page 13 of 28 66

Fig. 4 Selected results for two-step procedure with the logistic sur-
rogate in case study 2. Left: For NT = {5, 7, 10} the training data
set DT (black dots) and the mean of the T-posterior predictive dis-
tribution (red lines) is shown. Right: For each underlying true input

ω∗
I ∈ {−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)

et al. 2023). Concretely, we perform the adapted SBC pro-
cedure for surrogate-based inference as detailed in Sect. 2.3.
For the I-Step trials, we simulated the true values of the inputs
ω∗
I and measurement error σ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 graphi-
cally, we show the empirical cumulative distribution function
(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 NT = {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 NT = 8 all methods show good calibration.

As an additional continuous calibration metric, we cal-
culate the log(γ )-statistic (Modrák et al. 2022) from our
SBC results. We vary NT = {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 (NT = 5) slightly miscalibrated, for
NT > 5 proper calibration is achieved. In contrast, E-Log-
Lik and Point are both miscalibrated for NT < 7, but with

more trainingdata, as theT-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 ω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.

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 surro-
gate 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):

˜M(ω; c) =
d
∑

i=0

ciψi (ω), (45)

with the vector of surrogate coefficients c = (c0, ..., cd), the
maximum degree of polynomials d, and Legendre polyno-
mials ψi (ω) (see Sudret (2008) for a detailed definition). We
consider wide, independent normal priors for all surrogate
coefficients: ci ∼ N (0, 5). In the following, we fix the max-

123



66 Page 14 of 28 Statistics and Computing (2025) 35 :66

Fig. 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 ω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 distri-

bution function (ECDF) for two different number of simulation points
NT = {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
NT = {5, 6, 7, 8, 9, 10}

imum polynomial degree to d = 5. By doing so, we create a
scenario in which our surrogate is unable to fit the underlying
truemodel appropriately such that an approximation error eA
is induced. This creates a scenario in which it is important to
propagate both T-epistemic and T-aleatoric uncertainty.

Posterior distributionsWe perform the T- and I-Step, as pre-
viously described in Sect. 3.2.1. During the T-Step, we vary
the number of training points NT = {10, 20} to control for
the T-epistemic uncertainty. For the I-Step, we set the true
input parameters to ω∗

I = −0.05 and set σS = 0.01 to avoid
sampling issues. We show the result of the T-Step in the

123



Statistics and Computing (2025) 35 :66 Page 15 of 28 66

Fig. 6 Selected results for two-step procedure with PCE surrogate in
case study 2. Left: For NT = {10, 20} the training data set DT (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 ω∗

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 θ = c and in the right column we propagate both T-
epistemic and T-aleatoric uncertainty via θ = (c, σA)

left column of Fig. 6, where we present both the T-posterior
predictive distribution (T-aleatoric and T-epistemic uncer-
tainty) and themean of the T-posterior predictive distribution
(T-epistemic uncertainty only). The T-epistemic uncertainty
becomes smaller with increasing number of training data
points NT , but the T-aleatoric uncertainty stays approxi-
mately 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. θ = c.
For NT = 10, the Point I-posterior and the E-Log-Lik I-
posterior are narrow despite high T-epistemic uncertainty. In
contrast, E-Lik andE-Post producewider 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 prop-
agating both T-epistemic and T-aleatoric uncertainty by also
propagating posterior draws of the surrogate approximation
error: θ = (c, σA) (seeSect. 2.1.2). In general, all I-posteriors
now tend to be wider than for θ = c and are more similar
to each other. However, in the presence of substantial T-
epistemic uncertainty (as for NT = 10), E-Post produces the
widest I-posterior, followed by E-Lik, Point, and E-Log-Lik.
In Appendix C.3.1 we provide further results for different
true parameter values.
CalibrationWe performed SBC for our PCE surrogate setup
and show the results in Fig. 7. The top part shows the ECDF-
difference plots for NT = 10 for two cases: θ = c and
θ = (c, σA). When only T-epistemic uncertainty is propa-

gated, 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 calibra-
tions 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(γ )-statistic of SBC under
varying NT = {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 σA is modelled as
constant over ωI , we cannot achieve proper calibration with
either method as epistemic uncertainty reduces.

3.3 Case study 3: uncertainty propagation in an SIR
model

Finally, we apply our two-step procedure to a real-world case
study in epidemiology. The SIR model and its variants are
often used to mathematically describe the spread of infec-
tious diseases (Hethcote 2000; Giordano et al. 2020). By
considering this model, we demonstrate the applicability of
themethod to complex real-world problems that lack analytic
solutions and require computationally expensive numerical
methods. Our approach is particularly useful in such scenar-
ios, as it allows to replace the complex simulation model
while propagating relevant surrogate uncertainties.
Setup The SIR simulation model M is defined through the
following system of differential equations:

123



66 Page 16 of 28 Statistics and Computing (2025) 35 :66

Fig. 7 Calibration and sharpness of the I-posteriors using PCE surro-
gate in case study 2. Top: ECDF difference plots for the I-posterior
distributions of ωI resulting from the four different methods. We set the
number training points NT = 10 and show the results for T-epistemic

uncertainty propagation (θ = c) and T-epistemic and T-aleatoric uncer-
tainty propagation (θ = (c, σA)). Center: log-gamma-statistics of SBC.
Bottom: sharpness (90 % CI) of I-posterior for four different I-Steps
(colored dots/lines) for NT = {10, 20, 30, 40, 50}

dS(t)

dt
= −βS(t)

I (t)

N
dI (t)

dt
= βS(t)

I (t)

N
− γ I (t)

dR(t)

dt
= γ I (t),

(46)

where S(t) describes the number of susceptibles, I (t) the
number of infectives, and R(t) the number of recovered
individuals at time t . Furthermore, β describes the constant
contact rate, γ the constant recovery rate, and N denotes
the constant population. In the following, we set the con-
stant population to N = 763, and fix the initial conditions
to the exemplary values I0 = 1, S0 = N − I0, R0 = 0. To

123



Statistics and Computing (2025) 35 :66 Page 17 of 28 66

solve the differential equation defined in Eq. (46), we use
the Dormand-Prince algorithm (Dormand and Prince 1980),
a 4th/5th order Runge–Kutta method as implemented in Stan
(Stan Development Team 2024).

Typically, measurement data is given for the number of
infected individuals, i.e. y = I (t). The unknown parameters
to be inferred are ω = (β, γ ). As a surrogate model, we con-
sider again a PCE (see Sect. 3.2.2), this timewithmultivariate
Legendre polynomials for the 3-dimensional input (t, β, γ ).
The maximum degree of the polynomials is set to 4, which
leads to 34 unknown coefficients c by the standard truncation
scheme (see Sudret (2008)). This setup is chosen to demon-
strate the applicability of the method to arbitrary complex
simulation and surrogate models, as long as samples can be
drawn from the posteriors.

For the T-Step, we generate the training data set using a
3-dimensional Sobol sequence Sobol’ (1967) with NT = 38
input points {(t (i), β(i), γ (i))}NT

i=1. The bounds of the input
parameters are t ∈ [1, 14], β ∈ [1, 3], and γ ∈ [0.1, 0.9].
We scale all input parameters linearly to [−1, 1], i.e. the
standard scaling of the Legendre polynomials. The output
y(i)
T = I (t (i)) is then obtained for a given time t (i), β(i),

and γ (i) by solving Eq. (46). This creates a setup with a low
simulation budget, leading to a high T-epistemic uncertainty.

The surrogate parameter priors are chosen in the same
way as in 3.2.2. To learn the simulationmodel as efficiently as
possible while still enforcing the constraint of a non-negative
infection count, we choose a log-normal T-/I-likelihood. The
T-model is fitted using NUTS with 4 chains of 1000 warmup
and 250 post-warmup sampling iterations, resulting in a total
of S = 1000 samples to propagate.

For the I-Step, we generate NI = 50 measurements by
sampling output responses y(i)

T ∼ NegBin(I (t), φ) given
evenly spaced t , true input parameters ω∗

I , and the shape
parameter φ = 9.6. As priors we set β ∼ N[1,3](2, 0.52)
and γ ∼ N[0.1,0.9](0.5, 0.252) in order to stay within the
parameter bounds that were used to train the surrogate model
and further set σI ∼ Half-Normal(0.52). The densities and
parameterizations of the used distributions are provided in
A.1.

In this setup, we propagate epistemic uncertainty con-
tained in c through the S = 1000 T-posterior samples. For
the Point, E-Lik, and E-Log-Lik methods, we sample from
the I-posterior ofωI and σI usingNUTSwith 32 chains, each
running for 1000 warmup and 1000 post-warmup iterations,
resulting in a total of 32000 I-posterior draws. Instead, for
E-Post, where we fit S = 1000 separate models (one for each
T-posterior draw), we use NUTS with 4 chains, each running
for 250 iterations, resulting in a total of 106 posterior draws.
ResultsWe set the underlying true input parameters to exem-
plary values ω∗

I = (ω∗
I ,1, ω

∗
I ,2) = (β∗, γ ∗) = (1.6, 0.4)

and compare the I-posteriors resulting from the four differ-
ent methods in Fig. 8. In the marginal I-posterior distribution

plots, we observe that both the Point and E-Log-Lik methods
produce an I-posterior that is overconfident in one dimension.
We also observe that both the E-Lik and E-Post methods
have substantial probability mass around the true values. In
the scatter plots showing the joint I-posterior distribution,
we observe that Point and E-Log-Lik produce samples that
do not cover the true value ω∗

I , while E-Post and E-Lik do.
The large uncertainty produced by these methods (E-Post
and E-Lik) resembles the uncertainty in the T-posterior, as
the surrogate was only trained on very limited data.

4 Conclusion

We introduced a general two-step procedure for surrogate-
based Bayesian inference. Within our approach, we prop-
agate 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 propagationmeth-
ods developedwithin 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 sim-
ulator. To evaluate the uncertainty calibration of the resulting
inference, we proposed an extension of simulation-based cal-
ibration suitable for our two-step surrogate approach.

As we demonstrate in our case studies, even in seem-
ingly simple setups, complex behavior occurs in terms of
posterior shape and calibration when propagating surrogate
uncertainty. In particular two uncertainty propagation meth-
ods (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 com-
plete 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 uncer-
tainty 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 stud-
ies, the produced inference was very similar. They differ in
the computational requirements and ease of parallelization
though (see Sect. 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). In our case studies, we focused on simula-
tors with up to three-dimensional input parameters in order
to simplify the presentation and establish a better intuition

123



66 Page 18 of 28 Statistics and Computing (2025) 35 :66

Fig. 8 Pairs plot of the I-posteriors using a PCE surrogate in case study
3. On the diagonals we depict for the true input ω∗

I = (β∗, γ ∗) =
(1.6, 0.4) (black vertical line) the marginal I-posterior distributions for

the Point, E-Lik, E-Post, and E-Log-Lik methods (colored lines). The
off-diagonals show the scatter plots of 5000 randomly subsampled I-
posterior draws. We propagate T-epistemic uncertainty via θ = c

about the overall approach. To study the applicability of our
methods on high-dimensional 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 approx-
imation 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 condition-
ing it on the input parameters, for example, as suggested in
Kohlhaas et al. (2023). Thiswill 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 pro-
hibitively so, as the number of (surrogate) 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 et al. 1993;
Graves 2011) and partial UQ on the last neural network layer
(Kristiadi et al. 2020; Fiedler and Lucia 2023; Harrison et al.
2024) could be sensible alternatives. Our uncertainty prop-
agation methods only require the ability to draw samples

from an (approximate) posterior in the training step that can
be subsequently passed to the inference step. This flexibil-
ity ensures that our approach is still applicable in scenarios
where a full BayesianUQ is infeasible, but the implication on
inference validity, in particular uncertainty calibration, need
to be further studied.

Appendix A. Notation

• ω: input
• y: output
• M: simulator
• ˜M: surrogate model
• c: surrogate approximation parameters
• d: number of surrogate parameters
• T-Step: surrogate training step

– NT : number of simulation training pairs
– ωT : simulation input
– yT : simulation output
– eS : simulation noise
– σS : simulation noise hyperparamters
– eA: approximation error
– σA: surrogate approximation error hyperparameters
– θ = (c, σA): trainable surrogate parameters
– S: number of T-posterior samples

• I-Step: surrogate inference step

123



Statistics and Computing (2025) 35 :66 Page 19 of 28 66

– NI : number of measurements
– ωI : quantity of interest (QoI)
– yI : measurement data
– eI : measurement noise
– σI : measurement noise hyperparameters
– K : number of I-posterior samples

A.1 Negative binomial, log-normal, and half-normal
distributions

Below, we show the Negative binomial, Log-normal and
Half-normal distributions. The probability mass function of
the Negative binomial distribution for scalar count n ∈ N

with the two positive parameters μ ∈ R
+ and φ ∈ R

+ is
given by:

pNegBinom(n | μ, φ) =
(

n + φ − 1

n

)(

μ

μ + φ

)n (
φ

μ + φ

)φ

.
(A1)

In this parameterization the mean and variance are given by

E[n] = μ and Var[n] = μ + μ2

φ
. (A2)

The probability density function of the Log-normal distribu-
tion for a positive scalar y ∈ R

+ with the parameters μ ∈ R

and σ ∈ R
+ is given by:

pLog-Normal(y | μ, σ)

= 1√
2πσ

1

y
exp

(

−1

2

(

log y − μ

σ

)2
)

.
(A3)

In this parameterization the mean and variance are given by

E[y] = exp

(

μ + σ 2

2

)

and (A4)

Var[y] = [exp(σ 2) − 1] exp(2μ + σ 2). (A5)

The probability density function of the Half-normal distribu-
tion for a positive scalar y ∈ R

+ with the parameter σ ∈ R
+

is given by:

pHalf-Normal(y | σ 2) =
√
2√

πσ
exp

(

− y2

2σ 2

)

. (A6)

Appendix B: Derivations

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

In the following, we show the inequality of E-Post and E-
Lik (Main Sect. 2.2) using a counterexample with discrete
randomvariables.Wedenote theT-posterior as p(θ) and omit
the dependence on the training dataDT to simplify notation.

Let ω ∈ {0, 1}, y ∈ {0, 1}, θ ∈ {0, 1}. We set the proba-
bility values

p(ω = 0) 1/2
p(ω = 1) 1/2

p(θ = 0) 1/2
p(θ = 1) 1/2

ω = 0, θ = 0 ω = 1, θ = 0
p(y = 0 | ω, θ) 1/4 1/2
p(y = 1 | ω, θ) 3/4 1/2

ω = 0, θ = 1 ω = 1, θ = 1
p(y = 0 | ω, θ) 1/2 1/2
p(y = 1 | ω, θ) 1/2 1/2

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

p(ω | y, u = E-Lik)

= p(ω)
∑

i p(y | ω, θ = i)p(θ = i)
∑

i
∑

j p(y | ω = j, θ = i)p(ω = j)p(θ = i)

p(ω | y, u = E-Post)

=
∑

i

p(ω)p(y | ω, θ = i)p(θ = i)
∑

j p(y | ω = j, θ = i)p(ω = j)

× p(ω = 0)p(y = 0 | ω = 0, θ = 1)p(θ = 1)
∑

j p(y | ω = j, θ = 1)p(ω = j)

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

p(ω = 0 | y = 0, u = E-Post)

=
∑

i

p(ω = 0)p(y = 0 | ω = 0, θ = i)p(θ = i)
∑

j p(y = 0 | ω = j, θ = i)p(ω = j)

= 1/2 · 1/4 · 1/2
3/8

+ 1/2 · 1/2 · 1/2
1/2

= 1/16

3/8
+ 1/8

1/2

= 5

12
≈ 0.417

Next, we calculate E-Lik for ω = 0 and y = 0:

p(ω = 0 | y = 0, u = E-Lik)

= p(ω = 0)
∑

i p(y = 0 | ω = 0, θ = i)p(θ = i)
∑

i
∑

j p(y = 0 | ω = j, θ = i)p(ω = j)p(θ = i)

123



66 Page 20 of 28 Statistics and Computing (2025) 35 :66

= 1/2 · (1/4 · 1/2 + 1/2 · 1/2)
7/16

= 3

7
≈ 0.429.

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

The normalization constants for E-Post were given by:

∑

j

p(y | ω = j, θ = 0)p(ω = j)

= p(y = 0 | ω = 0, θ = 0)p(ω = 0)

+ p(y = 0 | ω = 1, θ = 0)p(ω = 1)

= 1/4 · 1/2 + 1/2 · 1/2 = 3/8

and

∑

j

p(y | ω = j, θ = 1)p(ω = j)

= p(y = 0 | ω = 0, θ = 1)p(ω = 0)

+ p(y = 0 | ω = 1, θ = 1)p(ω = 1)

= 1/2 · 1/2 + 1/2 · 1/2 = 1/2.

The normalization constant for E-Lik was given by:

pnorm(y = 0)

=
∑

i

∑

j

p(y = 0 | ω = j, θ = i)p(ω = j)p(θ = i)

=p(y = 0 | ω = 0, θ = 0)p(ω = 0)p(θ = 0)

+ p(y = 0 | ω = 1, θ = 0)p(ω = 1)p(θ = 0)

+ p(y = 0 | ω = 0, θ = 1)p(ω = 0)p(θ = 1)

+ p(y = 0 | ω = 1, θ = 1)p(ω = 1)p(θ = 1)

=1/4 · 1/2 · 1/2 + 1/2 · 1/2 · 1/2
+ 1/2 · 1/2 · 1/2 + 1/2 · 1/2 · 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 (Main Sect. 2.2) are equivalent. The E-Log-Lik is
defined as:

log(p(ωI | yI , u = E-Log-Lik))

∝ log(p(ωI )) +
∫

log(p(yI | ωI , θ))p(θ | DT )dθ

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

log(p(ωI | yI , u = E-Log-Post))

=
∫

log(p(ωI | yI ))p(θ | yT )dθ

=
∫

[

log(p(ωI )) + log(p(yI | ωI , θ) − log(C(θ))
]

× p(θ | yT )dθ

= log(p(ωI ))

∫

p(θ | yT )dθ

+
∫

log(p(yI | ωI , θ))p(θ | yT )dθ

−
∫

log(C(θ))p(θ | yT )dθ

∝ log(p(ωI )) +
∫

log(p(yI | ωI , θ))p(θ | yT )dθ

In the last equationweused that the integral over a probability
distribution is one, i.e.

∫

p(θ | yT )dθ = 1 and since the
posterior is a function of ωI we can define a new constant
C1(θ) := − ∫

log(C(θ))p(θ | yT )dθ .

B.3 Case study 1: slope intercept model

In this sectionwe derive the E-Log-Lik for the slope intercept
model stated in Main Sect. 3.1. Let c = [c1, c2]T ,

�T =

⎡

⎢

⎢

⎢

⎢

⎣

1 ω
(1)
T

1 ω
(2)
T

...
...

1 ω
(NT )
T

⎤

⎥

⎥

⎥

⎥

⎦

,

ω̂I = [1, ωI ]T , μT 1 = [μ(1)
T 1, μ

(2)
T 2]T ,

�T1 =
[

�
(1,1)
T 1 �

(1,2)
T 1

�
(2,1)
T 1 �

(2,2)
T 1

]

,

B.3.1 E-Log-Lik

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

p(ωI | yI , u = E-Log-Lik)

∝ p(ωI ) exp

{∫

log(p(yI | ωI , c))p(c | DT )dc

}

= p(ωI ) exp

{∫

log(N (yI | ω̂

I c, σ

2
I ))

×N (c | μT 1, �T 1)dc}

= p(ωI ) exp

{

∫

− 1

2σ 2
I

(ω̂

I c − yI )

2

×N (c | μT 1, �T 1)dc}

= p(ωI ) exp

{

− 1

2σ 2
I

∫

(y2I − 2yI ω̂


I c + c
ω̂I ω̂



I c)

×N (c | μT 1, �T 1)dc}

123



Statistics and Computing (2025) 35 :66 Page 21 of 28 66

= p(ωI ) exp

{

− 1

2σ 2
I

(y2I − 2yI ω̂


I E[c]

+E[c
ω̂I ω̂


I c])

}

= p(ωI ) exp

{

− 1

2σ 2
I

(y2I − 2yI ω̂


I μT 1

+Tr(ω̂I ω̂


I �T 1) + μ


T 1ω̂I ω̂


I μT 1)

}

∝ p(ωI ) exp

{

− 1

2σ 2
I

(y2I − 2yI ω̂


I μT 1 + ω2

I�
(2,2)
T 1

+ ωI (�
(2,1)
T1 + �

(1,2)
T1 ) + μ

(2)
T 1μ

(2)
T 1ω

2
I

+2μ(1)
T 1μ

(2)
T 1ωI )

}

∝ p(ωI ) exp

{

− 1

2σ 2
I

(μ
(2)
T 1μ

(2)
T 1 + �

(2,2)
T 1 )

×(ωI − yμ(2)
T 1 − �

(1,2)
T 1 − μ

(1)
T 1μ

(2)
T 1

μ
(2)
T 1μ

(2)
T 1 + �

(2,2)
T1

)2

}

∝ N (ωI | μI1, σ
2
I1),

with

σ 2
I1 = (σ−2

I0 + σ−2
I (μ

(2)
T 1μ

(2)
T 1 + �

(2,2)
T 1 ))−1,

μI1 = σ 2
I1(σ

−2
I0 μI0 + σ−2

I (μ
(2)
T 1yI − �

(1,2)
T 1 − μ

(1)
T 1μ

(2)
T 1))

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

Tr(ω̂I ω̂


I �T1)

= ω2
I�

(2,2)
T 1 + ωI (�

(2,1)
T 1 + �

(1,2)
T 1 ) + �

(1,1)
T 1

and

μ

T 1ω̂I ω̂



I μT 1

= μ
(2)
T 1μ

(2)
T 1ω

2
I + 2μ(1)

T 1μ
(2)
T 1ωI + μ

(1)
T 1μ

(1)
T 1

Appendix C. Further results

C.1 Case study 1: slope only model

In Main Sect. 3.1 we showed the I-posteriors for the lin-
ear model with a slope and an intercept. In this section, we
consider an even simpler setup, where the simulator and sur-
rogate model are both a slope only model:

y = M(ω) = ˜M(ω, c) = c · ω, (C7)

where ω ∈ R is the input, y ∈ R the output and c ∈ R is
the unknown parameter. We consider to propagate only the

Fig. 9 Two-Step-Procedure for linear model with normal likelihoods
and normal priors

surrogate approximation parameter, i.e. we set θ = c. We set
normal priors and normal likelihoods. T-Step As our simu-
lation data set we use one input–output pair DT = {ωT , yT }
and NT = 1. We set a normal prior on c with mean μT 0

and variance σ 2
T 0. We calculate the posterior using the con-

jugate prior relation for a normal-normal model (Murphy
2007, 2012) (Fig. 9):

p(c | yT ) = p(c)p(yT | c)
p(yT )

(C8)

= N (c | μT 1, σ
2
T 1),

σ 2
T 1 = (σ−2

T 0 + σ−2
A ω2

T )−1,

μT 1 = σ 2
T 1(σ

−2
T 0 μT 0 + σ−2

A ωT yT )

(C9)

I-Step: Point
We set a normal prior on ωI :

p(ωI ) = N (ωI | μI0, σ
2
I0). (C10)

We use the following normal likelihood:

p(yI | ωI , c) = N (yI | c · ωI , σ
2
I ) (C11)

We compute the mean of the T-Step posterior: c̄ = μT 1 Now
we can compute the Point I-posterior using the conjugate
prior relation (Murphy 2012):

p(ωI | yI , u = Point)

∝ p(ωI )p(yI | ωI , c̄)

= N (ωI | μI0, σ
2
I0)N (yI | μT 1 · ωI , σ

2
I )

= N (ωI | μI1, σ
2
I1)

, with:

σ 2
I1 = (σ−2

I0 + σ−2
I μ2

T 1)
−1,

μI1 = σ 2
I1(σ

−2
I0 μI0 + σ−2

I μT 1yI )

123



66 Page 22 of 28 Statistics and Computing (2025) 35 :66

I-Step: E-Lik Derivation

p(ωI | yI , u = E-Lik)

∝ p(ωI )

∫

p(yI | ωI , c)p(c | DT )dc

= N (ωI | μI0, σ
2
I0)

×
∫

N (yI | ωI · c, σ 2
I )N (c | μT 1, σ

2
T 1)dc

= N (ωI | μI0, σ
2
I0)N (yI | ωI · μT 1, ω

2
I · σ 2

T 1 + σ 2
I )

I-Step: E-Log-Lik

p(ωI | yI , u = E-Log-Lik)

∝ p(ωI ) exp

{∫

log(p(yI | ωI , c))p(c | DT )dc

}

= N (ωI | μI0, σ
2
I0) exp

{∫

log(N (yI | ωI · c, σ 2
I ))

×N (c | μT 1, σ
2
T 1)dc

}

∝ N (ωI | μI0, σ
2
I0) exp

{

∫

(

1

2σ 2
I

(yI − ωI c)
2

)

×N (c | μT 1, σ
2
T 1)dc

}

= N (ωI | μI1, σ
2
I1)

, with:

σ 2
I1 = (σ−2

I0 + σ−2
I (μ2

T 1 + σ 2
T 1))

−1,

μI1 = σ 2
I1(σ

−2
I0 μI0 + σ−2

I μT 1yI )

I-Step: E-Post

p(ωI | yI , u = E-Post)

= p(ωI )

∫

p(yI | ωI , c)p(c | DT )
∫

p(yI | ωI , c)p(ωI )dωI
dc

= N (ωI | μI0, σ
2
I0)

×
∫ N (yI | ωI · c, σ 2

I )N (c | μT 1, σ
2
T 1)

N (yI | c · μI0, c2 · σ 2
I0 + σ 2

I )
dc

In Fig. 10 we compare the four uncertainty propaga-
tion methods for exemplary parameters. To control for
T-epistemic uncertainty we vary σA. In general, we notice
similar I-posteriors to the slope-intercept model (see Main
Sect. 3.1). Additionally, for large σA we observe that both
E-Lik and E-Post produce bimodal I-posteriors.

C.1.1 Derivation E-Log-Lik

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

p(ωI | yI , u = E-Log-Lik)

∝ p(ωI ) exp

{∫

log(p(yI | ωI , c))p(c | DT )dc

}

∝ N (ωI | μI1, σ
2
I1)

, with:

σ 2
I1 = (σ−2

I0 + σ−2
I (μ2

T 1 + σ 2
T 1))

−1,

μI1 = σ 2
I1(σ

−2
I0 μI0 + σ−2

I μT 1yI )

where we used:

∫

log(p(yI | ωI , c))p(c | DT )dc

=
∫

log(N (yI | ωI c, σ
2
I ))N (c | μT 1, σT 1)dc

=
∫

− 1

2σ 2
I

(ωI c − yI )
2N (c | μT 1, σT 1)dc

= − 1

2σ 2
I

∫

(y2I − 2yIωI c + c2ω2
I )

× N (c | μT 1, σT 1)dc

= − 1

2σ 2
I

(y2I − 2yIωIE[c] + ωIE[c2])

= − 1

2σ 2
I

(y2I − 2yIωIμT 1 + μ2
T 1 + σ 2

T 1)

∝ − 1

2σ 2
I

(μ2
T 1 + σ 2

T 1)(ωI − μT 1yI
μ2
T 1 + σ 2

T 1

)2

C.2 Case study 2: logistic model

C.3 T-posterior distribution

In Fig. 11 show the T-posterior pairs plots using the logis-
tic surrogate model ((same setup as in Main Sect. 3.2.1) for
NT = 7.

C.3.1 PCE posterior distributions

We show I-posterior density plots using the PCE surro-
gate to approximate the logistic function (same setup as in
Main Sect. 3.2) for additional true input parameters ω∗

I =
{−0.05, 0.1, 0.4}. In Fig. 12 we propagate only epistemic
uncertainty and in Fig. 13 we propagate both epistemic and
aleatoric uncertainty.

123



Statistics and Computing (2025) 35 :66 Page 23 of 28 66

Fig. 10 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 σA

Fig. 11 Pairs plot of the T-posterior draws of the logistic surrogate in Case study 2 for NT = 7

C.3.2 Computational complexity

To evaluate the computational complexity of the differ-
ent UP methods, we estimate their sampling times using a
PCE surrogate approximating a logistic function (see Main
Sect. 3.2). For each UP method, we sample 4 chains with
1000 warmup iterations and 2000 total iterations. How-
ever, 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 polynomi-
als d ∈ {2, 5, 9, 25} and the number of clusters K ∈
{2, 5, 25, 100, 250, 500, 750, 1000}. We plot the times on a
logarithmic scale in Fig. 14.

123



66 Page 24 of 28 Statistics and Computing (2025) 35 :66

Fig. 12 Additional results for T-epistemic UP (θ = c) with PCE surro-
gate in Case study 2. Left: For NT = {10, 30, 50} the training data set
DT (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 ω∗
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)

Fig. 13 Additional results for T-epistemic and T-aleatoric UP (θ =
(c, σA)) with PCE surrogate in Case study 2. Left: For NT =
{10, 30, 50} the training data set DT (black dots), the T-posterior
predictive distribution and the mean of the T-posterior predictive dis-

tribution (dark red and red lines) is shown. Right: For the true inputs
ω∗
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)

123



Statistics and Computing (2025) 35 :66 Page 25 of 28 66

Fig. 14 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 ∈ {2, 5, 9, 25} and the number
of clusters K ∈ {2, 5, 25, 100, 250, 500, 750, 1000}

Fig. 15 Effect of parallelization on the warmup and sampling times.
The runtimes for E-Log-Lik, E-Post, and Point are shown on a logarith-
mic scale

C.3.3 Parallelization

To examine the effectiveness of parallelization, we measure
the runtime for different numbers of workers, comparing the
warmup and sampling times across the Point, E-Log-Lik, and
E-Post methods. We use a PCE surrogate to approximate the
logistic simulator, as described in Section 3.2.2. The Point
method serves as a baseline, does not propagate uncertainty,
and is not parallelized. To compute the I-posteriors using the
E-Log-Lik and the E-Post method, we propagate S = 500
T-posterior samples.

We repeat the experiment three times and report the mean
sampling time alongwith the standard deviation. The number
of workers is varied across {1, 4, 16}. For each uncertainty
propagation method, we run a single MCMC chain with

1000 warm-up and 1000 post-warmup iterations. Figure15
presents the results on a logarithmic scale.

This plot confirms that the runtime for both E-Log-Lik
and E-Post decreases as the number of workers increases.
Additionally, we observe that E-Post benefits more from par-
allelization than E-Log-Lik, which aligns with the expected
degree of parallelizability of these methods as discussed in
Section 2.2.6.Due to implementation constraints, we exclude
E-Lik from this analysis, as its use of the log-sum-exp
trick for numerical stability currently limits straightforward
parallelization. However, given the similarities in implemen-
tation between E-Lik and E-Log-Lik, we expect that the
conclusions drawn for E-Log-Lik also apply to E-Lik.

Acknowledgements Partially funded by Deutsche Forschungsgemein-
schaft (DFG, German Research Foundation) under Germany’s Excel-
lence Strategy EXC 2075 - 390740016 and DFG Project 500663361.
We acknowledge the support by the Stuttgart Center for Simulation
Science.

Author Contributions P.B. and A.G. conceived of the idea. P.R., A.G.
and P.B. developed the idea. P.R. designed and carried out the experi-
ments and case studies. P.R. implemented the majority of the software
code. J.A. helped with the derivations for case study 1. P.R. wrote the
manuscript and J.A., A.G. and P.B. reviewed and edited the manuscript.

Funding Open Access funding enabled and organized by Projekt
DEAL.

Open Access This article is licensed under a Creative Commons
Attribution 4.0 International License, which permits use, sharing, adap-
tation, distribution and reproduction in any medium or format, as
long as you give appropriate credit to the original author(s) and the
source, provide a link to the Creative Commons licence, and indi-
cate if changes were made. The images or other third party material
in this article are included in the article’s Creative Commons licence,
unless indicated otherwise in a credit line to the material. If material
is not included in the article’s Creative Commons licence and your
intended use is not permitted by statutory regulation or exceeds the

123



66 Page 26 of 28 Statistics and Computing (2025) 35 :66

permitted use, youwill need to obtain permission directly from the copy-
right holder. To view a copy of this licence, visit http://creativecomm
ons.org/licenses/by/4.0/.

References

Alden, K., Cosgrove, J., Coles,M., Timmis, J.: Using emulation to engi-
neer andunderstand simulations of biological systems. IEEE/ACM
Trans. Comput. Biol. Bioinf. 17(1), 302–315 (2020)

Austin, P.C., White, I.R., Lee, D.S., van Buuren, S.: Missing data in
clinical research:A tutorial onmultiple imputation.Can. J.Cardiol.
37(9), 1322–1331 (2021)

Bayarri,M.J., Berger, J.O., Liu, F.:Modularization inBayesian analysis,
with emphasis on analysis of computer models. Bayesian Anal.
4(1), 119–150 (2009)

Ben-Gal, I.: Bayesian Networks, (2008)
Blomstedt, P., Mesquita, D., Lintusaari, J., Sivula, T., Corander, J.,

Kaski, S.: Meta-analysis of Bayesian analyses, (2019)
Brandstetter, J., van den Berg, R., Welling, M., Gupta, J. K.: Clifford

neural layers for PDE modeling. In: The Eleventh International
Conference on Learning Representations, ICLR 2023, Kigali,
Rwanda, May 1-5, 2023. OpenReview.net, (2023)

Breiman, L.: Bagging predictors. Mach. Learn. 24, 123–140 (2004)
Bühlmann, P.: Discussion of Big Bayes Stories and BayesBag. Stat. Sci.

29(1), 91–94 (2014)
Bürkner, P.-C., Kröker, I., Oladyshkin, S., Nowak,W.: A fully Bayesian

sparse polynomial chaos expansion approach with joint priors on
the coefficients and global selection of terms. J. Comput. Phys.
488, 112210 (2023)

Bürkner, P.-C., Scholz, M., Radev, S.T.: Some models are useful, but
how do we know which ones? Towards a unified Bayesian model
taxonomy. Statistics Surveys 17, 216–310 (2023)

Bärkner, P.-C.: brms: An R package for Bayesian multilevel models
using Stan. J. Statistical Software 80(1), 1–28 (2017)

Bärkner, P.-C.: Advanced Bayesian Multilevel Modeling with the R
Package brms. R J. 10(1), 395–411 (2018)

Bürkner, P.-C., Gabry, J., Weber, S., Johnson, A., Modrak, M., Badr,
H. S.,Weber, F., Ben-Shachar, M. S., and Rabel, H.: Running brms
models with within-chain parallelization. (2023). URL https://
CRAN.R-Project.org/package=brms.Vignette included inRpack-
age brms, version 2.17.0

Carpenter, B., Gelman, A., Hoffman, M.D., Lee, D., Goodrich, B.,
Betancourt, M., Brubaker, M., Guo, J., Li, P., Riddell, A.: Stan:
A probabilistic programming language. J. Stat. Softw. 76(1), 1–32
(2017)

Cleary, E., Garbuno-Inigo, A., Lan, S., Schneider, T., Stuart, A.M.:
Calibrate, emulate, sample. J. Comput. Phys. 424, 109716 (2021)

Cook, S.R., Gelman, A., Rubin, D.B.: Validation of software for
Bayesianmodels using posterior quantiles. J. Comput. Graph. Stat.
15(3), 675–692 (2006)

Dempster, A.P., Laird, N.M., Rubin, D.B.: Maximum likelihood from
incomplete data via the em algorithm. J. Royal Statist. Soc. Series
B (Methodological) 39(1), 1–38 (1977)

Dormand, J., Prince, P.: A family of embedded runge-kutta formulae.
J. Comput. Appl. Math. 6(1), 19–26 (1980)

Douady, C.J., Delsuc, F., Boucher, Y., Doolittle, W.F., Douzery, E.J.P.:
Comparison of Bayesian and maximum likelihood bootstrap mea-
sures of phylogenetic reliability. Mol. Biol. Evol. 20(2), 248–254
(2003)

Efron, B.: Bootstrap methods: another look at the jackknife. Ann. Stat.
7(1), 1–26 (1979)

Fiedler, F., Lucia, S.: Improved uncertainty quantification for neural net-
works with Bayesian last layer. IEEE Access 11, 123149–123160
(2023)

Gelman, A., Carlin, J., Stern, H., Dunson, D., Vehtari, A., and Rubin,
D.: Bayesian Data Analysis, Third Edition. Chapman&Hall/CRC
Texts in Statistical Science. Taylor & Francis, (2013)

Geyer, C. J.: Markov chain monte carlo maximum likelihood, (1991)
Giordano, G., Blanchini, F., Bruno, R., Colaneri, P., Filippo, A., Di

Matteo, A., Colaneri, M.: Modelling the covid-19 epidemic and
implementation of population-wide interventions in italy. Nat.
Med. 26, 1–6 (2020)

Gneiting, T., Balabdaoui, F., Raftery, A.E.: Probabilistic Forecasts, Cal-
ibration and Sharpness. J. R. Stat. Soc. Ser. B StatMethodol. 69(2),
243–268 (2007)

Goodfellow, I., Bengio, Y., Courville, A.: Deep learning. MIT Press,
Cambridge (2016)

Gorinova, M. I., Moore, D., Hoffman, M. D.: Automatic reparameteri-
sation of probabilistic programs, (2019)

Gramacy, R. B.: Surrogates: Gaussian Process Modeling, Design and
Optimization for the Applied Sciences. Chapman Hall/CRC, Boca
Raton, Florida. (2020). http://bobby.gramacy.com/surrogates/

Graves, A.: Practical variational inference for neural networks. In:
Shawe-Taylor, J., Zemel, R., Bartlett, P., Pereira, F., Weinberger,
K. (eds.)Advances inNeural Information Processing Systems, vol.
24. Curran Associates Inc (2011)

Gruber, C., Schenk, P. O., Schierholz, M., Kreuter, F., Kauermann, G.:
Sources of uncertainty in machine learning – a statisticians’ view,
(2023)

Halton, J.H.: On the efficiency of certain quasi-random sequences of
points in evaluating multi-dimensional integrals. Numer. Math.
2(1), 84–90 (1960)

Harrison, J., Willes, J., Snoek, J.: Variational Bayesian last layers.
CoRR, abs/2404.11599, (2024)

Hethcote, H.W.: The mathematics of infectious diseases. SIAM Rev.
42(4), 599–653 (2000)

Hinton, G. E., van Camp, D.: Keeping the neural networks simple by
minimizing the description length of the weights. In: Proceedings
of the Sixth Annual Conference on Computational Learning The-
ory, COLT ’93, page 5-13, New York, NY, USA. Association for
Computing Machinery, (1993)

Hoffman,M.D., Gelman, A.: The no-u-turn sampler: Adaptively setting
path lengths in hamiltonian monte carlo. J. Mach. Learn. Res.
15(47), 1593–1623 (2014)

Huggins, J. H.,Miller, J.W.: Robust inference andmodel criticismusing
bagged posteriors, (2020)

Hüllermeier, E., Waegeman, W.: Aleatoric and epistemic uncertainty in
machine learning: an introduction to concepts andmethods.Mach.
Learn. 110(3), 457–506 (2021)

Izmailov, P., Vikram, S., Hoffman, M. D., Wilson, A. G.: What are
bayesian neural network posteriors really like?. In: Meila, M. and
Zhang, T., editors, Proceedings of the 38th International Confer-
ence on Machine Learning, ICML 2021, 18-24 July 2021, Virtual
Event, volume 139 of Proceedings ofMachine Learning Research,
pp. 4629–4640. PMLR, (2021)

Jacob, P.E.,Murray,L.M.,Holmes,C.C.,Robert,C. P.:Better together?
statistical learning in models made of modules, (2017)

Jordan,M.I. (ed.): Learning in graphicalmodels.MITPress, Cambridge
(1999)

Kallioinen, N., Paananen, T., Bürkner, P.-C., Vehtari, A.: Detecting and
diagnosing prior and likelihood sensitivity with power-scaling,
(2022)

Kennedy,M.C.,O’Hagan,A.:Bayesian calibration of computermodels.
J. Royal Statist. Soc. Series B (Statistical Methodology) 63(3),
425–464 (2001)

Kim, S., Moon, H., Modrák, M., Säilynoja, T.: Simulation Based Cal-
ibration for rstan/cmdstanr models. (2023). https://hyunjimoon.
github.io/SBC/, https://github.com/hyunjimoon/SBC/

Kohlhaas, R., Kröker, I., Oladyshkin, S., Nowak, W.: Gaussian active
learning on multi-resolution arbitrary polynomial chaos emulator:

123

http://creativecommons.org/licenses/by/4.0/
http://creativecommons.org/licenses/by/4.0/
https://CRAN.R-Project.org/package=brms
https://CRAN.R-Project.org/package=brms
http://bobby.gramacy.com/surrogates/
https://hyunjimoon.github.io/SBC/
https://hyunjimoon.github.io/SBC/
https://github.com/hyunjimoon/SBC/


Statistics and Computing (2025) 35 :66 Page 27 of 28 66

concept for bias correction, assessment of surrogate reliability and
its application to the carbon dioxide benchmark. Comput. Geosci.
27, 1–21 (2023)

Kristiadi, A., Hein, M., Hennig, P.: Being bayesian, even just a bit, fixes
overconfidence in relu networks. In: Proceedings of the 37th Inter-
national Conference on Machine Learning, ICML’20. JMLR.org,
(2020)

Kucukelbir, A., Tran, D., Ranganath, R., Gelman, A., Blei, D.M.: Auto-
matic differentiation variational inference. J.Mach. Learn.Res. 18,
14:1-14:45 (2017)

Kuehnert, J., McGlynn, D., Remy, S. L., Walcott-Bryant, A., Jones,
A.: Surrogate Ensemble Forecasting for Dynamic Climate Impact
Models. (2022). arXiv:2204.05795 [physics]

Laloy, E., Rogiers, B., Vrugt, J., Mallants, D., Jacques, D.: Efficient
posterior exploration of a high-dimensional groundwater model
from two-stagemcmc simulation and polynomial chaos expansion.
Water Resources Research, 49 (2013)

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., Pfeffer, A.: Simulation intelligence: Towards a new
generation of scientific methods. CoRR, abs/2112.03235, (2021)

Li, J., Marzouk, Y.M.: Adaptive construction of surrogates for the
Bayesian solution of inverse problems. SIAM J. Sci. Comput.
36(3), A1163–A1186 (2014)

Li, Z., Kovachki, N. B., Azizzadenesheli, K., Liu, B., Bhattacharya,
K., Stuart, A. M., Anandkumar, A.: Fourier neural operator for
parametric partial differential equations. In: 9th International Con-
ference on Learning Representations, ICLR 2021, Virtual Event,
Austria, May 3-7, 2021. OpenReview.net, (2021)

Little, R.J., Rubin, D.B.: Statistical analysis withmissing data, vol. 793.
John Wiley & Sons (2019)

Lloyd, S.: Least squares quantization in pcm. IEEE Trans. Inf. Theory
28(2), 129–137 (1982)

MacQueen, J. B.: Some methods for classification and analysis of mul-
tivariate observations. In: Cam, L. M. L. and Neyman, J., (eds),
Proc. of the fifth Berkeley Symposium on Mathematical Statistics
and Probability (1967)

Martino, S., Riebler, A.: Integrated Nested Laplace Approximations
(INLA). arXiv e-prints, page arXiv:1907.01248 (2019)

Marzouk, Y., Xiu, D.: A stochastic collocation approach to Bayesian
inference in inverse problems, p. 6. NNSACenter for Prediction of
Reliability, Integrity and Survivability of Microsystems, PRISM
(2009)

Marzouk, Y.M., Najm, H.N., Rahn, L.A.: Stochastic spectral methods
for efficient Bayesian solution of inverse problems. J. Comput.
Phys. 224(2), 560–586 (2007)

McElreath, R.: Statistical Rethinking: A Bayesian Course with Exam-
ples in R and STAN. Chapman & Hall/CRC Texts in Statistical
Science. CRC Press, (2020)

Medina-Aguayo, F.J., Christen, J.A.: Penalised t-walk mcmc. J. Statis-
tical Plann. Inference 221, 230–247 (2022)

Meyer, L., Pottier, L., Ribés, A., Raffin, B.: Deep surrogate for direct
time fluid dynamics. CoRR, abs/2112.10296, (2021)

Mitra, E.D., Hlavacek, W.S.: Parameter estimation and uncertainty
quantification for systems biology models. Curr. Opinion Syst.
Biol. 18, 9–18 (2019)

Modrák,M.,Moon, A. H., Kim, S., Bärkner, P., Huurre, N., Faltejsková,
K., Gelman, A., Vehtari, A.: Simulation-based calibration check-
ing for bayesian computation: The choice of test quantities shapes
sensitivity, (2022)

Mohammadi, F., Kopmann, R., Guthke, A., Oladyshkin, S., Nowak,
W.: Bayesian selection of hydro-morphodynamic models under
computational time constraints. Adv. Water Resour. 117, 53–64
(2018)

Murphy,K. P.: ConjugateBayesian analysis of the gaussian distribution,
(2007)

Murphy, K.P.: Machine learning: a probabilistic perspective. The MIT
Press, Cambridge (2012)

Neal, R.: Handbook of markov chain monte carlo. Chapman and
Hall/CRC, Boca Raton (2011)

Oladyshkin, S., Nowak, W.: Data-driven uncertainty quantification
using the arbitrary polynomial chaos expansion. Reliab. Eng. Syst.
Safety 106, 179–190 (2012)

Petersen, K. B., Pedersen, M.S.: The matrix cookbook. Version
20121115, (2012)

Piironen, J., Paasiniemi, M., Vehtari, A.: Projective inference in high-
dimensional problems: Prediction and feature selection. Electronic
J. Statist. 14(1), 2155–2197 (2020)

Plummer, M.: Cuts in Bayesian graphical models. Stat. Comput. 25,
37–43 (2014)

Psaros,A.F.,Meng,X., Zou, Z.,Guo, L.,Karniadakis,G.E.:Unc