UQ[PY]LAB USER MANUAL

BAYESIAN INFERENCE FOR MODEL

CALIBRATION AND INVERSE PROBLEMS

P.-R. Wagner, J. Nagel, S. Marelli, B. Sudret

CHAIR OF RISK, SAFETY AND UNCERTAINTY QUANTIFICATION

STEFANO-FRANSCINI-PLATZ 5

CH-8093 Z

URICH

Risk, Safety &

Uncertainty Quantification

How to cite UQ[PY]LAB

C. Lataniotis, S. Marelli, B. Sudret, Uncertainty Quantiﬁcation in the cloud with UQCloud, Proceedings of the 4th

International Conference on Uncertainty Quantiﬁcation in Computational Sciences and Engineering (UNCECOMP

2021), Athens, Greece, June 27–30, 2021.

How to cite this manual

P.-R. Wagner, J. Nagel, S. Marelli, B. Sudret, UQ[py]Lab user manual – Bayesian inference for model calibration

and inverse problems, Report UQ[py]Lab -V1.0-113, Chair of Risk, Safety and Uncertainty Quantiﬁcation, ETH

Zurich, Switzerland, 2024

BIBT

X entry

@TechReport{UQdoc_10_113,

author = {Wagner, P.-R. and Nagel, J. and Marelli, S. and Sudret, B.},

title = {{UQ[py]Lab user manual -- Bayesian inversion for model calibration and

validation }},

institution = {Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich,

Switzerland},

year = {2024},

note = {Report UQ[py]Lab - V1.0-113}

}

List of contributors:

Name Contribution

A. Hlobilov

a Translation from the UQLab manual

Document Data Sheet

Document Ref. UQ[PY]LAB-V1.0-113

Title: UQ[PY]LAB user manual – Bayesian inference for model calibration

and inverse problems

Authors: P.-R. Wagner, J. Nagel, S. Marelli, B. Sudret

Chair of Risk, Safety and Uncertainty Quantiﬁcation, ETH Zurich,

Switzerland

Date: 27/05/2024

Doc. Version Date Comments

V1.0 27/05/2024 Initial release

Abstract

Bayesian inference is a powerful tool for probabilistic model calibration and inversion. It

provides a comprehensive framework for combining information about parameters prior to

observations with information obtained from experiments. In Bayesian inference, this com-

bined information is expressed in a so-called posterior distribution of the parameters.

The UQ[PY]LAB Bayesian inference module offers an easy way to setup a Bayesian in-

verse problem and to compute its posterior distribution. It makes use of other available

UQ[PY]LAB modules (UQ[PY]LAB User Manual – the INPUT module, UQ[PY]LAB User Man-

ual – the MODEL module) to deﬁne the forward model and the prior distribution. For the

computation of the posterior distribution, state-of-the-art algorithms are supplied. The man-

ual for the Bayesian inversion module is divided into three parts:

• A brief introduction to the main ideas and theoretical foundations of Bayesian inver-

sion and discussions on Markov chain Monte Carlo, spectral likelihood expansion, and

stochastic spectral embedding algorithms;

• An example-based guide with an explanation of the available options and methods;

• A comprehensive reference list detailing all available functionalities of the Bayesian

inversion module.

Keywords: UQ[PY]LAB, Bayesian inversion, model calibration, inverse problems, Markov

chain Monte Carlo, spectral likelihood expansion, stochastic spectral embedding

Contents

1 Theory 1

1.1 Bayesian inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Bayesian model calibration and inverse problems . . . . . . . . . . . . . . . . 4

1.2.1 A wide class of problems sharing the same methods . . . . . . . . . . . 4

1.2.2 Simple problems with known discrepancy parameters . . . . . . . . . 5

1.2.3 General case: discrepancy with unknown parameters . . . . . . . . . . 6

1.2.4 Multiple data groups . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.2.5 Inverse solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.2.6 Model predictions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.3 Markov chain Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.3.1 Metropolis–Hastings algorithm . . . . . . . . . . . . . . . . . . . . . . 10

1.3.2 Adaptive Metropolis algorithm . . . . . . . . . . . . . . . . . . . . . . 11

1.3.3 Hamiltonian Monte Carlo algorithm . . . . . . . . . . . . . . . . . . . 12

1.3.4 Afﬁne invariant ensemble algorithm . . . . . . . . . . . . . . . . . . . 13

1.3.5 Assessing convergence in MCMC simulations . . . . . . . . . . . . . . 14

1.4 Spectral likelihood expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

1.5 Stochastic spectral likelihood embedding . . . . . . . . . . . . . . . . . . . . . 19

1.5.1 Sequential partitioning algorithm . . . . . . . . . . . . . . . . . . . . . 20

2 Usage 23

2.1 Reference problem: calibration of a simply supported beam model . . . . . . 23

2.2 Problem setup and solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.2.1 Initialize UQ[PY]LAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.2.2 Specify a prior distribution . . . . . . . . . . . . . . . . . . . . . . . . 25

2.2.3 Create a forward model . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.2.4 Provide measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.2.5 Perform the Bayesian inverse analysis . . . . . . . . . . . . . . . . . . 26

2.2.6 Advanced options: discrepancy model . . . . . . . . . . . . . . . . . . 28

2.3 Multiple model outputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.3.1 Create a forward model . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.3.2 Provide measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.3.3 Perform the inverse analysis . . . . . . . . . . . . . . . . . . . . . . . . 32

2.3.4 Advanced options: discrepancy model . . . . . . . . . . . . . . . . . . 33

2.4 Advanced options: solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2.4.1 MCMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2.4.2 Spectral likelihood expansion . . . . . . . . . . . . . . . . . . . . . . . 44

2.4.3 Stochastic spectral likelihood embedding . . . . . . . . . . . . . . . . . 45

2.4.4 No solver: posterior point by point evaluation . . . . . . . . . . . . . . 46

2.5 Advanced feature: multiple forward models . . . . . . . . . . . . . . . . . . . 47

2.5.1 Specify a prior distribution . . . . . . . . . . . . . . . . . . . . . . . . 48

2.5.2 Create a forward model . . . . . . . . . . . . . . . . . . . . . . . . . . 48

2.5.3 Provide measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

2.5.4 Deﬁne a discrepancy model . . . . . . . . . . . . . . . . . . . . . . . . 50

2.5.5 Perform the inverse analysis . . . . . . . . . . . . . . . . . . . . . . . . 50

3 Reference List 55

3.1 Create a Bayesian inverse analysis . . . . . . . . . . . . . . . . . . . . . . . . . 58

3.1.1 Data dictionary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

3.1.2 Forward model dictionary . . . . . . . . . . . . . . . . . . . . . . . . . 59

3.1.3 Discrepancy model options . . . . . . . . . . . . . . . . . . . . . . . . 59

3.1.4 Solver options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.2 Accessing analysis results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

3.3 Post-processing results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

3.3.1 Markov chain Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . 66

3.3.2 Spectral likelihood expansions . . . . . . . . . . . . . . . . . . . . . . 69

3.3.3 Stochastic spectral likelihood embedding . . . . . . . . . . . . . . . . . 70

3.4 Printing/Visualizing results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

3.4.1 Printing the results: uq.print . . . . . . . . . . . . . . . . . . . . . . . 72

3.4.2 Graphically display the results: uq.display . . . . . . . . . . . . . . . . 72

Chapter 1

Theory

This section contains a short introduction to Bayesian methods (Gelman et al., 2014) with

a focus on inverse problems (Tarantola, 2005; Kaipio and Somersalo, 2005). An inverse

problem arises when unknown parameters that cannot be directly measured are estimated

based on experimental data that is only indirectly related to the parameters through a com-

putational model. The problem is called inverse, because instead of propagating information

about input parameters through a computational model (so-called forward approach), the

goal is to propagate information about the observations backwards to obtain insight on the

model inputs. This formulation encompasses a range of problems in the engineering and

natural sciences (Hadidi and Gucunski, 2008; Beck, 2010; Yuen and Kuok, 2011). This is a

very important extension to the Bayesia manual.

1.1 Bayesian inference

Statistics is generally described as the science that allows one to build models of com-

plex phenomena based on data. Statistical inference usually considers that this data X

def

, . . . , x

} is made of independent realizations of an underlying random vector with an

associated probability density function (PDF) π(·), whose properties have to be established

from that data. Parametric statistical models make an assumption on the shape of this PDF

(e.g. Gaussian, Weibull, lognormal in the one-dimensional case), and the goal of inference is

to estimate the parameters θ of this PDF given the data or more formally:

Θ|X ∼ π(θ|X ) (1.1)

where Θ is the random vector associated with the parameters θ and the vertical line | is used

to denote a conditional dependence of the quantity on the left on the quantity on the right.

When a sufﬁcient amount of data exists, classical estimators can be used: for instance, if a

Gaussian distribution X |θ ∼ N (x|µ, σ

) is to be ﬁtted to a sufﬁciently large data set X , the

UQ[PY]LAB user manual

empirical mean and standard deviation of the sample may be used as estimators

θ of the

parameters θ

def



µ, σ



. Such a direct estimation is, however, not reliable when there is

only a handful of data points. In technical terms, the statistical uncertainty of the estimator

θ becomes too large in this case.

In this context, Bayesian statistics allows one to ﬁt a statistical model by combining some prior

knowledge on the parameters with the (possibly few) observed data points, using Bayes’

theorem

. Before considering the data, in the Bayesian paradigm, the parameters of the

parametric distribution π(x|θ) are considered as a random vector denoted by Θ which is

assumed to follow the so-called prior distribution (with support D

Θ ∼ π(θ). (1.2)

This subjective choice should reﬂect the level of information existing on the parameters θ

before any measurement of X is carried out. From Bayes’ theorem, the posterior distribution

of the parameters, denoted as π(θ|x), is obtained by:

π(θ|x) =

π(x|θ) π(θ)

π(x)

. (1.3)

Consider now a data set of measured values X = {x

, . . . , x

}, whose points are viewed as

independent realizations of X |θ ∼ π(x|θ). With these measurements the likelihood function

L(θ; X ), a function of the parameters θ, can be deﬁned:

L : θ 7→ L(θ; X )

def

k=1

π(x

|θ). (1.4)

This implicitly assumes independence between individual measurements in X . Intuitively

the likelihood function for a given θ returns the relative likelihood of observing the data at

hand, under the assumption that it follows the prescribed parametric distribution π(x|θ).

Following Bayes’ theorem, the posterior distribution π(θ|X ) of the parameters θ given the

observations in X can now be written as:

π(θ|X ) =

L(θ; X ) π(θ)

, (1.5)

where the normalizing factor Z, known as the evidence or marginal likelihood, shall ensure

that this distribution integrates to 1:

def

L(θ; X ) π(θ) dθ. (1.6)

The posterior distribution in Eq. (1.5) summarizes the information inferred about the param-

Bayes’ theorem is an elementary result of probability theory that reads as follows: for two events A and B

with non-zero probabilities, the following equation holds: P (B |A) =

P (A|B) P (B)

P (A)

, where P (A|B) denotes

the conditional probability of A given B.

UQ[PY]LAB-V1.0-113 - 2 -

Bayesian inference for model calibration and inverse problems

eters by combining the prior knowledge and the observed data. In this sense, the posterior

π(θ|X ) is an “update” of the prior distribution π(θ).

The practical computation of posterior distributions is not trivial. Particular analytical solu-

tions exist only for speciﬁc combinations of prior distributions on Θ and likelihood functions,

the so-called conjugate distributions (Gelman et al., 2014). In the general case though, sam-

pling methods shall be used such as Markov chain Monte Carlo simulation (see Section 1.3

for details).

In practical applications related to uncertainty quantiﬁcation, the purpose of Bayesian infer-

ence may not just be to ﬁnd the posterior distribution π(θ|X ), but also to propose “the best

distribution” for the parameters X |X given the information and data at hand. One possi-

bility is to select a point estimator

θ, i.e. a particular value from the posterior distribution of

π(θ|X ). Then the point posterior distribution of X |X simply reads:

π(x|X )

def

= π(x|

θ). (1.7)

Popular choices for

θ are the posterior mean, which is the mean value of the posterior dis-

tribution (1.5) and the posterior mode, a.k.a. maximum a posteriori (MAP), which is the

mode of this posterior distribution. Such choices disregard the estimation uncertainty in the

parameters θ.

In contrast, it is also possible to incorporate the uncertainty in θ into the prior and posterior

assessment of X and X |X respectively. This results in the so-called predictive distributions.

The prior predictive distribution of X is obtained by “averaging” the parametric distribution

π(x|θ) over the prior distribution π(θ):

π(x)

def

π(x|θ) π(θ) dθ. (1.8)

The posterior predictive distribution π(x|X ) is obtained by “averaging” the parametric distri-

bution π(x|θ) over the posterior distribution π(θ|X ) in Eq. (1.5):

π(x|X )

def

π(x|θ) π(θ|X ) dθ. (1.9)

It is noted here that for a single measurement X = x

, the prior predictive distribution at X equals the

normalization constant Z deﬁned in Eq. (1.6).

UQ[PY]LAB-V1.0-113 - 3 -

UQ[PY]LAB user manual

1.2 Bayesian model calibration and inverse problems

1.2.1 A wide class of problems sharing the same methods

Let us consider a computational model M that allows the analyst to predict certain quantities

of interest gathered in a vector y ∈ R

out

as a function of input parameters x:

M : x ∈ D

⊂ R

7→ y = M(x) ∈ R

out

. (1.10)

Such models are commonly established based on ﬁrst principles in engineering sciences (e.g.

mechanics, electromagnetism, ﬂuid dynamics), but also natural sciences (e.g. geophysics,

wave propagation). Although analytical models with closed-form equations may be used,

the vast majority of computational models are black-box computer codes that solve the un-

derlying differential equations that govern the system of interest.

When input parameters {x

, i = 1, . . . , M} are not measurable directly, one resorts to mea-

suring the quantities of interest. Let us consider N independent measurement y

∈ R

out

gathered in a data set Y

def

= {y

, . . . , y

}. In the context of computational modelling and

uncertainty quantiﬁcation, two main classes of applications beneﬁt from Bayesian inference,

namely model calibration and inverse problems. The two classes are very much related to each

other and they share the same problem statement and solution techniques, but they differ in

their ﬁnal focus.

On the one hand, Bayesian inversion focuses on the identiﬁcation of the values of the input

parameters x, rather than on the model used to infer them. For this reason, it is also known

as Bayesian inverse modelling: instead of predicting a model response from a set of input

parameters, the latter are inferred from a set of observed model responses Y. This is the

typical usage scenario in tomographic imaging applications, where the goal is to identify

the set of input parameters that caused a speciﬁc set of observations. The resulting inferred

input parameters can then be directly used to identify anomalies (e.g. position and length of

cracks in pressure vessels), or to identify properties of interest (e.g. location and volume of

subsurface oil reservoirs).

On the other hand, Bayesian model calibration focuses on identifying the input parameters of

a computational model to allow one to recover the observations in Y . A common scenario in

this respect is identifying unknown properties of key components of a complex system, based

on their observed response to controlled external loads in a laboratory experiment. Through

this procedure, known as calibration, the inferred values (and possibly the uncertainty of the

estimation) can then be used to predict the response of the same system to different external

loads, or even to design different systems sharing the same calibrated model component. This

approach is at the basis of the so-called veriﬁcation and validation under uncertainty paradigm

(VVUQ) that is gaining momentum in the engineering practice worldwide (Oberkampf et al.,

2004; Oberkampf and Roy, 2010; Hu and Orient, 2016).

UQ[PY]LAB-V1.0-113 - 4 -

Bayesian inference for model calibration and inverse problems

1.2.2 Simple problems with known discrepancy parameters

Regardless of the speciﬁc context (model calibration or inversion), all Bayesian inverse prob-

lems share the same ingredients: a computational forward model M, a set of input parame-

ters x ∈ D

that need to be inferred, and a set of experimental data Y.

The forward model x 7→ M(x) is a mathematical representation of the system under con-

sideration. All models are always simpliﬁcations of the real world. Thus, to connect model

predictions to the observations Y, a discrepancy term shall be introduced. We consider the

following well-established format:

y = M(x) + ε, (1.11)

where ε ∈ R

out

is the term that describes the discrepancy between an experimental obser-

vation y and the model prediction. For the sake of simplicity, we consider it as an additive

Gaussian discrepancy

with zero mean value and given covariance matrix Σ in this introduc-

tion:

ε ∼ N (ε |0, Σ). (1.12)

This discrepancy term represents in practice the effects of measurement error (on y

∈ Y)

and model inaccuracy. In the above equation, again for the sake of simplicity, this term is

supposed to have a zero mean, but it could more generally include a model bias term.

In the context of model inversion or calibration, the goal is to ﬁnd the optimal values of the

input parameters x that allow one to ﬁt the model predictions to the observations. In this

respect the epistemic uncertainty (lack of knowledge) on the input parameters is modelled by

considering the input parameters as a random vector X ∼ π(x), with given prior distribution

as in Eq. (1.2).

Note: In this section, the measurement data is denoted by Y, which plays the role of X

in Section 1.1. In contrast, the parameters to infer are the input parameters X of

the computational model M, which plays the role of parameters Θ in Section 1.1.

From Eqs. (1.11), (1.12) a particular measurement point y

∈ Y is a realization of a Gaussian

distribution with mean value M(x) and covariance matrix Σ. This distribution is called the

error or discrepancy model π(y|x) given by:

π(y|x) = N (y|M(x), Σ). (1.13)

If N independent measurements y

are available and gathered in the data set Y = {y

, . . . , y

It is noted here that this simple Gaussian discrepancy assumption is only one out of many possible models. In

a more general setting, other distributions for the discrepancy are used as well (Schoups and Vrugt, 2010). Due

to the widespread use of the additive Gaussian models in engineering disciplines, the discussion is limited to this

discrepancy type.

UQ[PY]LAB-V1.0-113 - 5 -

UQ[PY]LAB user manual

the likelihood function can thus be written as

L(x; Y) =

i=1

N (y

|M(x), Σ)

i=1

(2π)

out

det(Σ)

exp



−



− M(x)



⊺

−1



− M(x)





(1.14)

Combining the prior π(x) and the above likelihood L(x; Y), the posterior distribution in

Eq. (1.5) establishes the solution of the inverse problem:

π(x|Y) =

π(x)

i=1

N (y

|M(x), Σ). (1.15)

It summarizes the collected information about the unknowns x after conditioning on the

data. In this sense, the data are inverted through the forward model M.

1.2.3 General case: discrepancy with unknown parameters

In many practical situations, it is unrealistic to assume that the residual covariance matrix

Σ in Eq. (1.12) is perfectly known. However, by parametrizing the matrix as Σ(x

), one

may treat its parameters x

as additional unknowns that can be inferred jointly with the

input parameters of M. In this setting the parameter vector is deﬁned by x = (x

, x

i.e. a combined vector of forward model parameters x

and discrepancy parameters x

For the sake of simplicity, consider a diagonal covariance matrix of the form Σ = σ

out

with unknown residual variances σ

= Var[ε

], i = 1, . . . , N

out

. This way, the discrepancy

parameter vector reduces to a single scalar, i.e. x

≡ σ

Assuming that one can elicit a prior distribution π(x

) for the unknown variance σ

, and by

treating the uncertain model and the discrepancy parameters as being priorly independent,

one gets the joint prior distribution

π(x) = π(x

)π(σ

). (1.16)

The likelihood is then given as

L(x

, σ

; Y) =

i=1

(2πσ

)

out

exp



−

2σ



− M(x

)



⊺



− M(x

)





. (1.17)

With the prior distribution in Eq. (1.16) and likelihood function in Eq. (1.17), the corre-

sponding posterior distribution can then again be computed as in Eq. (1.15):

π(x

, σ

|Y) =

π(x

)π(σ

)L(x

, σ

; Y). (1.18)

This posterior distribution summarizes the updated information about the unknowns (x

, x

≡

UQ[PY]LAB-V1.0-113 - 6 -

Bayesian inference for model calibration and inverse problems

) after conditioning on the data Y. One may then extract the marginals π(x

M,i

|Y) of indi-

vidual forward model inputs or the marginal π(x

|Y) ≡ π(σ

|Y) of the residual variance.

More generally, any parametrization Σ(x

) of the positive deﬁnite covariance matrix can

be incorporated into the Bayesian analysis. This just requires the speciﬁcation of a prior

distribution π(x

) and construction of a likelihood function of the more general form:

L(x

, x

; Y) =

i=1

(2π)

out

det(Σ(x

))

exp



−



− M(x

)



⊺

−1

)



− M(x

)





(1.19)

1.2.4 Multiple data groups

In practice, it occurs frequently that the measurements Y = {y

, . . . , y

} stem from various

measurement devices or experimental conditions with different discrepancy properties. In

these cases, it is necessary to arrange the elements of Y in disjoint data groups and deﬁne

different likelihood functions for each data group.

Denoting the g-th data group by G

(g)

= {y

}

i∈u

, where u ⊆ {1, . . . , N}, the full data set can

be combined by

Y =

[

g=1

(g)

. (1.20)

Each of the N

data groups contains measurements collected with the same instruments un-

der similar measurement conditions. With this, it is clear that each data group requires a

different likelihood function L

(g)

describing the experimental conditions that led to measur-

ing G

(g)

. Possible choices for L

(g)

are presented in Section 1.2.2 and Section 1.2.3. Under the

assumption of independence between the N

measurement conditions, the full likelihood

function can then be written as

L(x

, x

; Y) =

g=1

(g)

, x

(g)

; G

(g)

), (1.21)

where x

(g)

are the parameters of the g-th discrepancy group.

1.2.5 Inverse solution

The posterior distribution of the parameters computed by Eq. (1.15), is often characterized

through its ﬁrst statistical moments. The posterior mean vector is given as

E [X |Y] =

x π(x|Y) dx. (1.22)

UQ[PY]LAB-V1.0-113 - 7 -

UQ[PY]LAB user manual

It can be considered as a point estimate of the unknown parameter values. The estimation

uncertainty can be quantiﬁed through the posterior covariance matrix

Cov[X |Y] =

(x − E [X |Y])(x − E [X |Y])

⊺

π(x|Y) dx. (1.23)

One may also be interested in the posterior marginals. The univariate posterior marginal of

a speciﬁc parameter x

with i ∈ {1, . . . , M} can be computed by integration over the other

components (sometimes called nuisance parameters):

|Y) =

∼i

π(x|Y) dx

∼i

, (1.24)

where x

∼i

refers to the parameter vector x excluding the i-th parameter x

More generally, it is also possible to deﬁne multivariate posterior marginals by integrating

the posterior over all but the parameters of interest. To this end, we split the random vector

X into two vectors X

with components {X

}

i∈u

∈ D

and X

with components {X

}

i∈v

∈

, where u and v are two non-empty disjoint index sets such that u ∪ v = {1, . . . , M}. The

multivariate posterior marginals then read:

|Y) =

π(x|Y) dx

. (1.25)

In practical inverse problems, the posterior distribution π(x|Y) can also be an intermediate

quantity that is further used for computing the conditional expectation of a certain quantity of

interest (QoI) h: D

→ R. This can be anything from a simple analytical function to complex

secondary models. This conditional expectation is simply the expectation of h(X |Y) and is

computed by the integral

E [h(X |Y)] =

h(x) π(x|Y) dx. (1.26)

1.2.6 Model predictions

To assess the predictive capabilities of a computational model, the Bayesian inference frame-

work offers the possibility to compute predictive distributions, as seen in Section 1.1. Using

the previously deﬁned discrepancy model, the prior predictive distribution from Eq. (1.8) can

be written as

π(y) =

π(y|x)π(x) dx. (1.27)

UQ[PY]LAB-V1.0-113 - 8 -

Bayesian inference for model calibration and inverse problems

The posterior predictive distribution (see Eq. (1.9)) can be similarly written as

π(y|Y) =

π(y|x)π(x|Y) dx. (1.28)

Owing to the considered additive Gaussian discrepancy model in Eq. (1.13), a sample from

those predictive distributions can be obtained by using samples from X and X |Y respec-

tively, propagating them through the model M and adding an independently sampled dis-

crepancy term ε.

1.3 Markov chain Monte Carlo

Posterior distribution as in Eq. (1.15) do not have a closed-form solution in practice. One

widespread option to solve inverse problems relies upon Markov chain Monte Carlo (MCMC)

simulations (Robert and Casella, 2004; Liu, 2004).

The basic idea of MCMC simulations is to construct a Markov chain (X

(1)

, X

(2)

, . . .) over

the prior support D

with an invariant distribution that equals the posterior distribution of

interest. Markov chains can be uniquely deﬁned by their transition probability K(x

(t+1)

(t)

)

from the step x

(t)

of the chain at iteration t to the step x

(t+1)

at the subsequent iteration t+1.

Then, the posterior is the invariant distribution of the Markov chain if the speciﬁed transition

probability fulﬁls the so-called detailed balance condition:

π(x

(t)

|Y) K(x

(t+1)

(t)

) = π(x

(t+1)

|Y) K(x

(t)

(t+1)

). (1.29)

This condition ensures that the Markov chain is reversible, i.e., that the probability to be at

(t)

and move to x

(t+1)

is equal the probability to be at x

(t+1)

and move to x

(t)

. By integrating

this condition over dx

(t)

, it can be shown that the invariant distribution of the Markov chain

is the posterior distribution:

π(x

(t+1)

|Y) =

π(x

(t)

|Y) K(x

(t+1)

(t)

) dx

(t)

. (1.30)

A Markov chain constructed this way can be used to approximate the expectation in Eq. (1.26)

as the iteration average of the T + 1 generated sample points x

(t)

E [h(X)|Y] ≈

t=1

h(x

(t)

). (1.31)

A prototypical technique to ensure that this equation is fulﬁlled is the Metropolis–Hastings

(MH) algorithm (Metropolis et al., 1953; Hastings, 1970) that is based on proposing and

subsequently accepting or rejecting candidate points. In the following, the original MH al-

UQ[PY]LAB-V1.0-113 - 9 -

UQ[PY]LAB user manual

gorithm and three other popular MCMC algorithms are discussed as techniques to efﬁciently

sample from the posterior distribution.

1.3.1 Metropolis–Hastings algorithm

In the Metropolis–Hastings algorithm (MH), a chain is initialized at a certain seed point x

(0)

∈

from the admissible domain. At iteration t from the current point x

(t)

, one then draws a

candidate point x

(⋆)

from a proposal distribution p(x

(⋆)

(t)

). Subsequently, the candidate is

accepted (i.e., x

(t+1)

= x

(⋆)

) with probability:



(⋆)

, x

(t)



= min

(

π(x

(⋆)

|Y) p(x

(t)

(⋆)

)

π(x

(t)

|Y) p(x

(⋆)

(t)

)

, (1.32)

and rejected otherwise (i.e., x

(t+1)

= x

(t)

). With this procedure, the transition probability

fulﬁlls Eq. (1.30) and the chain of sample points eventually follows the posterior distribution.

These sample points can then, for example, be used to approximate expectations under the

posterior distribution as shown in Eq. (1.31).

It is advantageous that the model evidence cancels out from the acceptance probability in

Eq. (1.32). Hence, the MH algorithm only calls for pointwise evaluations of the unnormalized

posterior density π(x|Y) ∝ L(x; Y)π(x). This avoids the calculation of the often intractable

integral in Eq. (1.6).

The original Metropolis algorithm is based on a symmetrical proposal distribution with p(x

(⋆)

(t)

) =

p(x

(t)

(⋆)

). In this case, the acceptance probability in Eq. (1.32) reduces to:



(⋆)

, x

(t)



= min

(

π(x

(⋆)

|Y)

π(x

(t)

|Y)

)

. (1.33)

A commonly used symmetric proposal is the Gaussian distribution p(x|x

(t)

) = N (x|x

(t)

, Σ

)

centered around the current step x

(t)

with a covariance matrix Σ

. This proposal corresponds

to the classical random walk Metropolis (RWM) sampler. Note that with a symmetric proposal

distribution, a candidate x

(⋆)

is always accepted if one has π(x

(⋆)

|Y) ≥ π(x

(t)

|Y), i.e., if it

is more likely to belong to the posterior distribution than x

(t)

. For π(x

(⋆)

|Y) < π(x

(t)

|Y),

however, the proposed candidate is not rejected, but accepted only with probability α =

π(x

(⋆)

|Y)/π(x

(t)

|Y).

In practice, in order to accept and reject the proposed candidates with the probability in

Eq. (1.32) or Eq. (1.33), one usually samples a random variate u ∈ [0, 1] according to a

standard uniform distribution U ∼ U(u|0, 1) and compares it to the ratio α from Eq. (1.32).

If α ≥ u, then the proposed candidate is accepted. Otherwise, in the event that α < u, the

proposed candidate is rejected.

UQ[PY]LAB-V1.0-113 - 10 -

Bayesian inference for model calibration and inverse problems

1.3.2 Adaptive Metropolis algorithm

A practical weak point of the standard Metropolis–Hastings algorithm is the need to choose a

proposal distribution p(x

(⋆)

(t)

). Ideally, this distribution should be as similar to the poste-

rior distribution as possible. In most applications, however, the posterior shape is not known

a priori. Moreover, a badly chosen proposal distribution signiﬁcantly affects the MCMC per-

formance up to the point where the MCMC algorithm fails because it does not accept any

proposed candidates. This typically occurs in high dimensions with strongly correlated pos-

terior distributions.

A workaround was proposed in Haario et al. (2001). In this approach known as adaptive

Metropolis algorithm (AM), the Gaussian proposal distribution of the classic Metropolis algo-

rithm (see Section 1.3.1) is tuned during the sampling procedure based on previously gen-

erated samples. The algorithm starts as a standard random walk Metropolis algorithm with

an initial proposal covariance C

. Following a starting period t

, the proposal distribution

covariance matrix is updated to:

C(t + 1) =

(

, t + 1 ≤ t

(M)

C(t), t + 1 > t

(1.34)

where s

(M) =

2.38

is a tuning parameter that depends only on the dimension of the prob-

lem (Gelman et al., 1996). The empirical covariance

C(t) is estimated based on available

sample points generated up to step t and can be computed through:

C(t) =

t − 1

i=1

(i)

−

(t)

)(x

(i)

−

(t)

)

⊺

, where

(t)

i=1

(i)

. (1.35)

In practical applications, an iterative approach can be used to compute the empirical covari-

ance matrices

C(t) with a negligible computational burden (i.e., updated from one step t to

the next) (Haario et al., 2001). To avoid singularity of this estimated covariance matrix, a

small constant ϵ is added to the diagonal of its correlation matrix. This modiﬁed empirical

covariance matrix C

⋆

(t) is then used in the deﬁnition of the Gaussian proposal distribution

centered at the current step of the Markov chain:

p(x|x

(1)

, . . . , x

(t)

) = N (x|x

(t)

, C

⋆

(t)). (1.36)

After drawing a candidate point x

(⋆)

from the proposal distribution p(x|x

(1)

, . . . , x

(t)

), this

candidate is accepted with probability:



(⋆)

, x

(t)



= min

(

π(x

(⋆)

|Y)

π(x

(t)

|Y)

)

, (1.37)

which is the acceptance probability already deﬁned in Eq. (1.32) for symmetric proposal

UQ[PY]LAB-V1.0-113 - 11 -

UQ[PY]LAB user manual

distributions. As the transition probability K(x

(t+1)

(0)

, . . . , x

(t)

) at each step t depends on

all previous steps through the proposal distribution, the generated chain is non-Markovian

and it does not fulﬁl the symmetry condition from Eq. (1.29). Nonetheless, it was shown in

Haario et al. (2001) that the generated sample points can be used to approximate posterior

properties by Eq. (1.31).

1.3.3 Hamiltonian Monte Carlo algorithm

Instead of purely relying on a random walk to sample from the posterior distribution, Hamil-

tonian Monte Carlo algorithms (HMC) exploit the gradient of the posterior distribution to

construct a Markov chain using Hamiltonian dynamics. The connection between Hamilto-

nian dynamics and MCMC algorithms was originally established in Duane et al. (1987) and

a more detailed description of the HMC algorithm can be found in Neal (2011); Nagel and

Sudret (2016a).

The core idea of the algorithm lies in randomly assigning a momentum to a particle and

letting it travel over a potential surface. This can be formalized by deﬁning a potential U (x)

and a kinetic energy function K(p):

U(x) = − log (π(x)L(x; Y)) , K(p) =

⊺

−1

, (1.38)

where p = (p

, . . . , p

)

⊺

is the momentum vector and M is a mass matrix for the particle,

which is often assumed to be a diagonal or simply a scaled identity matrix. This mass property

can be considered a tuning parameter of the algorithm.

The Hamiltonian Monte Carlo algorithm then uses the energy functions from Eq. (1.38) to

deﬁne the Hamiltonian:

H(x, p) = U(x) + K(p). (1.39)

The Hamiltonian captures the total energy of a particle at position x with a given momentum

p. According to Hamiltonian dynamics, the movement of a particle in this system can be

calculated by (where the dot denotes the time derivative):

˙x

∂H(x, p)

∂p

, ˙p

∂H(x, p)

∂x

, for i = 1, . . . , M. (1.40)

These equations can be solved by the well-known leapfrog integration algorithm (Neal, 2011).

It starts out at the current position x(0) = x

(t)

of the particle and a given momentum p(0)

drawn from the distribution:

p(0) ∼ N (p|0, I

). (1.41)

The leapfrog algorithm then solves the Hamiltonian equations for a total duration τ , using a

UQ[PY]LAB-V1.0-113 - 12 -

Bayesian inference for model calibration and inverse problems

discrete timestep size of

. The result is the proposal position x(τ ) and momentum p(τ) of

the particle at time τ.

Given the new location and momentum of the particle, one again computes the Hamiltonian

and accepts the new candidate with probability:

α (x(τ), p(τ), x(0), p(0)) = min {1, exp (H(x(0), p(0)) − H(x(τ ), p(τ)))} . (1.42)

If the new candidate is accepted, the next position of the Markov chain is set to x

(t+1)

= x(τ);

if it is rejected, it is set to x

(t+1)

= x

(t)

. It was shown in Neal (2011) that the invariant

distribution of the generated chain is the posterior distribution.

As the Hamiltonian is generally invariant between the initial point H(x(0), p(0)) and last

point H(x(τ ), p(τ)) of the dynamic simulation, the acceptance probability is theoretically

always one (i.e., no proposal points are rejected). Due to the numerical integration carried

out by the leapfrog method this is only true approximately. Nevertheless, the acceptance

probability of Hamiltonian Monte Carlo is typically close to one.

1.3.4 Afﬁne invariant ensemble algorithm

Most MCMC algorithms perform poorly when the target (i.e., posterior) distribution shows

strong correlation between the parameters. The performance of these algorithms can typ-

ically only be improved by considerable amount of tuning. The afﬁne invariant ensemble

algorithm (AIES) originally presented in Goodman and Weare (2010) alleviates this prob-

lem. It has the desirable property of being invariant to afﬁne transformations of the target

distribution. This means that if there exists an afﬁne transformation of the difﬁcult-to-sample

(by standard MCMC methods) target distribution to an easier-to-sample target distribution,

AIES samples both distributions equally easily without explicitly requiring this afﬁne trans-

formation.

The algorithm simultaneously runs an ensemble of C Markov chains {X

, . . . , X

}, where

each chain is called a walker. The Markov chain locations x

are updated walker by walker.

One such update consists of picking randomly a conjugate walker x

(t)

from the set of walkers

excluding the current i-th walker (j ̸= i).

The afﬁne invariance property is achieved by generating proposals according to a so-called

stretch move. This refers to proposing a new candidate by:

(⋆)

= x

(t)

+ Z ·



(

− x

(t)



, (1.43)

where

t = t + 1 if j < i and

t = t otherwise, i.e.it denotes the latest state of the j-th walker.

UQ[PY]LAB-V1.0-113 - 13 -

UQ[PY]LAB user manual

Z is randomly drawn from the PDF

p(z|a) =







√

z(2

√

a−

√

)

if z ∈ [1/a, a] ,

0 otherwise,

(1.44)

The candidate x

(⋆)

is then accepted as the new location of the i-th walker with probability:



(⋆)

, x

(t)

, z



= min

(

1, z

M−1

π(x

(⋆)

|Y)

π(x

(t)

|Y)

)

. (1.45)

This is repeated for all C walkers in the ensemble. The resulting chains fulﬁll the detailed

balance condition and the generated sample can thus be combined to estimate expectations

under the posterior distribution using Eq. (1.31). A practical advantage of the AIES algorithm

is that it only has a single scalar tuning parameter a, which is often set to a = 2 (Goodman

and Weare, 2010; Allison and Dunkley, 2013; Wicaksono, 2017). On the other hand, due to

its sequential nature, the algorithm cannot be parallelized which makes it comparably slow.

1.3.5 Assessing convergence in MCMC simulations

All MCMC algorithms produce chains of sample points that will eventually follow the pos-

terior distribution. In practice, however, one is forced to make decisions about convergence

based on a ﬁnite number of sample points. As MCMC algorithms lack a convergence crite-

rion, numerous heuristics have been developed to allow practitioners to assess the quality of

the produced Markov chains.

1.3.5.1 Acceptance rate

The acceptance rate r

gives a quantitative indication of how many proposed sample points

were accepted. It can be simply computed as the ratio between the number of accepted

points and the total number of iterations T + 1.

In MCMC algorithms, the acceptance rate depends mostly on their tuning parameters. For

the Metropolis algorithm with a Gaussian proposal, the optimal acceptance rate is shown

to approach r

= 0.23 as M → ∞ (Roberts et al., 1997). For Hamiltonian Monte Carlo

algorithms, it was already mentioned that r

is typically close to one. It is difﬁcult to assess

the quality of a generated MCMC chain purely based on the computed acceptance rate, but

it can serve as an indicator of a badly tuned algorithm.

In practical applications, acceptance rates close to one (except in the Hamiltonian Monte

Carlo algorithm) typically indicate that the proposal distribution does not sufﬁciently explore

the target distribution. Acceptance rates close to zero indicate instead that the proposed

UQ[PY]LAB-V1.0-113 - 14 -

Bayesian inference for model calibration and inverse problems

(a) After 500 steps

(b) After 10

steps

Figure 1: Trace plot and corresponding KDE at two different iterations of the chain.

candidate points are in low probability regions. The most common reasons for this are too

wide proposal distributions or proposal distributions that do not sufﬁciently resemble the

target distributions.

1.3.5.2 Trace and density plots

Trace plots show the evolution of a Markov chain. As chains are typically initialized at ran-

dom points, the evolution of an MCMC chain can give valuable insights about convergence.

Trace plots are typically assessed visually for each dimension individually.

A sample generated by the chain should eventually be distributed according to the posterior

distribution. A kernel density estimation (KDE) scheme (Wand and Jones, 1995) can thus be

employed to obtain an approximation of the posterior marginal. If the chain has reached its

steady state, this KDE of the posterior marginal should not change considerably with further

iterations.

An example of a trace plot with a corresponding KDE is displayed in Figure 1. It can be clearly

seen that the chain has not reached its steady state after 500 steps (Figure 1a) whereas it has

after 10

steps (Figure 1b).

UQ[PY]LAB-V1.0-113 - 15 -

UQ[PY]LAB user manual

1.3.5.3 Gelman-Rubin diagnostics

A quantitative approach to assess convergence was introduced by Gelman and Rubin (1992)

and later generalized by Brooks and Gelman (1998). The idea presented there is to compare

a set of C independent Markov chains that were initiated at different seed points. If the

chains are converged, the empirical second moments computed from the individual chains

should be the same as the empirical second moments computed from combining the samples

from all C chains.

Figure 2: Convergence of the MPSRF

More formally, let {X

, . . . , X

} be the C chains run in parallel from different seed points

(0)

, where each chain X

= (x

(0)

, . . . , x

(T )

) contains T + 1 sample points with x

(t)

∈ R

The Gelman-Rubin diagnostic requires the computation of two covariance matrices. The

covariance matrix of the i-th chain is estimated by:

t=0



(t)

−



(t)

−



⊺

T + 1

t=0

(t)

. (1.46)

The matrices for all C chains are then averaged to obtain the within-sequence covariance

W ∈ R

M×M

W =

i=1

. (1.47)

The second matrix required is the so-called between-sequence variance B ∈ R

M×M

. It cap-

tures the covariance between the individual MCMC chains and is estimated as:

B =

C − 1

i=1

(

−

x) (

−

⊺

, where: (1.48)

x =

C(T + 1)

i=1

t=0

(t)

(1.49)

is the average of all (T + 1) states of the C chains. To estimate the difference between the

UQ[PY]LAB-V1.0-113 - 16 -

Bayesian inference for model calibration and inverse problems

within-sequence covariance estimate W and the between-sequence covariance estimate B, the

following multivariate potential scale reduction factor (MPSRF) is proposed in Brooks and

Gelman (1998):

T + 1



C + 1



, (1.50)

where λ

is the largest eigenvalue of the symmetric positive deﬁnite matrix W

−1

B. This

approaches 1 (from above) with increasing convergence of the MCMC algorithm. This

convergence is showcased for a sample MCMC chain in Figure 2. This method requires a set

of independent, parallel MCMC chains.

1.3.5.4 Burn-in

Once an MCMC chain has reached its steady state, the generated sample follows the posterior

distribution. When, however, a ﬁnite number of sample points is used to estimate posterior

properties (e.g. moments), the sample points generated prior to convergence can pollute the

estimation.

It is therefore common practice in practical MCMC applications to discard sample points that

were generated prior to convergence (Brooks et al., 2011). This discarded fraction is called

burn-in.

1.4 Spectral likelihood expansion

A different way to solving Bayesian inverse problem named spectral likelihood expansion

(SLE) was proposed in Nagel and Sudret (2016b). This sampling-free approach uses spectral

expansion techniques such as polynomial chaos expansion (PCE, UQ[PY]LAB User Manual –

Polynomial Chaos Expansions) to approximate the likelihood function at the core of Bayesian

inverse problems.

Likelihood functions can be seen as scalar functions of the input random vector X ∼ π(x)

with ﬁnite output variance (Nagel and Sudret, 2016b). Assuming independent priors, i.e.

π(x) =

i=1

), their spectral expansion reads:

L(X) ≈ L

SLE

(X)

def

α∈A

(X), (1.51)

where Ψ

are basis functions (polynomials in the case of PCE) that are orthogonal w.r.t. the

prior distribution π(x) and c

are the corresponding coefﬁcients. Following this approxi-

mation, it becomes possible to exploit the orthogonality of the spectral basis functions to

post-process the expansion coefﬁcients and extract the following important posterior quanti-

ties of interest:

UQ[PY]LAB-V1.0-113 - 17 -

UQ[PY]LAB user manual

Evidence The evidence emerges as the coefﬁcient of the constant polynomial a

Z =

L(x)π(x) dx ≈ E [L

SLE

(X)] = a

. (1.52)

Posterior Upon computing the evidence Z, the posterior can be evaluated directly through

π(x|Y) ≈

SLE

(x)π(x)

π(x)

α∈A

(x). (1.53)

Posterior marginals An approximation of the univariate posterior marginals deﬁned in Eq. (1.24)

can then also be derived analytically:

|Y) =

∼i

π(x|Y) dx

∼i

≈

)

α∈A

∼i=0

), (1.54)

where A

∼i=0

= {α ∈ A: α

= 0 ⇔ j ̸= i}.

Similarly, an expression for the multivariate posterior marginal deﬁned in Eq. (1.25)

can be derived. Denote by π

)

def

i∈u

) and π

)

def

i∈v

) the prior

marginal density functions of X

and X

respectively, the posterior marginal then

reads:

|Y) =

π(x|Y) dx

≈

)

α∈A

v=0

), (1.55)

where A

v=0

= {α ∈ A: α

= 0 ⇔ i ∈ v}. The series in Eq. (1.54) is a subexpansion

that contains non-constant polynomials only along the dimensions i ∈ u.

Quantities of interest It is also possible to analytically compute posterior expectations of

functions that admit a polynomial chaos expansion on the same basis, of the form

h(X) ≈

α∈A

(X). Eq. (1.26) then reduces to the spectral product:

E [h(X)|Y] =

α∈A

. (1.56)

This expression can be used to compute posterior moments like mean, variance or

covariance.

The quality of these results depends only on the approximation error introduced in Eq. (1.51).

The latter, in turn, depends mainly on the chosen PCE truncation strategy (Nagel and Su-

dret, 2016b; L

uthen et al., 2020) and the number of points used to compute the coefﬁcients

(i.e. the experimental design). It is known that informative likelihood functions have quasi-

compact supports (i.e. L(X) ≈ 0 on a majority of D

). Such functions require a very high

polynomial degree to be approximated accurately, which in turn can lead to the need for

prohibitively large experimental designs.

UQ[PY]LAB-V1.0-113 - 18 -

Bayesian inference for model calibration and inverse problems

1.5 Stochastic spectral likelihood embedding

An extension of SLE (see Section 1.4), namely stochastic spectral likelihood embedding (SSLE),

was recently proposed in Wagner et al. (2021). This approach, similarly to SLE, directly ap-

proximates the likelihood function, but replaces the global spectral expansion in Eq. (1.51)

for an SSE metamodel ((Marelli et al., 2021), see also UQ[PY]LAB User Manual – Stochas-

tic spectral embedding). Due to their local approximation strength, SSE metamodels are

more suitable for approximating functions with quasi-compact supports than purely global

approximation approaches.

Again, viewing the likelihood as a function of a random vector X with independent compo-

nents, we can write its SSLE representation as:

L(X) ≈ L

SSLE

(X)

def

k∈K

(X)

(X), (1.57)

where

(X) =

α∈A

(X) (1.58)

are residual expansions.

The variable X is distributed according to the prior distribution π(x) and, consequently, the

local basis Ψ

is orthonormal w.r.t. that distribution.

Due to the local spectral properties of the residual expansions, the SSLE representation of

the likelihood function retains all of the post-processing properties of SLE (see Section 1.4):

Evidence The normalization constant Z emerges as the sum of the constant polynomial

coefﬁcients weighted by the prior mass:

Z =

k∈K

α∈A

(x)π(x) dx =

k∈K

, where V

π(x) dx.

(1.59)

Posterior This allows us to write the posterior density as

π(x|Y) ≈

SSLE

(x)π(x)

π(x)

k∈K

(x)

(x). (1.60)

Posterior marginals Utilizing the disjoint sets u and v from Eq. (1.54) it is also possible to

analytically derive posterior marginal PDFs as

|Y) =

π(x|Y) dx

≈

)

k∈K

)

S,u

(1.61)

UQ[PY]LAB-V1.0-113 - 19 -

UQ[PY]LAB user manual

where

S,u

) =

α∈A

v=0

) and V

) dx

. (1.62)

S,u

) is a subexpansion of

(x) that contains only non-constant polynomials along

the dimensions i ∈ u. Note that, as we assumed that the prior distribution has inde-

pendent components, the constants V

and V

are obtained as products of univariate

integrals which are available analytically from the prior marginal cumulative distribu-

tion functions (CDFs).

Quantities of interest Posterior expectations of a function h(x) =

α∈A

(x) for k ∈

K can be approximated by:

E [h(X)|Y] =

h(x)π(x|Y) dx

k∈K

α∈A

h(x)Ψ

(x)π(x) dx

k∈K

α∈A

(1.63)

where b

are the coefﬁcients of the PCE of h in the card(K) bases {Ψ

}

α∈A

. The same

expression can also be used for computing posterior moments like mean, variance or

covariance.

These expressions can be seen as a generalization of the ones for SLE detailed in Section 1.4.

For a single-level global expansion (i.e. card(K) = 1) and consequently V

(0,1)

= 1, they are

identical.

1.5.1 Sequential partitioning algorithm

The algorithm to construct SSEs implemented in UQ[PY]LAB is called sequential partitioning

algorithm. It sequentially partitions selected reﬁnement domains and constructs local expan-

sions of the residual to ultimately produce a likelihood approximation of the form shown in

Eq. (1.57). This algorithm is explained in detail in the UQ[PY]LAB User Manual – Stochastic

spectral embedding. In Wagner et al. (2021), modiﬁcations to this algorithm were proposed

that were shown to improve the SSE approximation accuracy for likelihood functions in SSLE.

These modiﬁcations pertain to the partitioning and sample enrichment strategies.

1.5.1.1 Partitioning strategy

The partitioning strategy determines how a selected reﬁnement domain is split. For likelihood

functions, it was proposed in Wagner et al. (2021) to pick the split direction along which a

UQ[PY]LAB-V1.0-113 - 20 -

Bayesian inference for model calibration and inverse problems

(a) Reﬁnement domain (b) Split along d = 1 (c) Split along d = 2 (d) Selected pair

Figure 3: Partitioning strategy for a 2D example visualized in the quantile space U . The

reﬁnement domain D

ℓ,p

is split into two subdomains D

ℓ+1,s

and D

ℓ+1,s

split yields a maximum difference in the residual empirical variance between the two candi-

date subdomains created by the split. This can easily be visualized with an example given by

the M = 2 dimensional domain D

ℓ,p

in Figure 3a. Assume this subdomain was selected as

the reﬁnement domain. To decide along which dimension to split, we construct the M candi-

date subdomain pairs {D

i,1

split

, D

i,2

split

}

i=1,...,M

and estimate the corresponding {E

split

}

i=1,...,M

in those subdomains deﬁned by

split

def



Var

ℓ+1

i,1

split

)

− Var

ℓ+1

i,2

split

)



. (1.64)

In this expression, X

i,1

split

and X

i,2

split

denote subsets of the experimental design X that lie within

the subdomains D

i,1

split

and D

i,2

split

respectively. The corresponding variances can be estimated

with the empirical variance of the residuals in the respective candidate subdomains.

After computing the residual variance differences, the split is carried out along dimension

d = arg max

i∈{1,...,M}

split

, (1.65)

i.e. to keep the subdomains D

d,1

split

and D

d,2

split

that introduce the largest difference in variance.

For d = 1, the resulting split can be seen in Figure 3d.

1.5.1.2 Sample enrichment

Likelihood functions are typically not equally informative at every point of the input space.

On the contrary, it is often the case that large input regions are close to zero, while only a

small subdomain of the input space yields likelihood responses that are multiple orders of

magnitude larger.

It was shown in Wagner et al. (2021) that it is therefore more efﬁcient in terms of total

number of likelihood evaluations, to sequentially sample the experimental design in SSLE.

The rationale is to enrich the experimental design only in domains that have been selected

for reﬁnement. Informative regions of the likelihood function (that are more difﬁcult to ap-

UQ[PY]LAB-V1.0-113 - 21 -

UQ[PY]LAB user manual

proximate) will then receive a greater share of likelihood evaluations, because the reﬁnement

domain selection in the sequential partitioning algorithm is based on the local approximation

error.

UQ[PY]LAB-V1.0-113 - 22 -

Chapter 2

Usage

In this chapter the implementation of Bayesian inversion as discussed in Chapter 1 is de-

scribed. A simple engineering inverse problem is solved with UQ[PY]LAB to exemplify the

usage of the Bayesian inversion module. This simple example is then extended to treat more

complex problems.

2.1 Reference problem: calibration of a simply supported beam

model

We consider a simply-supported beam such as the one shown in Figure 4. The beam has a

known rectangular cross-section of width b and height h and a known span of length L.

Figure 4: Simple beam bending test.

A set of N = 5 independent experiments are carried out with this beam, with the goal of

inferring the unknown material stiffness, i.e. its Young’s modulus E. In the experiments,

the beam is subject to a constant distributed load p and the mid-span deﬂection V

mid

is mea-

sured. The measurements are reported in Table 1. Due to measurement error, the measured

deﬂections vary across experiments.

Table 1: Beam deﬂection: measured deﬂections.

Experiment 1 2 3 4 5

mid

(mm) 12.84 13.12 12.13 12.19 12.67

UQ[PY]LAB user manual

The parameters (b, h, L, p) are considered known and their values are given in Table 2.

Table 2: Beam experiments: nominal values of the beam properties.

Variable Nominal Value

b (m) 0.15

h (m) 0.3

L (m) 5

p (N/m) 12 000

The analytical expression for the mid-span deﬂection V

mid

of the beam according to the stan-

dard beam theory is:

mid

Ebh

. (2.1)

This simple equation serves as the forward model and relates the unknown Young’s modulus

to the measurable mid-span deﬂection.

Additionally, it is known from prior experiments that the Young’s modulus of the material

follows a lognormal distribution:

E ∼ LN (λ, ζ), with µ

= 30 (GPa) and σ

= 4.5 (GPa). (2.2)

In Bayesian inversion terms, the prior information on the model parameter x

≡ E is E ∼

LN (λ, ζ). Due to a lack of more information, an unknown additive Gaussian experimental

discrepancy model is assumed. As a weakly informative prior on the positive discrepancy

variance x

≡ σ

, a uniform distribution σ

∼ U(0, µ

mid

) with µ

mid

equal to the empirical

mean of the observations given in Table 1.

2.2 Problem setup and solution

Solving an inverse problem with the Bayesian inverse module of UQ[PY]LAB typically re-

quires the speciﬁcation of a prior distribution π(x), N independent observations Y = {y

, . . . , y

where y

∈ R

out

and a forward model M. All these ingredients are brieﬂy discussed in this

section.

2.2.1 Initialize UQ[PY]LAB

The ﬁrst step is to initialize UQ[PY]LAB and ﬁxing a random seed for reproducibility:

# Package imports

from uqpylab import sessions

import numpy as np

# Start the session

mySession = sessions.cloud()

UQ[PY]LAB-V1.0-113 - 24 -

Bayesian inference for model calibration and inverse problems

# (Optional) Get a convenient handle to the command line interface

uq = mySession.cli

# Reset the session

mySession.reset()

# Set the random seed for reproducibility

uq.rng(100, 'twister');

2.2.2 Specify a prior distribution

The prior distribution of the model parameters is deﬁned as an INPUT object by:

PriorOpts = {

"Marginals": [

{

"Name": "b", # beam width

"Type": "Constant",

"Parameters": [0.15] # (m)

{

"Name": "h", # beam height

"Type": "Constant",

"Parameters": [0.3] # (m)

{

"Name": "L", # beam length

"Type": "Constant",

"Parameters": [5] # (m)

{

"Name": "E", # Young's modulus

"Type": "LogNormal",

"Moments": [30e9,4.5e9] # (N/mˆ2)

{

"Name": "p", # constant distributed load

"Type": "Constant",

"Parameters": [12000] # (N/m)

}

]

}

myPriorDist = uq.createInput(PriorOpts)

As the prior distribution is speciﬁed as an INPUT object, all the features of the UQ[PY]LAB

INPUT module (UQ[PY]LAB User Manual – the INPUT module) can be used, i.e., constant and

marginal distributions of any kind (including user-deﬁned ones). Dependence may also be

speciﬁed using copulas.

Note: The known parameters from Table 2 are deﬁned as Constant input marginals

and will not be considered during the calibration procedure, except when evalu-

ating the forward model.

UQ[PY]LAB-V1.0-113 - 25 -

UQ[PY]LAB user manual

2.2.3 Create a forward model

The computational model, given in Eq. (2.1), is deﬁned as a Python file called

SimplySupportedBeam.py which you can download on https://uqpylab.uq-cloud.

io/examples. A UQ[PY]LAB MODEL is then created as:

ModelOpts = {

'Type': 'Model',

'Name': 'Forward model',

'ModelFun': 'SimplySupportedBeam.model'

}

myForwardModel = uq.createModel(ModelOpts)

For more details about the conﬁguration options available for a MODEL object, please refer to

the UQ[PY]LAB User Manual – the MODEL module.

2.2.4 Provide measurements

The measurements Y are stored in an N × N

out

list of lists V_mid, where N

out

is the number

of model outputs:

V_mid = np.array([12.84, 13.12, 12.13, 12.19, 12.67])/1000 # (m)

myData = {

'y': V_mid.tolist(),

'Name': 'Beam mid-span deflection'

}

Note: In the present case where N

out

= 1, this list of lists reduces to a list.

2.2.5 Perform the Bayesian inverse analysis

The options are then gathered in a Python dictionary, here called BayesOpts:

BayesOpts = {

'Type': 'Inversion',

'Data': myData

}

The BayesOpts dictionary contains all information required to solve the inverse problem.

If not explicitly speciﬁed by the user, by default the Bayesian inversion module uses the last

created INPUT object (in this case myPriorDist) as a prior distribution and the last created

MODEL object (in this case myForwardModel) as the forward model. Therefore, to perform

the analysis, it is sufﬁcient to create the corresponding ANALYSIS object:

myBayesianAnalysis = uq.createAnalysis(BayesOpts)

UQ[PY]LAB-V1.0-113 - 26 -

Bayesian inference for model calibration and inverse problems

Note: Without an explicitly-speciﬁed discrepancy model, UQ[PY]LAB assumes by de-

fault an unknown additive Gaussian discrepancy term, with a single unknown

residual parameter x

≡ σ

. The prior distribution of σ

is a weakly informative

uniform distribution σ

∼ U(0, µ

) with µ

equal to the empirical mean of the

provided data Y.

Note: If not otherwise speciﬁed UQ[PY]LAB uses the afﬁne invariant ensemble sampler

(Section 1.3.4) with C = 100 parallel chains and initial points drawn from the

prior distribution and performs T = 300 iterations.

The sample set generated by MCMC algorithms often needs to be post-processed before it

can be used as a true posterior sample (e.g., to remove burn-in, Section 1.3.5.4). By default

in UQ[PY]LAB the ﬁrst half of the sample points generated by all chains are removed as

burn-in (see Section 1.3.5.4) and the empirical parameter mean E [X |Y] along with the

2.5th to 97.5th percentile based on the post-processed sample are estimated. Additionally,

samples from the prior distribution and from the posterior predictive distribution are drawn

(see Section 1.2.6) and stored in the myBayesianAnalysis object. For all available post-

processing options see Section 2.4.1.6 and Section 3.3.

A brief report of the analysis can then be generated by:

uq.print(myBayesianAnalysis)

which produces:

%----------------------- Inversion output -----------------------%

Number of calibrated model parameters: 1

Number of non-calibrated model parameters: 4

Number of calibrated discrepancy parameters: 1

%------------------- Data and Discrepancy

% Data-/Discrepancy group 1:

Number of independent observations: 5

Discrepancy:

Type: Gaussian

Discrepancy family: Scalar

Discrepancy parameters known: No

Associated outputs:

Model 1:

Output dimensions: 1

%------------------- Solver

Solution method: MCMC

Algorithm: AIES

Duration (HH:MM:SS): 00:00:55

Number of sample points: 3.00e+04

%------------------- Posterior Marginals

---------------------------------------------------------------------

UQ[PY]LAB-V1.0-113 - 27 -

UQ[PY]LAB user manual

---------------------------------------------------------------------

| E | 2.4e+10 | 1.7e+09 | (2.2e+10 - 2.9e+10) | Model |

| Sigma2 | 3.7e-06 | 1.1e-05 | (1.1e-07 - 2.8e-05) | Discrepancy |

---------------------------------------------------------------------

%------------------- Point estimate

----------------------------------------

| Parameter | Mean | Parameter Type |

----------------------------------------

| E | 2.4e+10 | Model |

| Sigma2 | 3.7e-06 | Discrepancy |

----------------------------------------

The results can also be visualized by:

uq.display(myBayesianAnalysis)

which produces the images in Figure 5.

2.2.6 Advanced options: discrepancy model

The discrepancy model deﬁnes the connection between the supplied data and the forward

model. The Bayesian module of UQ[PY]LAB currently supports the option to specify a

Gaussian additive discrepancy as deﬁned in Eq. (1.12).

In real applications, it is often beneﬁcial to specify more accurately the discrepancy model.

The options are to either specify a known residual variance σ

(e.g. tabulated measurement

discrepancies from the instrument supplier) or an unknown x

= σ

which can be inferred

together with the model parameters x

2.2.6.1 Known residual variance

If the variance is known a priori, as detailed in Section 1.2.2, it can be directly deﬁned in an

DiscrepancyOpts dictionary (assuming that σ

= 10

−7

)):

DiscrepancyOpts = {

'Type': 'Gaussian',

'Parameters': 1e-7 # (mˆ2)

}

Which is then passed to BayesOpts by:

BayesOpts['Discrepancy'] = DiscrepancyOpts

UQ[PY]LAB-V1.0-113 - 28 -

Bayesian inference for model calibration and inverse problems

(a) Scatterplots of prior and posterior sample with the posterior mean point estimator from (1.31).

(b) Posterior predictive distribution from (1.28), data, and model at mean pre-

diction obtained by propagating the mean point estimator through the forward

model.

Figure 5: Visualize analysis: The results of the Bayesian inverse analysis on the input and the

model predictions.

UQ[PY]LAB-V1.0-113 - 29 -

UQ[PY]LAB user manual

2.2.6.2 Unknown residual variance

If σ

is not known a priori, as detailed in Section 1.2.3, the Bayesian framework can infer the

distribution of the discrepancy parameter x

≡ σ

. This requires the initial speciﬁcation of a

prior distribution of the discrepancy parameter π(σ

) (see Eq. (1.2)).

The prior distribution of the parameter π(σ

) can be deﬁned as a UQ[PY]LAB INPUT object

and then assigned to the DiscrepancyOpts dictionary’s Prior key:

DiscrepancyPriorOpts = {

'Name': 'Prior of discrepancy parameter',

'Marginals': {

'Name': 'Sigma2',

'Type': 'Uniform',

'Parameters': [0, np.mean(V_mid)

}

myDiscrepancyPrior = uq.createInput(DiscrepancyPriorOpts)

DiscrepancyOpts = {

'Type': 'Gaussian',

'Prior': myDiscrepancyPrior['Name']

}

Which is then passed to BayesOpts by:

BayesOpts['Discrepancy'] = DiscrepancyOpts

If multiple UQ[PY]LAB INPUT objects are deﬁned, the prior distribution must be speciﬁed by:

BayesOpts['Prior'] = myPriorDist['Name']

Note: Here a uniform prior on σ

with bounds U(0, µ

mid

) is implemented. This is the

default (see above Section 2.2.5) when no discrepancy options are provided.

Note: As the prior π(σ

) is deﬁned for the variance σ

, only distributions with positive

support can be used here.

2.3 Multiple model outputs

Models with multiple outputs are often encountered in real calibration problems. These mul-

tiple outputs can be different measurable quantities (e.g., temperature and displacement),

quantities at different locations (e.g., deformations at different physical points), or at differ-

ent times (e.g., time series).

To show how such problems can be treated in UQ[PY]LAB, the reference problem from

UQ[PY]LAB-V1.0-113 - 30 -

Bayesian inference for model calibration and inverse problems

Section 2.1 is slightly extended. It is assumed that additionally to measurements of the

deﬂection at the beam mid-span at L/2, measurements are also available at L/4 as shown in

Figure 6.

Figure 6: Simple beam bending test.

The N = 5 measurements of the deﬂections V

mid

and V

L/4

are given in Table 3.

Table 3: Beam deﬂection: measured deﬂections.

Experiment 1 2 3 4 5

L/4

(mm) 8.98 8.66 8.85 9.19 8.64

mid

(mm) 12.84 13.12 12.13 12.19 12.67

Similarly to (2.1), the quarter deﬂection at L/4 can be computed analytically by the standard

beam theory:

L/4

512

Ebh

. (2.3)

Again due to a lack of information on the discrepancy, two independent unknown experimen-

tal discrepancies ε

and ε

are assumed for the measured displacements V

mid

and V

L/4

, re-

spectively. As a weakly informative prior on the positive discrepancy variances x

= (σ

, σ

two independent uniform distributions π(σ

) = U(0, µ

) are chosen with µ

equal to the

mean of the observations V

mid

and V

L/4

respectively (see Table 3).

2.3.1 Create a forward model

The equations for the beam mid-span (L/2) and quarter-span (L/4) deﬂection, given in

Eqs. (2.1) and (2.3), are implemented as a Python file called

SimplySupportedBeamTwo.py which you can download on https://uqpylab.uq-cloud.

io/examples. This function returns the two beam deﬂections for a single model parameter

realization x

in a row list of length N

out

= 2.

This is added to UQ[PY]LAB with the following commands:

ModelOpts = {

'Type': 'Model',

'Name': 'Forward model',

'ModelFun': 'SimplySupportedBeamTwo.model'

}

myForwardModel = uq.createModel(ModelOpts)

UQ[PY]LAB-V1.0-113 - 31 -

UQ[PY]LAB user manual

2.3.2 Provide measurements

The measurements are stored in an N × N

out

list of lists V and assigned to the BayesOpts

object:

V = np.array([[8.98, 8.66, 8.85, 9.19, 8.64], # L/4 (m)

[12.84, 13.12, 12.13, 12.19, 12.67]])/1000 # L/2 (m)

V = V.T # NxN_out array, where N_out is the number of model outputs

myData = {

'y': V.tolist(),

'Name': 'Beam quarter and midspan deflection',

}

Note: By default, UQ[PY]LAB assigns the i-th list of the list to the i-th output of the

forward model. For more advanced options, refer to Section 2.3.4.3.

Note: Without an explicitly-speciﬁed discrepancy model, the Bayesian inversion module

by default assumes unknown additive Gaussian discrepancies for all N

out

outputs

of the forward model, with the residual parameter vector x

= (σ

, . . . , σ

out

The prior distribution is by default taken as independent uniform distributions

π(x

) =

out

i=1

π(σ

) where π(σ

) = U(0, µ

) with µ

equal to the empirical

mean of measurements available for the i-th output dimension.

2.3.3 Perform the inverse analysis

The analysis can be run with:

myBayesianAnalysis = uq.createAnalysis(BayesOpts)

To generate a prior predictive sample, call:

myBayesianAnalysis = uq.postProcessInversion(myBayesianAnalysis,

'priorPredictive',1000)

see Section 2.4.1.6 and Section 3.3 for all available post-processing options.

After the myBayesianAnalysis object is created, the results can be summarized with:

uq.print(myBayesianAnalysis)

which produces the report:

%----------------------- Inversion output -----------------------%

Number of calibrated model parameters: 1

Number of non-calibrated model parameters: 4

Number of calibrated discrepancy parameters: 2

%------------------- Data and Discrepancy

% Data-/Discrepancy group 1:

UQ[PY]LAB-V1.0-113 - 32 -

Bayesian inference for model calibration and inverse problems

Number of independent observations: 5

Discrepancy:

Type: Gaussian

Discrepancy family: Row

Discrepancy parameters known: No

Associated outputs:

Model 1:

Output dimensions: 1

%------------------- Solver

Solution method: MCMC

Algorithm: AIES

Duration (HH:MM:SS): 00:00:46

Number of sample points: 3.00e+04

%------------------- Posterior Marginals

---------------------------------------------------------------------

---------------------------------------------------------------------

| E | 2.3e+10 | 6.6e+08 | (2.2e+10 - 2.5e+10) | Model |

| Sigma2 | 2e-06 | 6.3e-06 | (3.3e-08 - 1.8e-05) | Discrepancy |

| Sigma2 | 4.5e-06 | 1.2e-05 | (1.2e-07 - 3.6e-05) | Discrepancy |

---------------------------------------------------------------------

%------------------- Point estimate

----------------------------------------

| Parameter | Mean | Parameter Type |

----------------------------------------

| E | 2.3e+10 | Model |

| Sigma2 | 2e-06 | Discrepancy |

| Sigma2 | 4.5e-06 | Discrepancy |

----------------------------------------

%------------------- Correlation matrix (discrepancy parameters)

-------------------------------

| | Sigma2 Sigma2 |

-------------------------------

| Sigma2 | 1 -0.045 |

| Sigma2 | -0.045 1 |

-------------------------------

and visualized graphically with:

uq.display(myBayesianAnalysis)

which produces a set of plots similar to those in Figure 7.

2.3.4 Advanced options: discrepancy model

In the case of models with multiple outputs, the full covariance matrix Σ of the residual

discrepancy vector ε = (ε

, . . . , ε

out

) has to be deﬁned for the likelihood function. In this

UQ[PY]LAB-V1.0-113 - 33 -

UQ[PY]LAB user manual

(a) Scatterplots of prior and posterior sample with the posterior mean point estimator from (1.31).

(b) Prior and posterior predictive distributions from (1.27) and (1.28), data,

and model at mean prediction obtained by propagating the mean point estima-

tor through the forward model.

Figure 7: Advanced discrepancy options: The results of the Bayesian inverse analysis on the

input and the model predictions.

UQ[PY]LAB-V1.0-113 - 34 -

Bayesian inference for model calibration and inverse problems

section, all options that are available in the UQ[PY]LAB Bayesian inversion module to deﬁne

this covariance matrix are presented:

1. Known residual variance: Section 2.3.4.1, see also Section 1.2.2

2. Unknown residual variance: Section 2.3.4.2, see also Section 1.2.3

3. Data and discrepancy groups: Section 2.3.4.3

2.3.4.1 Known residual variance

In special cases, the residual variance may be known. This can happen when the forward

model is supposed to perfectly represent the experimental setup, and when the discrepancy

term reduces to measurement error. If the variance of this error is provided by the instrument

supplier, it can be directly used to specify the covariance matrix Σ of the residual vector

ε = (ε

, . . . , ε

out

). This known covariance matrix can be speciﬁed in three different ways in

the Bayesian module of UQ[PY]LAB.

Independent and identically distributed ε

: In the case when all elements ε

of the resid-

ual vector ε independently follow the same distribution N (0, σ

) with a known variance σ

(e.g., σ

= 10

−7

)), this can speciﬁed by a DiscrepancyOpts dictionary like

DiscrepancyOpts = {

'Type': 'Gaussian',

'Parameters': 1e-7 # single scalar

}

This dictionary is then passed to BayesOpts as follows:

BayesOpts['Discrepancy'] = DiscrepancyOpts

Independent ε

: If each element of the residual vector ε

follows a Normal distribution

N (0, σ

) with a speciﬁc known residual variance σ

, but independence can still be assumed,

the DiscrepancyOpts dictionary can be deﬁned as:

DiscrepancyOpts = {

'Type': 'Gaussian',

'Parameters': [[1e-7, 5e-7]] # row vector of length N_out

}

where the length of the Parameters list is equal to N

out

Dependent ε

: In the general case where a known Gaussian distribution can be assumed

for the discrepancy term, the covariance matrix can be passed to UQ[PY]LAB as follows:

UQ[PY]LAB-V1.0-113 - 35 -

UQ[PY]LAB user manual

DiscrepancyOpts = {

'Type': 'Gaussian',

'Parameters': [[1e-7, -5e-8],[-5e-8, 5e-7]] # N_out x N_out matrix

# passed as a nested list

}

This covariance matrix introduces negative correlation between the ﬁrst and second discrep-

ancy parameter σ

and σ

. Any positive-deﬁnite matrix can be used as a covariance matrix.

2.3.4.2 Unknown residual variance

In most practical applications, the parameters σ

are not known a priori. As detailed in Sec-

tion 1.2.3, the Bayesian framework can infer the distribution of the discrepancy parameters

gathered in x

. This requires the initial speciﬁcation of a prior distribution of the discrep-

ancy parameters π(x

) (see Eq. (1.2)). Similarly to the previous section, there are different

ways of specifying an unknown variance parameter σ

, depending on the distribution of the

residuals ε

Note: In contrast to the known residual variance case (see Section 2.3.4.1), the def-

inition of dependent unknown discrepancy terms ε

(see Section 2.3.4.1) is not

currently supported.

If an unknown residual variance is used, it becomes necessary to explicitly assign the INPUT

object deﬁning the prior of the model parameters π(x

) to the BayesOpts dictionary. This

is necessary to avoid confusions between the model parameter INPUT and error parameter

INPUT objects. It can be assigned by:

BayesOpts['Prior'] = myPriorDist['Name']

Independent and identically distributed ε

: If the residuals of all observations are inde-

pendent and identically distributed, a single unknown variance parameter σ

can be used in

the distribution of all residuals ε

. The prior distribution of the parameter π(σ

) can be de-

ﬁned as a UQ[PY]LAB INPUT object and then assigned to the DiscrepancyOpts dictionary:

DiscrepancyPriorOpts = {

'Name': 'Prior of discrepancy parameter',

'Marginals': {

'Name': 'Sigma2',

'Type': 'Uniform',

'Parameters': [0, 1e-4] # (mˆ2)

}

myDiscrepancyPrior = uq.createInput(DiscrepancyPriorOpts)

DiscrepancyOpts = {

'Type': 'Gaussian',

UQ[PY]LAB-V1.0-113 - 36 -

Bayesian inference for model calibration and inverse problems

'Prior': myDiscrepancyPrior['Name']

}

Here a uniform prior π(σ

) = U(0, 10

−4

) (m

) is chosen.

Note: As the prior π(σ

) is deﬁned for the variance σ

, only distributions with positive

support can be used here.

Independent ε

: If the residuals are independent but not identically distributed, the user

can specify a dedicated discrepancy distribution for each σ

. In the present case of N

out

= 2,

two independent prior distributions π(σ

) for the discrepancy parameters have to be speci-

ﬁed. For the sake of illustration a lognormal prior is chosen for σ

and a uniform prior for σ

in the following code:

DiscrepancyPriorOpts = {

'Name': 'Prior of discrepancy parameter',

'Marginals': [

{

'Name': 'Sigma2_1',

'Type': 'Lognormal',

'Moments': [1e-5, 5e-6] # (mˆ2)

{

'Name': 'Sigma2_2',

'Type': 'Uniform',

'Parameters': [0, 1e-4] # (mˆ2)

}

]

}

myDiscrepancyPrior = uq.createInput(DiscrepancyPriorOpts)

DiscrepancyOpts = {

'Type': 'Gaussian',

'Prior': myDiscrepancyPrior['Name']

}

2.3.4.3 Data and discrepancy groups

It often occurs in inverse problems that different types or number of data are collected for in-

dividual model outputs. In this case, it is often also necessary to deﬁne dedicated discrepancy

options for individual outputs y

(see Section 1.2.4).

In UQ[PY]LAB this is achieved through so-called data- and discrepancy groups. The groups

are deﬁned by specifying the DiscrepancyOpts and Data as list of dictionaries, rather than

simple dictionaries. All options that were discussed in the previous sections can then be

assigned to each dictionary in this array independently.

Consider the following case: for the ﬁrst residual ε

the variance σ

is known to be σ

UQ[PY]LAB-V1.0-113 - 37 -

UQ[PY]LAB user manual

−7

), while for the second residual ε

the variance σ

is unknown and assigned a

uniform distribution σ

∼ U(0, 10

−4

) (m

). These discrepancy options can be passed to

UQ[PY]LAB by deﬁning N

= 2 dedicated DiscrepancyOpts and Data dictionaries in the

following way:

# group 1

V_quart = np.array([10.51, 9.60, 10.22, 8.16, 7.47])/1000 # L/4 (m)

Data = [

{

'y': V_quart.tolist(),

'Name': 'Deflection measurements at L/4',

'MOMap': [1] # Model Output Map

}

]

DiscrepancyOpts = [

{

'Type': 'Gaussian',

'Parameters': 1e-7 # (mˆ2)

}

]

# group 2

V_mid = np.array([12.59, 11.23, 15.28, 12.45, 13.21])/1000; # L/2 (m)

Data.append(

{

'y': V_mid.tolist(),

'Name': 'Deflection measurements at L/2',

'MOMap': [2] # Model Output Map

}

)

DiscrepancyPriorOpts = {

'Name': 'Prior of sigma',

'Marginals': {

'Name': 'Sigma2_2',

'Type': 'Uniform',

'Parameters': [0, 1e-4], # (mˆ2)

}

DiscrepancyPrior = uq.createInput(DiscrepancyPriorOpts)

DiscrepancyOpts.append(

{

'Type': 'Gaussian',

'Prior': DiscrepancyPrior['Name']

}

)

To link the deﬁned group pairs with the model outputs, every Data dictionary requires a

MOMap List that maps the output indices i ∈ {1, . . . , N

out

} to the respective group.

Through the use of data and discrepancy groups, it also becomes possible to address problems

where the number of measurements is not the same for each model output.

A simple example on how to use the MOMap array is given in Example 4 of the Bayesian

inversion module (04-PredPrey.ipynb).

UQ[PY]LAB-V1.0-113 - 38 -

Bayesian inference for model calibration and inverse problems

Note: The output IDs speciﬁed in the MOMap list are not unique. By giving the same

index in the MOMap vectors of different data groups, it becomes possible to cali-

brate the same computational model using measurements gathered in different

experiments.

2.4 Advanced options: solver

In the previous sections, the solver was not explicitly speciﬁed. By default UQ[PY]LAB uses

an MCMC algorithm with the afﬁne invariant ensemble sampler (see Section 1.3.4) with

a = 2.

The Bayesian module offers four different solvers that are described in more detail next:

• MCMC: Section 2.4.1

• SLE: Section 2.4.2

• SSLE: Section 2.4.3

• None: Section 2.4.4

2.4.1 MCMC

Currently, four MCMC samplers are shipped with the inversion module of UQ[PY]LAB:

Metropolis Hastings (MH), Adaptive-Metropolis (AM), Hamiltonian Monte Carlo (HMC) and

afﬁne invariant ensemble sampler (AIES). Their theoretical foundations are detailed in Sec-

tion 1.3. An exhaustive list of all available options is given in Section 3.1.4.

To select an MCMC sampler, the following keys have to be speciﬁed:

Solver = {

'Type': 'MCMC',

'MCMC': {

'Sampler': 'MH', # AM, HMC, AIES

}

The number of iterations done by the sampler is given as a scalar in the Steps key:

Solver['MCMC']['Steps'] = 200

Note that the cost per iteration depends on the speciﬁc sampler (see Section 1.3). All MCMC

samplers require a set of initial seeds for the individual chains. These seeds are speciﬁed

through an M × C List of lists Seed, where C is the number of desired parallel chains:

Solver['MCMC']['Seed'] = Seed

UQ[PY]LAB-V1.0-113 - 39 -

UQ[PY]LAB user manual

Alternatively it is also possible to just pass the number of chains C to the solver by passing a

scalar value to the NChains key. For C = 20 this can be done by:

Solver['MCMC']['NChains'] = 20

The Bayesian module then automatically samples seeds from the prior distribution π(x).

The key Sampler speciﬁes the sampling algorithm. The options are MH (Metropolis-Hastings,

Section 1.3.1), AM (adaptive Metropolis, Section 1.3.2), HMC (Hamiltonian Monte Carlo, Sec-

tion 1.3.3), and AIES (afﬁne invariant ensemble sampler, Section 1.3.4). Depending on the

sampler, different options can be speciﬁed that are discussed in more detail next.

2.4.1.1 Metropolis-Hastings algorithm

The only parameter of the Metropolis-Hastings algorithm is the proposal distribution (see

Section 1.3.1 and Table 12), which can be speciﬁed by deﬁning a myProposal dictionary

(see also Table 14).

If the algorithm is to use a Gaussian proposal distribution centered at the previous sample

(standard random walk algorithm), the myProposal dictionary should only contain one key

PriorScale that can, for instance, be set to:

myProposal = {'PriorScale': 0.1}

This PriorScale is then used to deﬁne a covariance matrix Σ

as a diagonal matrix propor-

tional to the M prior marginal variances. Alternatively this covariance matrix can also be

fully speciﬁed:

myProposal = {'Cov': Cov}

where Cov is an M × M positive-deﬁnite matrix.

Alternatively, if more advanced (e.g. non-Gaussian) proposal distributions are required, the

myProposal dictionary can also contain two keys:

myProposal = {

'Distribution': myProposalDistribution['Name'],

'Conditioning': 'Previous' # Other valid option: 'Global'

}

where myProposalDistribution is a UQ[PY]LAB INPUT object (UQ[PY]LAB User Manual

– the INPUT module). The key Conditioning can either be set to 'Global' (default), to

draw proposals from the distribution speciﬁed by myProposalDistribution independently

of the previous sample point, or to 'Previous', which sets the mean value of the proposal

distribution to the previous sample point in every step.

UQ[PY]LAB-V1.0-113 - 40 -

Bayesian inference for model calibration and inverse problems

Note: A proposal distribution speciﬁed with the 'Previous' option can be very slow

compared to 'Global' proposals. Unless there is a justiﬁed reason to use this

option it is thus not recommended.

Finally, the proposal dictionary needs to be assigned to the Solver['MCMC']['Proposal']

key:

Solver['MCMC']['Proposal'] = myProposal

2.4.1.2 Adaptive Metropolis algorithm

The adaptive Metropolis algorithm takes the keys Proposal, T0, and Epsilon (see also

Table 13). They can be, for example, set to:

SolverMCMC = {

'Proposal': myProposal,

'MCMC': {

'T0': 1e2,

'Epsilon': 1e-4

}

where the myProposal dictionary can be speciﬁed as detailed in Section 2.4.1.1 and T0 is

the number of iterations t

during which the sampler uses the supplied proposal distribu-

tion myProposal before switching to the Gaussian distribution with the empirical covariance

matrix as detailed in Eq. (1.34). The small number Epsilon speciﬁes the ϵ added to the

empirical correlation matrix to avoid singularity (see Section 1.3.2). If it is not speciﬁed, it is

automatically set to ϵ = 10

−6

2.4.1.3 Hamiltonian Monte Carlo algorithm

The parameters of the Hamiltonian Monte Carlo algorithm are (see also Table 15):

SolverMCMC = {

'MCMC': {

'LeapfrogStep': 40,

'LeapfrogSize': 0.1,

'Mass': 1

}

with the number of leapfrog steps LeapfrogSteps (N

), the leapfrog step size LeapfrogSize

(

), and the mass matrix M (see Section 1.3.3).

The mass matrix can be passed as a scalar value m, in which case the mass matrix takes the

form M = mI

or directly as an M × M matrix.

UQ[PY]LAB-V1.0-113 - 41 -

UQ[PY]LAB user manual

2.4.1.4 Afﬁne invariant ensemble algorithm

The only parameter of the afﬁne invariant ensemble algorithm is the scalar a (see also Ta-

ble 16). It can be set to any scalar parameter a > 1, for example for a = 3:

SolverMCMC = {'MCMC': {'a': 3}}

This deﬁnes the parameter used for the stretch move distribution in Eq. (1.43). If this param-

eter is not set, a value of a = 2 is assumed in accordance with Goodman and Weare (2010);

Allison and Dunkley (2013); Wicaksono (2017).

2.4.1.5 Visualization

UQ[PY]LAB currently does not offer the option to enable live trace plots during runtime.

However, it is possible to display trace plots to assess convergence after the analysis is per-

formed by using the following command:

uq.display(myBayesianAnalysis,trace='all')

This command displays all available traces. You can also specify desired traces using the

trace parameter. For instance, you can use trace=1 to display only the trace for the ﬁrst

marginal or trace=[1,3] to display traces for the ﬁrst and the third marginals. The resulting

plot will be similar to the one shown in Figure 8.

Figure 8: Trace plot and corresponding KDE after execution of the MCMC algorithm.

2.4.1.6 Post-processing

Following the analysis, the sample points generated by the MCMC algorithm are stored in

the Results key of the myBayesianAnalysis dictionary:

{

UQ[PY]LAB-V1.0-113 - 42 -

Bayesian inference for model calibration and inverse problems

'Sample': [[[...], [...]], [[...], [...]], ..., [[...], [...]]],

'Acceptance': [...],

'Time': ...,

'ForwardModel': {'evaluation': [[...], [...], ..., [...]]},

'LogLikeliEval': [[...], [...], ..., [...]],

'PostProc': {...}

}

The sample points are stored in the 3D list Sample with a regular shape T × M × C. The

associated forward model and log likelihood evaluations are stored in the ForwardModel

dictionary and LogLikeliEval list of lists with a regular shape T × C respectively.

Sample points generated by MCMC algorithms typically require post-processing before they

can be used as a true posterior sample. In the Bayesian module of UQ[PY]LAB this post-

processing is automatically done with the uq.postProcessInversionMCMC function that is

called for MCMC analyses by the wrapper function:

myBayesianAnalysis = uq.postProcessInversion(myBayesianAnalysis)

This function is called automatically after every analysis and performs a set of default post-

processing procedures: the ﬁrst half of all sample points generated by the MCMC chains are

removed, the empirical parameter mean is estimated from these remaining sample points

along with the 2.5th and 97.5th percentiles, the covariance matrix is estimated, and samples

are drawn from the prior distribution and posterior predictive distribution.

These post-processing results are then stored in the PostProc dictionary inside the Results

dictionary. If uq.postProcessInversion is called with all possible options (see Sec-

tion 3.3), it contains the following keys of myBayesianAnalysis['Results']['PostProc']

dictionary

{

'PostSample': [[[...], [...]], [[...], [...]], ..., [[...], [...]]],

'PostLogLikeliEval': [[...], [...], ..., [...]],

'PostModel': {'evaluation': [[...], [...], ..., [...]]},

'PointEstimate': {'ForwardRun': {...}, 'X': [...], 'Type': ...},

'Dependence': {'Corr': [[...],[...]], 'Cov': [[...],[...]]},

'Percentiles': {'Values': [[...],[...]], 'Probabilities': [...],

'Mean': [...], 'Var': [...]},

'PriorSample': [[...], [...], ...,[...]],

'PostPredSample': {'ModelEvaluations': [...], 'Sample': [...],

'Discrepancy': [...]},

'ChainsQuality': {'BadChains': [...], 'GoodChains': [...]},

'MPSRF': ...,

'PriorPredSample': {'ModelEvaluations': ..., 'Sample': ...,

'Discrepancy': ...}

}

where the value of the PostSample key is a regular nested list of size T"xMxC" and type ﬂoat,

where T" is the length of the MCMC chains without the burn in, C" is the number of chains

excluding the badChains and P is the number of drawn prior sample points. The value of

UQ[PY]LAB-V1.0-113 - 43 -

UQ[PY]LAB user manual

the PostLogLikeliEval key is a regular nested list of size T"xC" and type ﬂoat.

For more control over the post-processing operations, it is recommended to adapt the post-

processing options to the problem at hand and recall the uq.postProcessInversion func-

tion after the analysis is complete. For example, to draw 1,000 sample points from the prior

predictive distribution, the function needs to be called with:

myBayesianAnalysis = uq.postProcessInversion(myBayesianAnalysis,

'priorPredictive', 1000)

uq.postProcessInversion always operates on the original Sample list and overwrites the

contents of the respective PostProc key when called repeatedly. As an example, after calling

myBayesianAnalysis = uq.postProcessInversion(myBayesianAnalysis,

'burnIn', 0.6)

T" will be 40% of the original T . After calling

myBayesianAnalysis = uq.postProcessInversion(myBayesianAnalysis,

'burnIn', 0.7)

T" will be 30% of the original T , independent of the previous call.

To maintain previously computed content of PostProc’s keys, uq.postProcessInversion

only overwrites values of the keys if the call produces new content or an output is explicitly

turned off by the user. This does not happen with the default options of uq.postProcessInversion

(see Section 3.3). In order to delete a previously computed post-processing result, it needs

to be removed explicitly e.g.:

myBayesianAnalysis = uq.postProcessInversion(myBayesianAnalysis,

'priorPredictive', 0)

removes a previously computed prior predictive sample while any other call would leave it

untouched.

An exhaustive list of available post-processing options can be found in Section 3.3.

2.4.2 Spectral likelihood expansion

To solve the Bayesian problem with the spectral likelihood expansion (SLE) technique (see

Section 1.4), the user should specify the Solver dictionary as follows:

Solver = {'Type': 'SLE'}

The Bayesian module of UQ[PY]LAB uses the polynomial chaos expansion module (PCE,

UQ[PY]LAB User Manual – Polynomial Chaos Expansions) to construct the likelihood approx-

imation. Additional options for the PCE can be passed to UQ[PY]LAB through the following

syntax:

UQ[PY]LAB-V1.0-113 - 44 -

Bayesian inference for model calibration and inverse problems

Solver['SLE'] = {

'Degree': list(range(1,21)),

'ExpDesign': {

'NSamples': 1e4

}

which speciﬁes a degree adaptive PCE to be computed based on an experimental design of

size 10

. In general, SLE accepts all options that can be passed to the MetaOpts dictionary in

UQ[PY]LAB User Manual – Polynomial Chaos Expansions (e.g., TruncOptions). See Table 18

for more information.

To conduct an SLE-based inverse analysis, the following command should be executed (after

assigning the Solver dictionary to the BayesOpts['Solver'] key):

uq.createAnalysis(BayesOpts)

By default, the SLE-speciﬁc post-processing function uq.postProcessInversionSLE is

called after every SLE analysis by the wrapper function uq.postProcessInversion: it

computes the Bayesian evidence (see Eq. (1.52)) and the posterior mean (see Eq. (1.56)).

These values are stored in the myBayesianAnalysis['Results']['PostProc'] key that

reads

{

'Evidence': ...,

'Posterior': '...',

'Mean': [..., ...],

'PointEstimate': {'X': [..., ...], 'Type': '...'}

}

To additionally compute the posterior covariance and correlation matrices, the post-processing

function can be called again with the following arguments:

myBayesianAnalysis = uq.postProcessInversion(myBayesianAnalysis,

'dependence', True)

The results of the analysis can be inspected with the uq.print function.

2.4.3 Stochastic spectral likelihood embedding

To solve the Bayesian problem with the stochastic spectral likelihood embedding (SSLE)

technique (see Section 1.5), the user should specify the Solver dictionary as follows:

Solver = {'Type': 'SSLE'}

The Bayesian module of UQ[PY]LAB uses the stochastic spectral embedding module (SSE,

UQ[PY]LAB User Manual – Stochastic spectral embedding) to construct the likelihood ap-

proximation. Additional options for the SSE can be passed to UQ[PY]LAB through the

UQ[PY]LAB-V1.0-113 - 45 -

UQ[PY]LAB user manual

following syntax:

Solver['SSLE'] = {

# Expansion options

'ExpOptions': {

'Degree': list(range(1,5))

# Experimental design options

'ExpDesign': {

'NSamples': 1000,

'NEnrich': 100

}

which speciﬁes degree adaptive PCE for the residual expansions with an experimental design

of size 1000. By default, the experimental design is added sequentially (see Section 1.5.1.2),

where we now speciﬁed an enrichment rate of 100 sample points per reﬁnement step.

By default, SSLE uses a partitioning strategy based on the residual variance difference (see

Section 1.5.1.1)

uq.createAnalysis(BayesOpts)

By default, the SSLE-speciﬁc post-processing function uq.postProcessInversionSSLE is

called after every SLE analysis by the wrapper function uq.postProcessInversion: it

computes the Bayesian evidence (see Eq. (1.59)) and the posterior mean (see Eq. (1.63)).

These values are stored in the myBayesianAnalysis['Results']['PostProc'] key:

{

'Evidence': ...,

'Mean': [..., ...],

'PointEstimate': {'X': [..., ...], 'Type': '...'},

'Posterior': '...',

}

To additionally compute the posterior covariance and correlation matrices, the post-processing

function can be called again with the following arguments:

myBayesianAnalysis = uq.postProcessInversion(myBayesianAnalysis,

'dependence', True).

The results of the analysis can be inspected with the uq.print function.

2.4.4 No solver: posterior point by point evaluation

Sometimes it is not required to solve the inverse problem, but only to evaluate the prior/-

posterior PDFs, or the likelihood function for speciﬁc parameter points x

(e.g., maximum a

posteriori estimation). In this case, the solver type has to be set to 'None':

Solver['Type'] = 'None'

UQ[PY]LAB-V1.0-113 - 46 -

Bayesian inference for model calibration and inverse problems

After the analysis creation with uq.createAnalysis(), the analysis object contains the

following:

{

'results_idx': 2,

'Options': {'Data': {...}, 'Solver': {...}, 'Type': 'Inversion'},

'Results': 0,

'Internal': {'customLikeli': 0, ...},

'Name': 'Analysis 1',

'Type': 'uq_inversion',

'core_component': 'analysis',

'displayFun': 'None',

'printFun': 'None',

'Discrepancy': {'Type': 'Gaussian', ...},

'Data': {'Name': 'Beam mid-span deflection', 'y': [...], 'MOMap': [...]},

'PriorDist': 'Input 2',

'ForwardModel': {'Model': 'Model 1', 'PMap': [...]},

'LogPrior': '@(x)uq_evalLogPDF(x,Internal.FullPrior)',

'UnnormPosterior': '...',

'UnnormLogPosterior': '...',

'Class': 'uq_analysis',

'Prior': '@(x)uq_evalPDF(x,Internal.FullPrior)',

'Likelihood': "@(x)uq_inversion_likelihood(x,Internal,'Likelihood')",

'LogLikelihood': "@(x)uq_inversion_likelihood(x,Internal,'LogLikelihood')"

}

If the solver type is set to 'None', no results are generated. Instead the analysis generates

function handles to the prior distribution, the likelihood function and the unnormalized pos-

terior distribution. The handle to the unnormalized posterior distribution can be used to ﬁnd

the maximum a posteriori (MAP) parameter value as shown in the supplied Example 7 of

Bayesian inversion module (07-MAP.ipynb).

2.5 Advanced feature: multiple forward models

When data from multiple data sources are available (e.g., stresses and temperatures), dif-

ferent computational models may be needed. In the Bayesian module of UQ[PY]LAB, it is

possible to perform an inversion on multiple models with different output and discrepancy

options that depend on the same set, or subsets, of the parameters x

. This can be achieved

by specifying dictionaries for the ForwardModel key.

To show how such a problem is set up in UQ[PY]LAB, it is assumed now that additional

longitudinal tensile tests are carried out on a beam as shown in Figure 9.

The nominal value of the load P is 50 kN. In total N = 3 tensile tests are carried out,

resulting in the deformations given in Table 4.

Table 4: Beam elongation: measured deformations.

Experiment 1 2 3

U (mm) 0.235 0.236 0.229

UQ[PY]LAB-V1.0-113 - 47 -

UQ[PY]LAB user manual

Figure 9: Tensile test.

The elongation U of the specimen under a load P can be computed with (Sudret, 2018):

U =

P L

Ebh

. (2.4)

In this case, the measurements were carried out with a measuring device that has a known

measurement error distribution of ε ∼ N (0, 2 · 10

−11

) (m).

2.5.1 Specify a prior distribution

The model prior object from Section 2.2.2 can be extended to contain the point load P :

PriorOpts['Marginals'].append(

{

'Name': 'P',

'Type': 'Constant',

'Parameters': [50000] # (N)

}

)

myPriorDist = uq.createInput(PriorOpts)

2.5.2 Create a forward model

Each forward model has to be set up as a dedicated UQ[PY]LAB MODEL. To do this, Eq. (2.4)

is translated to a string UQ[PY]LAB MODEL and assigned together with the

uq.SimplySupportedBeam model:

# Forward model 1

ModelOpts1 = {

'Type' : 'Model',

'Name': 'Beam bending deflection',

'ModelFun':'SimplySupportedBeam.model',

}

myModel1 = uq.createModel(ModelOpts1)

# Forward model 2

ModelOpts2 = {

'Type' : 'Model',

'Name': 'Beam elongation',

'mString': 'X(:,5).

X(:,3)./(X(:,1).

X(:,2).

X(:,4))',

UQ[PY]LAB-V1.0-113 - 48 -

Bayesian inference for model calibration and inverse problems

}

myModel2 = uq.createModel(ModelOpts2)

# Forward models

ForwardModels = [

{

'Model': myModel1['Name'],

'PMap': [1,2,3,4,5]

{

'Model': myModel2['Name'],

'PMap': [1,2,3,4,6]

}

]

The parameter map PMap deﬁnes which parameters from the model parameter vector x

are used for the respective model. The tensile test in ForwardModels[1] takes the 6-th

parameter from the model prior distribution as a 5-th input ('X(:,5)' in the equation string

takes the 6-th parameter from the model prior distribution as a 5-th input ('X(:,5)' in the

equation string; note that mString follows the Matlab syntax).

Note: By default, the PMap list is set to address each parameter in the model parameter

list x

. However, in most realistic usage scenarios, it should be updated to list

the desired parameters properly.

2.5.3 Provide measurements

The results of the tensile test are stored in a list U and are assigned to the second data group

myData[1]:

U = np.array([0.235, 0.236, 0.229])/1000

myData = [

# Data group 1

{

'y': V_mid.tolist(),

'Name': 'Beam mid-span deflection',

'MOMap': [1, # Model ID

1] # Output ID

# Data group 2

{

'y': U.tolist(),

'Name': 'Beam elongation',

'MOMap': [2, # Model ID

1] # Output ID

}

]

UQ[PY]LAB-V1.0-113 - 49 -

UQ[PY]LAB user manual

Note: In the case of multiple models, the MOMap is mandatory and has to be deﬁned as

a list if each model has only one output or a nested list in case of several outputs

per model. In the second case scenario, each inner list contains an index of the

corresponding model and an index of the corresponding model output.

2.5.4 Deﬁne a discrepancy model

The discrepancy can then be speciﬁed separately for the ﬁrst model M

), where the

variance of the discrepancy term is inferred from data, and the second model M

where the measurement error variance is known, as follows:

# Prior options for Discrepancy group 1

DiscrepancyPriorOpts1 = {

'Name': 'Prior of sigma_1ˆ2',

'Marginals': {

'Name': 'Sigma2',

'Type': 'Uniform',

'Parameters': [0, 1e-4], #(mˆ2)

}

myDiscrepancyPrior1 = uq.createInput(DiscrepancyPriorOpts1)

DiscrepancyOpts = [

# Discrepancy group 1

{

'Type': 'Gaussian',

'Prior': myDiscrepancyPrior1['Name']

# Discrepancy group 2

{

'Type': 'Gaussian',

'Parameters': [2e-11], # (mˆ2) known discr. variance

}

]

2.5.5 Perform the inverse analysis

As in the previous case, all the options are then gathered in a single dictionary that contains

all the information to perform the Bayesian inversion:

BayesOpts = {

'Type': 'Inversion',

'Name': 'Bayesian multiple models',

'Prior': myPriorDist['Name'],

'ForwardModel': ForwardModels,

'Data': myData,

'Discrepancy': DiscrepancyOpts,

}

In this case, PriorDist has to be explicitly assigned to the BayesOpt dictionaries to avoid

confusion with the INPUT object in the discrepancy model. It is also necessary to explicitly

UQ[PY]LAB-V1.0-113 - 50 -

Bayesian inference for model calibration and inverse problems

assign the ForwardModels list of dictionaries, as there are multiple forward models now.

The inverse problem can then be solved by running the analysis:

myBayesianAnalysis = uq.createAnalysis(BayesOpts)

Run post-processing with non-default options by removing 30%, instead of the default 50%,

of the initially generated sample points:

myBayesianAnalysis = uq.postProcessInversion(myBayesianAnalysis,

'burnIn', 0.3)

The results of this analysis can be assessed with a brief report generated by:

uq.print(myBayesianAnalysis)

which returns:

%----------------------- Inversion output -----------------------%

Number of calibrated model parameters: 1

Number of non-calibrated model parameters: 5

Number of calibrated discrepancy parameters: 1

%------------------- Data and Discrepancy

% Data-/Discrepancy group 1:

Number of independent observations: 5

Discrepancy:

Type: Gaussian

Discrepancy family: Scalar

Discrepancy parameters known: No

Associated outputs:

Model 1:

Output dimensions: 1

% Data-/Discrepancy group 2:

Number of independent observations: 3

Discrepancy:

Type: Gaussian

Discrepancy family: Scalar

Discrepancy parameters known: Yes

Associated outputs:

Model 2:

Output dimensions: 1

%------------------- Solver

Solution method: MCMC

Algorithm: AIES

Duration (HH:MM:SS): 00:00:53

Number of sample points: 3.00e+04

%------------------- Posterior Marginals

---------------------------------------------------------------------

UQ[PY]LAB-V1.0-113 - 51 -

UQ[PY]LAB user manual

---------------------------------------------------------------------

| E | 2.4e+10 | 2.6e+08 | (2.3e+10 - 2.4e+10) | Model |

| Sigma2 | 2e-06 | 4.9e-06 | (1.4e-07 - 1.4e-05) | Discrepancy |

---------------------------------------------------------------------

%------------------- Point estimate

----------------------------------------

| Parameter | Mean | Parameter Type |

----------------------------------------

| E | 2.4e+10 | Model |

| Sigma2 | 2e-06 | Discrepancy |

----------------------------------------

A visualization of the analysis can be produced by:

uq.display(myBayesianAnalysis)

which produces the images shown in Figure 10.

A simple example on how to use the multiple forward model feature is given in Example 5 of

the Bayesian inversion module 05-MultipleModels.ipynb.

UQ[PY]LAB-V1.0-113 - 52 -

Bayesian inference for model calibration and inverse problems

(a) Scatterplots of prior and posterior sample with the posterior mean point estimator from (1.31).

(b) Posterior predictive distribution from (1.28), data, and model at mean prediction obtained by propa-

gating the mean point estimator through the forward model for both data groups.

Figure 10: Multiple forward models: The results of the Bayesian inverse analysis on the input

and the model predictions.

UQ[PY]LAB-V1.0-113 - 53 -

Chapter 3

Reference List

How to read the reference list

Python dictionaries play an important role throughout the UQ[PY]LAB syntax. They offer

a natural way to semantically group conﬁguration options and output quantities. Due to

the complexity of the algorithms implemented, it is not uncommon to employ nested dictio-

naries to ﬁne-tune the inputs and outputs. Throughout this reference guide, a table-based

description of the conﬁguration dictionaries is adopted.

The simplest case is given when a value of a dictionary key is a simple value or a list:

Table X: Input

Name String A description of the ﬁeld is put here

which corresponds to the following syntax:

Input = {

'Name' : 'My Input'

}

The columns, from left to right, correspond to the name, the data type and a brief description

of each key-value pair. At the beginning of each row a symbol is given to inform as to whether

the corresponding key is mandatory, optional, mutually exclusive, etc. The comprehensive

list of symbols is given in the following table:

Mandatory

□ Optional

UQ[PY]LAB user manual

⊕ Mandatory, mutually exclusive (only one of

the keys can be set)

⊞ Optional, mutually exclusive (one of them

can be set, if at least one of the group is set,

otherwise none is necessary)

When the value of one of the keys of a dictionary is a dictionary itself, a link to a table that

describes the structure of that nested dictionary is provided, as in the case of the Options

key in the following example:

Table X: Input

Name String Description

□ Options Table Y Description of the Options

dictionary

Table Y: Input['Options']

Key1 String Description of Key1

□ Key2 Double Description of Key2

In some cases, an option value gives the possibility to deﬁne further options related to that

value. The general syntax would be:

Input = {

'Option1' : 'VALUE1',

'VALUE1' : {

'Val1Opt1' : ... ,

'Val1Opt2' : ...

}

This is illustrated as follows:

Table X: Input

Option1 String Short description

'VALUE1' Description of 'VALUE1'

'VALUE2' Description of 'VALUE2'

⊞ VALUE1 Table Y Options for 'VALUE1'

⊞ VALUE2 Table Z Options for 'VALUE2'

Table Y: Input['VALUE1']

□ Val1Opt1 String Description

□ Val1Opt2 Float Description

UQ[PY]LAB-V1.0-113 - 56 -

Bayesian inference for model calibration and inverse problems

Table Z: Input['VALUE2']

□ Val2Opt1 String Description

□ Val2Opt2 Float Description

UQ[PY]LAB-V1.0-113 - 57 -

UQ[PY]LAB user manual

3.1 Create a Bayesian inverse analysis

Syntax

myBayesianAnalysis = uq.createAnalysis(BayesOpts)

Input

The dictionary BayesOpts contains the information for a Bayesian inverse analysis. N

mod

the number of computational forward models and N

is the number of data- and discrepancy

groups.

Table 5: BayesOpts

Type 'Inversion' Inverse modelling

Data Table 6 The data used for inversion.