Bayesian inference for the stochastic identification of elastoplastic material
parameters: Introduction, misconceptions and insights
H. Rappela,b , L.A.A. Beexa , J.S. Halea , S.P.A. Bordasa,c,d,∗
a Faculty

arXiv:1606.02422v4 [cs.CE] 15 Jan 2017

of Science, Technology and Communication, University of Luxembourg, Campus Kirchberg, 6, Rue
Coudenhove-Kalergi, L-1359 Luxembourg, Luxembourg.
b Computational & Multiscale Mechanics of Materials (CM3), Department of Aerospace and Mechanical Engineering,
University of Liège, Quartier Polytech 1, Allée de la Découverte 9, B-4000 Liège, Belgium.
c School of Engineering, Cardiff University, Queens Buildings, The Parade, Cardiff CF243AA, Wales, UK.
d Adjunct Professor, Intelligent Systems for Medicine Laboratory School of Mechanical and Chemical Engineering, The
University of Western Australia, 35 Stirling Highway, Crawley/Perth WA 6009, Australia.

Abstract
We discuss Bayesian inference (BI) for the probabilistic identification of material parameters. This contribution aims to shed light on the use of BI for the identification of elastoplastic material parameters. For
this purpose a single spring is considered, for which the stress-strain curves are artificially created. Besides
offering a didactic introduction to BI, this paper proposes an approach to incorporate statistical errors both
in the measured stresses, and in the measured strains. It is assumed that the uncertainty is only due to
measurement errors and the material is homogeneous. Furthermore, a number of possible misconceptions on
BI are highlighted based on the purely elastic case.
Keywords: Bayesian inference, Bayes’ theorem, stochastic identification, statistical identification,
parameter identification, elastoplasticity, plasticity

1. Introduction
The most commonly used approach to identify material parameters in mechanics is to formulate an error
function that measures the difference between the response of a model and experimental data [1]. This error
functional is then minimised with respect to the material parameters of interest to calculate them. Such an
approach provides a deterministic estimate of the unknown material parameters, unable to account for the
unavoidable uncertainties associated with any experimental observation. Accounting for such uncertainties
is of critical importance to all engineering disciplines. An alternative, and rather different approach is to use
Bayesian inference (BI). Using Bayes’ theorem a probability density function (PDF), the so-called posterior
distribution (or simply posterior), can be formulated as a function of the material parameters of interest.
Subsequently, the PDF is analysed to determine relevant values, such as the mean of the material parameters,
the material properties at which the PDF is maximum and the covariance. Only for simplistic cases can the
PDF be explored analytically. Hence, numerical methods must generally be employed, e.g. Markov chain
Monte Carlo (MCMC) techniques [2–5]. An alternative is to first approximate the PDF (e.g. by a Laplace
approximation) and then analyse it[6, 7].
In contrast to conventional error minimisation approaches, BI incorporates the statistical noises of experimental devices [1]. Consequently, the identified material parameters come with a level of uncertainty.
Additionally, BI, also provides an intrinsic statistical regularisation which makes inverse problems with
limited observations solvable [8]. On the other hand, applying Bayes’ theorem for material parameter identification does require the measurement noise to be known, i.e. the noise distributions and the parameters
∗ Corresponding

author
Email address: stephane.bordas@uni.lu (S.P.A. Bordas)

1

of these distributions must be established. The numerical techniques to analyse the PDF may furthermore
need careful attention.
The developments of BI in the field of parameter identification for mechanical models started with the
identification of elastic constants. Isenberg [9] proposed a Bayesian approach for the identification of elastic
parameters in 1979. This framework was subsequently used by various researchers to identify elastic material
parameters based on dynamic responses [6, 10, 11]. Lai and Ip [12] used BI to identify the elastic properties
of a thin composite plate. Daghia et al. [13] used the Bayesian framework for the dynamic identification
of the elastic constants of thick laminated composite plates. Koutsourelakis [14] used Bayesian inference
to identify spatially varying elastic material parameters. In 2010 Gogu et al. [1] presented an introduction
in the Bayesian approach for the identification of elastic constants, and compared the results with those
of error minimisation. The influence of the prior distributions was however not systematically studied. In
another study, Gogu et al. [15] used the Bayesian framework to identify elastic constants in multi-directional
laminates.
BI is also used for parameter identification of nonlinear constitutive models. Muto and Beck [16] and Liu
and Au [17] used the approach for hysteretic models, whereas Fitzenz et al. [18] used BI for a creep model of
quartz. Most [19] used the Bayesian updating procedure for the parameter identification of an elastoplastic
model without hardening (perfect plasticity). Rosić et al. [20] used linear Bayesian updating via polynomial
chaos expansion for an elastoplastic system. Hernandez et al. [21] used BI for viscoelastic materials.
Another study that uses Bayes’ theorem to identify material parameters is the work of Nichols et al. [22].
They employed the Bayesian approach to identify the nonlinear stiffness of a dynamic system. Furthermore,
Nichols et al. [22] used the method to find the location, size and depth of delamination in a composite beam.
Abhinav and Manohar [23] used BI to characterise the dynamic parameters of a structural system with
geometrical nonlinearities. The approach is also employed as a framework to assess the quality of different
models with respect to measured data (i.e. model selection): e.g. hyperelastic constitutive models for soft
tissue [24], phenomenological models for tumour growth [25], models for damage progression in composites
due to fatigue [26] and fatigue models for metals [27]. Sarkar et al. [28] used the Bayesian method to identify
thermodynamical parameters of cementitious materials. BI is also used in the field of mechanics for inverse
problems employing differential equations [3, 29].
Bayesian inversion relies on concept that are not generally intuitive to those who have focused upon
deterministic inversion methods. This is particularly true for nonlinear material models. Understanding the
role of prior knowledge on the results is essential.
We have not encountered any Bayesian inversion study within the field of mechanics where errors on
strain and stress measures are simultaneously taken into account. In other fields, a few studies can be found
that deal with two error sources [30–35]. Yet, the methods described in these papers are only investigated in
the context of linear regression, which is insufficient to tackle the problems we are interested in, in particular
nonlinear hardening.
Remark. In practice displacement and force are usually the only quantities which can be directly measured
(strains and stresses are deduced from those displacement and force measurements). In this work, we focus
on a single spring. In this case, the strain is proportional to the displacement and the force to the stress.
Consequently, in the remainder of this paper, we will abusively assume that the stress and strains are the
measured quantities, which will facilitate writing the constitutive relations directly in terms of measured data.
This contribution aims to show how Bayesian inference can be applied for the stochastic identification of
elastoplastic material parameters. The important issues focused on are
(1) a formulation to incorporate the split in the purely elastic part of the response and the elastoplastic
part,
(2) a formulation to incorporate the uncertainty in the measured stresses as well as the uncertainty in the
measured strains for elastic(-plastic) models and
(3) the presentation of possible misconceptions on BI. This is important to help newcomers to the field.
The elastic case is used as an illustration for this, because it is the simplest repetitive case to grasp.
2

Our contribution focuses on one dimensional stress-strain results of a single spring in order to be as
straightforward as possible. This single spring can be seen as an isotropic material model with no influence
of the Poisson’s ratio. We also propose a Bayesian framework to not only take the statistical errors in the
stress measurements into account, but also the statistical error in the strain measurements. The latter part
shows similarities with the study of Kelly [34] in 2007, by computing the likelihood function considering all
measured data, and the marginalisation integral.
The structure of the paper is as follows. Section 2 discusses the theoretical fundamentals behind Bayes’
theorem, the elastoplastic material models and MCMC methods as the numerical techniques to analyse the
posterior distribution. Section 3 describes a Bayesian approach for the stochastic identification of elastoplastic material parameters, when only the stress measurements include stochastic errors. Section 3 presents
how to split the experimental data in a purely elastic domain and an elastoplastic domain. In section 4
the Bayesian method to incorporate both the statistical errors in the stresses and strains is presented. In
section 5 a considerable number of relevant examples for all the aforementioned approaches is presented.
Furthermore, in this section some misconceptions about the Bayesian method are indicated. In section 6
conclusions are finally presented.
2. Concepts
A conventional procedure to identify material parameters in mechanics relies on error minimisation.
Often, a least-squared error functional is formulated that measures the squared difference between the
experimental data and the predicted model response. This error functional is subsequently minimised with
respect to the material parameters of interest. A measure for the quality of the model response with respect to
the measured data is then formed by the residual (i.e. the final least-squared difference between the measured
data and the model response). This formulation does not take into account the statistical information.
However there exist some advanced least-squares formulations that use the statistical information [36, 37].
Identification approaches based on BI provide an alternative that aims to construct a probability distribution function (PDF) based on the measured data and expert’s prior knowledge, as a function of the
material parameters of interest. In the case of one parameter, the PDF measures the probability that a
particular value of the parameter occurs. Hence, the PDF must be analysed to investigate which values
the parameter will most probably be. Important quantities of interest are the value of parameter at which
the PDF is maximum (i.e. ‘modes’ in general statistical terminology or ‘maximum-a-posteriori-probability’,
abbreviated to ‘MAP’, in Bayesian terminology), the mean and the variance (i.e. the average of the squared
differences from the mean).
BI may be considered as an approach that ‘accounts’ for the fact that only a limited number of observations are made. This is achieved by incorporating some assumed prior knowledge about the parameters
of interest. For inverse modelling, the prior knowledge (prior distribution, or simply ‘prior’ in Bayesian
terminology) regularises the inverse problem. This ensures that the systems to be solved are not ill-posed.
Note that conventional error minimisation approaches based on least-squares, are the same as the frequentist approach based on the maximum likelihood (ML) i.e. maximising the likelihood function see sub
section 2.1 with respect to the parameter of interest, if the statistical error is additive, uncorrelated and
normally distributed with a zero standard deviation [8]. Note that in the limit of an infinite number of
measurements, the frequentist approach based on maximum likelihood and the Bayesian approach tend to
the same results because the likelihood function will become dominant in the Bayesian formulation [38].
The remainder of this section presents the essential concepts of necessary to apply BI to the identification
of elastoplastic material parameters in the case of a single spring. The concepts involve Bayes’ theorem, the
material models of interest and the adaptive MCMC method to analyse the posterior distributions.
Note that in this section capital letters denote random variables and bold ones denote vectors and
matrices.
2.1. Bayesian inference
We start by considering events A and B, and the discrete probabilities of each event: P (A) and P (B).
The probability that events A and B both occur, is given by the joint probability, P (A, B), which can be
3

expanded as:
P (A, B) = P (A|B)P (B) = P (B|A)P (A),

(1)

where P (A|B) and P (B|A) are conditional probabilities. Conditional probability P (A|B) expresses the
probability that event A occurs, if it is certain that event B occurs. Using Eq. (1), the simplest form of
Bayes’ theorem can be written as:
P (A|B) =

P (A)P (B|A)
.
P (B)

(2)

If one regards two continuous random variables X ∈ Rn and Y ∈ Rk , instead of discrete variables, where
X is a random vector with n unknown parameters and Y a random vector with k measurements, Eq. (2)
can be rewritten in terms of the following probability distribution functions (where π denotes a PDF):
π(x|y) =

π(x)π(y|x)
,
π(y)

(3)

where π(x), π(y|x) and π(x|y) are referred to as the prior distribution (i.e. the PDF that includes the prior
knowledge), the likelihood function ( i.e. the PDF of the observed data y as a function of unknown parameters
x) and the posterior distribution (i.e. the PDF of the unknown parameters x given the observations y),
respectively.
Using the law of total probabilities [39] which relates the marginal probabilities (π(x) and π(y)) to the
conditional probabilities (π(y|x)), the denominator in Eq. (3) can be written as:
Z
π(y) =
π(x)π(y|x)dx.
(4)
Rn

Since the data (y) is already measured, the denominator in Eq. (3) is a constant number, C ∈ R+ . This
constant number can be regarded as a normalisation factor that ensures that the integral of the posterior
(π(x|y)) equals 1:
1
π(x)π(y|x).
(5)
C
In the case of parameter identification, one is generally interested in the maximum-a-posteriori-probability
point (the value of the unknown parameter at which the PDF has the maximum value), the mean and the
covariance matrix, which are independent of C. Hence, it suffices to rewrite Eq. (5) as:
π(x|y) =

π(x|y) ∝ π(x)π(y|x).

(6)

In order to obtain the posterior (i.e. the PDF of the vector of unknown parameters given the observations
π(x|y)), the likelihood function (π(y|x)) and the prior (π(x)) need to be formulated. First, the likelihood
function is considered.
In order to construct the likelihood function, a noise model has to be formulated and a noise distribution
(πnoise ) has to be determined and calibrated. For the moment, we assume that the noise distribution is known
(including its parameters). The noise model used in this study is additive, which is frequently employed,
amongst others in several inverse studies [8, 13]. The additive noise model can be written as follows:
Y = f(X) + Ω,
n

(7)
k

where X ∈ R is again the vector with unknown material parameters, Y ∈ R the vector with the measured
data and Ω ∈ Rk the noise vector. f : Rn → Rk denotes the model and is a function of the unknown material
parameters (X). Given the realisations X = x and Ω = ω, and assuming that the parameters (X) and the
error (Ω) are statistically independent, the likelihood function reads:
π(y|x) = πnoise (y − f(x)),
4

(8)

where πnoise (ω) is the PDF of the noise (which is assumed to be identified based on separate calibration
experiments, see subsections 3.1 and 4.1). Substitution of Eq. (8) in Eq. (6) yields:
π(x|y) ∝ π(x)πnoise (y − f(x)).

(9)

A critical aspect of the Bayesian framework is the selection of the prior distribution (π(x)) [8] in which
a-priori knowledge about the parameters is translated in terms of a PDF. If enough data is observed the
prior hardly has an influence [24]. The influence of the prior distribution will be considered in more detail
in section 5.
Once the posterior is formulated (Eq. (9)), the mean parameter values, MAP parameter values and the
covariance matrix have to be extracted from it in an analytical or numerical manner (e.g. using Markov
chain Monte Carlo methods). The adaptive Markov chain Monte Carlo is considered in subsection 2.3, as it
is the numerical approach employed in this contribution to evaluate the posteriors.
2.2. Material models
In this contribution, BI is developed to identify the parameters of four one-dimensional material models:
linear elasticity, linear elasticity with perfect plasticity, linear elasticity with linear hardening and linear
elasticity with nonlinear hardening. Hardening is considered to be kinematic and associative. For each
model, the identification is based on the results of monotonic uniaxial tensile tests. Below, the governing
equations of the models are given, as well as those for when monotonic tensile loading takes place.
2.2.1. Linear elastic
The linear elastic model assumes a linear relationship between the stress and the strains. In the one
dimensional case, this writes:
σ(, x) = E,

(10)

where σ is the stress,  is the strain, x is the material parameter vector (here x = E) and E is Young’s
modulus and assumed to be constant in the domain.
2.2.2. Linear elastic-perfectly plastic
The linear elastic-perfectly plastic model neglects the effect of work hardening, assuming that purely
plastic deformation occurs when the stress reaches its yield value. The total strain () in this contribution
is additively split into an elastic part e and a plastic part p :
 = e + p ,

(11)

and the stress is defined as a function of the elastic strain, e :
σ(, x) = Ee = E( − p ).

(12)

The yield condition at which plastic yielding occurs, is written as:
f (σ) = |σ| − σy0 ≤ 0,

where σy0 is the initial yield stress and f is the yield function. Consequently, x = E
Furthermore, the flow rule for the plastic strain can be written as:

(13)
σy0

T

.

∂f
= α̇ sgn(σ),
(14)
∂σ
where sgn(·) is the sign function and α is the cumulative plastic strain. Finally, the Kuhn-Tucker conditions
[40] must hold:
˙p = α̇

α̇ ≥ 0,

f (σ) ≤ 0, α̇f (σ) = 0.
5

(15)

The stress-strain response of the linear elastic-perfectly plastic model during monotonic tension can be
written as:
(
σ
E
if  ≤ Ey0
.
(16)
σ=
σ
σy0 if  > Ey0
Using the Heaviside step function (h(·)), Eq. (16) can be conveniently expressed as:




σy0 
σy0 
σ(, x) = E 1 − h  −
+ σy0 h  −
.
E
E

(17)

2.2.3. Linear elastic-linear hardening
The linear elastic-linear hardening model is identical to the linear elastic-perfectly plastic model, except
for the yield function, which writes:
f (σ) = |σ| − σy0 − Hα ≤ 0,

T
where H is the plastic modulus. Here, x = E σy0 H .
Consequently, the stress-strain response of the model during monotonic tension writes:
(
σ
E
if  ≤ Ey0
σ=
,
σ
σy0 + Hp if  > Ey0
which can again be expressed using the Heaviside step function:

 


σy0 
HE 
σy0  
σy0 
σ(, x) = E 1 − h  −
+ σy0 +
−
h −
.
E
H +E
E
E

(18)

(19)

(20)

2.2.4. Linear elastic-nonlinear hardening
The linear elastic-nonlinear hardening model only differs from the linear elastic-perfectly plastic model
through the yield function which writes:
f (σ) = |σ| − σy0 − Hαn ≤ 0,

where n is an additional plastic material parameter and hence, x = E σy0 H
For the monotonic loading case, the stress-strain response can be written as:
(
σ
E
if  ≤ Ey0
,
σ=
σ
σy0 + Hnp if  > Ey0

(21)
n

T

.

or using the Heaviside step function:

 



σy0 
σ(, x) n 
σy0 
σ(, x) = E 1 − h  −
+ σy0 + H  −
h −
.
E
E
E

(22)

(23)

It is important to note that Eq. (23) is an implicit function of the stress (σ(, x) appears both on the left
hand side and right hand side of Eq. (23) and cannot immediately be identified from Eq. (23)). This is in
contrast to the stress-strain expressions of the previous material models for monotonic tension (Eqs. (10),
(17) and (20)), which are all explicit functions (i.e. σ(, x) can be directly computed when one has ). This
changes the likelihood function construction procedure for the linear elastic-nonlinear hardening case which
is discussed in subsection 3.4.

6

2.3. Markov chain Monte Carlo method (MCMC)
Once a posterior has been constructed, it needs to be analysed to determine the unknown parameters.
Only in a few straightforward cases, the analysis can be performed analytically. A numerical approach
needs to be employed in the majority of the cases, e.g. when the model is nonlinear or when the model is
an implicit function. If the posterior is continuous one may use approximation approaches (e.g. a Laplace
approximation [6]) that require the derivatives of the posterior. If the posterior is C0 -continuous however,
numerical approaches which do not require any derivatives need to be employed. The posteriors for for the
aforementioned elastoplastic models are all C0 -continuous due to the abrupt transition between the purely
elastic part and the elastoplastic part of the response (see e.g. Eq. (22)).
Markov chain Monte Carlo (MCMC) techniques are frequently employed numerical approaches to investigate posteriors, that do not require any derivatives [41–43]. MCMC approaches are based on drawing
samples from the posterior to approximate the parameters’ statistical properties (e.g. the mean value and
variance). Below, the fundamental concepts of the Monte Carlo method are discussed as well as the adaptive
Metropolis algorithm to perform the sampling.
2.3.1. Monte Carlo method
The main purpose of the Monte Carlo method is to approximate integrals of the following form:
Z
I=
g(x)π(x)dx,

(24)

Rn

where π is the PDF of interest (in our case the posterior) and g : Rn → Rk is an integrable function over
Rn . This integral can be approximated using the following quadrature:
Î =

N
1 X
g(xi ),
N i=1

(25)

N

where {xi }i is the set of samples drawn from the PDF of interest (π) and the hat on Î represents the
numerically approximated equivalent of I. The drawing of samples from π implies that most of the samples
are in the domain in which numerical evaluations of π do not equal zero. By the law of large numbers (LLN)
[44], Î converges as follows:
N
1 X
g(xi ) = I.
N →+∞ N
i=1

lim

(26)

d g ) is [45]:
The numerical approximation of the components of the covariance matrix for g(x) (cov
(d
covg )jm =

N


1 X
gj (xi ) − Ij (gm (xi ) − Im , j = 1, 2, · · · , k, m = 1, 2, · · · , k.
N i=1

(27)

The mean of the posterior (µpost ) can be computed by substituting g(x) = x and π = πpost in Eq. (24),
which yields:
Z
µpost =

x πpost (x)dx =
Rn

N
1 X
xi .
N →+∞ N
i=1

lim

(28)

Furthermore, the components of the posterior’s covariance matrix are approximated as follows:

(d
covpost )jm =

N


1 X
(xi )j − (µpost )j (xi )m − (µpost )m , j = 1, 2, · · · , k, m = 1, 2, · · · , k.
N i=1

7

(29)

If we again assume that a large number of samples is taken (i.e. N is large), the MAP estimate can
furthermore be approximated as [44]:
\ = argmax π(xi ).
MAP

(30)

xi;i=1,...,N

The essential part of a Monte Carlo procedure is the drawing of admissible samples (xi ). Below, the
standard and adaptive Metropolis algorithms are discussed. The adaptive one is the algorithm used in
section 5.
2.3.2. The standard and the adaptive Metropolis algorithm
The standard Metropolis-Hastings approach is a frequently employed MCMC algorithm [44]. The basic
idea of the Metropolis-Hastings algorithm is to explore the PDF of interest by making a random walk
across the parameter space x. Considering sample xi and an evaluation of the PDF π(xi ), a new sample
xp is proposed by drawing from a proposal distribution (q in Algorithm 1). If the PDF evaluated at the
proposed sample (π(xp )) is larger than at the current sample (π(xi )), the proposed sample is always accepted.
However, if the PDF at the proposed sample is smaller than at the current sample, the proposed sample may
be accepted based on the ratio of the PDF evaluated at the proposed and current sample (r in Algorithm
1). The ratio is compared to a random number generated from a uniform distribution. If the ratio is greater
than the random number, the proposed sample is accepted. If the ratio is smaller than the random number,
the proposed sample is rejected, and we remain at the current sample. In the case the proposed sample
was accepted, the proposed sample becomes the current sample. The algorithm repeats until convergence is
achieved.
Critically, as the number of samples increases, it can be shown that the distribution of values approximates
the target distribution. For more information, see [45]. The standard Metropolis-Hastings algorithm is
presented in Algorithm 1 in more detail.
Algorithm 1 The standard Metropolis-Hastings algorithm
select the initial sample x0 ∈ Rn and γ
for i = 0, 1, 2, ..., N − 1 do
draw xp ∈ Rn from the proposal distribution
q(xp|xi ) in Eq. (33)

π(xp )q(xi |xp )
4:
calculate the ratio r(xi , xp ) = min 1, π(xi )q(xp |xi )
1:
2:
3:

5:
6:
7:
8:
9:
10:
11:

. π(·) denotes the target distribution (i.e. posterior).
draw u ∈ [0, 1] from uniform probability density
if r(xi , xp ) ≥ u then
xi+1 = xp
else
xi+1 = xi
end if
end for
In the case that the transition kernel is symmetric (as in this contribution), the following relation holds:
q(xi |xp ) = q(xp |xi ).

(31)

Consequently, step 4 in the algorithm simplifies to:
 π(x ) 
p
r(xi , xp ) = min 1,
.
π(xi )

(32)

The stability and convergence of the algorithm can be checked by tracing the generated samples and
analysing their characteristics. The evolution of the mean value and the standard deviation can for instance

8

be checked for convergence [28]. The distribution approximated by MCMC converges to the target distribution if the samples show a stationary statistical behaviour, i.e. if the mean and standard deviation remain
constant [28].
The efficiency of the algorithm is influenced by the initial sample (x0 ) and the proposal distribution (q)
[8]. The most common proposal distribution for the Metropolis-Hastings algorithm (as employed here) is of
the following Gaussian form:


1
2
q(xi |xp ) = q(xp |xi ) ∝ exp − 2 kxi − xp k ,
2γ

(33)

where γ is the parameter that determines the width of the proposal distribution and must be tuned to obtain
√
[46], where n is the number of
an efficient and converging algorithm. An efficient starting value is γ = 2.38
n
unknown parameters and hence, the dimension of the posterior.
To overcome the tuning of γ, Haario et al. [47] introduced the Adaptive Proposal (AP). The AP method
updates the width of the proposal distribution, based on the existing knowledge of the posterior. The existing
knowledge is based on the previous samples. For sample i + 1, the update employs the following formulation:
q(xp |xi ) ∼ N (xi , γ 2 Ri ),

(34)

where N (xi , γ 2 Ri ) denotes a normal distribution with mean xi and covariance matrix γ 2 Ri , of size n × n.
To establish Ri , all i previous samples are first stored in matrix K of size i × n. Ri , is then computed as:
Ri =

1 eT e
K K,
i−1

(35)

e = K − Kmean and Kmean reads:
where K

Kmean



kmean
kmean 


= . 
 .. 
kmean

,

(36)

i×n

and kmean is a row matrix of length n which is determined as follows:
 i

i
i
P
P
1 P
(K)j1
(K)j2 · · ·
(K)jn .
kmean =
i j=1
j=1
j=1

(37)

The following relation is used for N (xi , γ 2 Ri ) in this contribution:
N (xi , γ 2 Ri ) = xi + √

γ
e T N (0, Ii ),
K
i−1

(38)

where Ii is the identity matrix of size i × i and N (0, Ii ) is the i-dimensional normal distribution.
Note that it is computationally inefficient to update the proposal distribution after each new sample is
generated. In the numerical examples in this study therefore, updating takes place once per 1000 sample generations. Algorithm 2 shows the Metropolis-Hastings algorithm with the symmetric AP proposal (Eq. (38))
that is employed here.
3. Bayesian inference for uniaxil tensile tests with noise in the stress measurements
The current section focuses on the formulation of Bayesian identification approaches for the material
models of subsection 2.2, when only the stress measurements are polluted with statistical errors. The
measured strains are thus considered to be exact. In section 4, Bayesian updating is considered when also
the measured strains are also polluted with a statistical noise.

9

Algorithm 2 The Metropolis-Hastings algorithm with symmetric AP proposal
√
select the initial sample x0 ∈ Rn and set γ = 2.38
n
2: for i = 0, 1, 2, ..., N − 1 do
3:
draw xp ∈ Rn from the proposal distribution
 q(xp |xi ) in Eq. (38)


1:

4:
5:
6:
7:
8:
9:

10:
11:
12:
13:

calculate the ratio r(xi , xp ) = min 1,

π(xp )
π(xi )

. π(·) denotes the target distribution (i.e. posterior).
draw u ∈ [0, 1] from uniform probability density
if r(xi , xp ) ≥ u then
xi+1 = xp
else
xi+1 = xi
end if
per 1000 samples
e
update matrix K
end for

3.1. Noise distribution
To determine the noise distribution and its parameters, two sets of ‘calibration experiments’ can be
performed. First, a test is performed without any specimen. The stress-strain measurements of this test are
shown in Fig. 1(a). It shows that the PDF of the noise in the ‘stress measurements’ is a normal distribution
with a zero mean and a standard deviation of Snoise .
Second, the evolution of the noise distribution (including its parameters) must be determined. To this
purpose, a tensile test is performed on a calibration specimen (of which the Young’s modulus is known).
The assumed results are presented in Fig. 1(b). The mean stress value varies linearly with the strain. The
standard deviation Snoise however remains the same.
As the ‘calibration measurements’ indicate that an additive noise model can be used and the stresses are
polluted by a normal noise distribution with standard deviation Snoise , BI is first employed to identify the
Young’s modulus of the linear elastic model.
3.2. Linear elastic
The only unknown material parameter in the linear elastic model is the Young’s modulus (E). Based on
section 2, the additive noise model for a single stress measurement can be written as follows:
σ m = E + Ω,

(39)

m

where σ is the measured stress and Ω is the random variable representing the noise in the stress measurement. Based on the artificially generated calibration results, the noise distribution is normal and can be
written as:
πnoise (ω) = √


1
ω2 
exp −
.
2
2Snoise
2πSnoise

(40)

Using Eq. (8), the likelihood function is expanded as:
π(σ m |E) = πnoise (σ m − E) = √

 (σ m − E)2 
1
exp −
,
2
2Snoise
2πSnoise

(41)

Substitution of Eq. (41) in Eq. (9) then yields the following expression for the posterior:
 (σ m − E)2 
π(E|σ m ) ∝ π(E)exp −
.
2
2Snoise
10

(42)

σ

σ

πnoise (ω)




(a) Fitted noise distribution

(b) The shifted noise distribution

Figure 1: Schematic of the stress-strain measurements (red circles) of the ‘calibration experiments’, including
the noise distributions (dashed). The theoretical stress-strain relation (which is exact for the calibration
experiments) is presented as the bold straight line on the diagram on the right.
If we use a prior in the form of a modified normal distribution as follows:


(
2
exp − (E−E)
if E ≥ 0
2
2SE
π(E) ∝
,
0
otherwise

(43)

the posterior distribution reads:
m

π(E|σ ) ∝

 h
(
2
exp − (E−E)
+
2S 2
E

(σ m −E)2
2
2Snoise

i

0

if E ≥ 0

,

(44)

otherwise

where E and SE are the mean and standard deviation of the prior distribution, respectively. Note that
the Young’s modulus cannot be negative for an actual material, which is taken into account in the prior
distribution (Eq. (43)).
If we now consider the posterior distribution of the previous measurement to be the prior distribution of
the current measurement, the posterior for all k measurements can be expressed as:
m



π(E|σ ) ∝ exp −

h

k
P

(E−E)2
2
2SE

+

(σim −Ei )2 i

i=1

2
2Snoise

,

E≥0

where π(E|σ m ) = π(E|σ1m , σ2m , · · · , σkm ). Eq. (45) can now be written in the following form:


(E−µ
)2
π(E|σ m ) ∝ exp − 2S 2post , E ≥ 0

(45)

(46)

post

where µpost and Spost are the mean and standard deviation of the posterior distribution, which is again a
normal distribution (with the condition E ≥ 0). Both can be expressed as:
2
2
Snoise
E+SE

µpost =

k
P
i=1
k
P

2
2
Snoise
+SE

i σim

i=1

,
2i

Spost =

s

2
2
Snoise
SE
.
k
P
2
2
Snoise
+SE
2i

(47)

i=1

Hence, it is possible to analytically examine the posterior distribution for the linear elastic model, in case
the noise model is additive and the noise distribution as well as the prior distribution are (modified) normal
11

distributions. For the other cases below, MCMC approaches are required and the adaptive MCMC approach
used for those is verified based on the analytical expressions in Eq. (47) for the linear elastic model.
3.3. Linear elastic-perfectly plastic
The parameters to be identified for the linear elastic-perfectly plastic model are the Young’s modulus

T
and the initial yield stress which are stored in parameter vector x = E σy0 . Since the experimental
equipment and experimental condition both remain the same (i.e. the measured stresses are still polluted
by noise stemming from the same normal distribution and the measured strains are still exact), the same
additive noise model applies:
σ m = σ(, x) + Ω,

(48)

where σ(, x) is referred to as the theoretical stress and takes the form of the stress-strain relation for
monotonic tension in Eq. (16). Using Eq. (40) for the noise distribution, the likelihood function for a single
stress measurement reads:
π(σ m |x) = πnoise (σ m − σ(, x)) = √

 (σ m − σ(, x))2 
1
,
exp −
2
2Snoise
2πSnoise

(49)

or:

1
π(σ |x) = √
exp
2πSnoise
m

σ

m



− E 1 − h  −

−

σy0
E




− σy0 h  −

σy0
E

 2
!
.

2
2Snoise

(50)

Taking the physical constraints into account that the Young’s modulus and the initial yield stress must
be nonnegative, the following prior distribution is selected:


(
(x−x)T Γ−1
x (x−x)
exp −
if E ≥ 0 and σy0 ≥ 0
2
π(x) ∝
,
(51)
0
otherwise
where x is the mean value vector of the prior distribution and Γx is the covariance matrix of the prior.
Substitution of Eq. (50) and Eq. (51) in the reduced variant of Bayes’ formula of Eq. (9), yields the following
posterior distribution for k measurements:
k
P

(x − x)T Γ−1
x (x − x) i=1
π(x|σ ) ∝ exp −
+
2


m




m
σi − Ei 1 − h i −

σy0
E




− σy0 h i −

σy0
E

 2
!
, (52)

2
2Snoise

where the probability of obtaining a negative Young’s modulus and initial yield stress is zero thanks to the
selected prior distribution.
It is important to realise that the posterior distribution in this case is not of a form that allows an
analytical evaluation due to the presence of the Heaviside function. Hence, the adaptive MCMC approach
is employed to analyse the posterior for the linear elastic-perfectly plastic model in section 5.
3.4. Linear elastic-linear hardening

T
The parameter vector for the linear elastic-linear hardening model reads x = E σy0 H . Assuming
the same experimental equipment and condition (and hence, the same noise model and noise distribution),
the likelihood function for a single measurement reads:

σ
m

π(σ |x) ∝ exp

m



− E 1 − h  −

σy0
E




− σy0 +

−

2
2Snoise

12

HE
H+E



−

σy0
E

 
h −

σy0
E

 2
!
.

(53)

In addition to the physical constraints for the Young’s modulus and the initial yield stress, the plastic
modulus (H) must also be nonnegative. The following prior distribution is therefore selected:


(
(x−x)T Γ−1
x (x−x)
if E ≥ 0 and σy0 ≥ 0 and H ≥ 0
exp −
2
.
(54)
π(x) ∝
0
otherwise
Using Bayes’ formula, the posterior distribution for k observations reads:

m

π(x|σ ) ∝ exp

(x − x)T Γ−1
x (x − x)
+
2



k
P
m
σi − Ei 1 − h i −


−

i=1

σy0
E




− σy0 +

HE
H+E


i −

σy0
E

 
h i −

σy0
E

2
2Snoise

2
!
. (55)

3.5. Linear elastic-nonlinear hardening

T
The parameter vector for the linear elastic-nonlinear hardening model is x = E σy0 H n . Considering no change of experimental equipment (and hence, the same noise model and noise distribution), the
expression for the measured stress reads again:
σ m = σ(, x) + Ω,
which, together with Eq. (23), results in the following expression for the measured stress:

 



σy0 
σ(, x) n 
σy0 
m
σ = E 1 − h  −
+ σy0 + H  −
h −
+ Ω.
E
E
E

(56)

(57)

It is important to note that the expression for the measured stress is now a function of the theoretical stress
(σ(, x)). This is in contrast to the expressions of the measured stresses of the other material models. Consequently, the construction of the likelihood function changes. First, one needs to determine the probability
that the measured stress occurs, for a given theoretical stress and a set of given material parameters:
π(σ m |x, σ(, x)) = πnoise (σ m − σ(, x)).

(58)

The likelihood function required for Bayes’ theorem however, expresses the likelihood that measured stresses
occur for given parameters. This can be obtained by integrating Eq. (58), together with the theoretical stress
likelihood function, over the theoretical stress [48]:
Z +∞
π(σ m |x) =
π(σ m |x, σ(, x))π(σ(, x)|x)dσ.
(59)
0

In Eq. (59), the theoretical stress is integrated from 0 to infinity, because it cannot be smaller than zero
in monotonic tension (otherwise compression would occur). As the theoretical stress can in principle be
infinitely large, no true upper bound is included.
In the case that the parameters are given (as in MCMC approaches in which each sample is a realisation
of the parameters), the theoretical stress likelihood function can be expressed as follows:
π(σ(, x)|x) = δ(σ(, x) − w),

(60)

where δ is the Dirac delta function and w is used here to represent the right hand side of Eq. (23). Note that
in the case that all the parameters are given, the only possible value for the theoretical stress is the same

13

expression of the theoretical stress, but only as a function of the parameters. In Eq. (60) this is imposed by
the Dirac delta function. Substitution of Eq. (58) and Eq. (60) in Eq. (59) yields:
Z +∞
πnoise (σ m − σ(, x))δ(σ(, x) − w)dσ.
(61)
π(σ m |x) =
0

The integral in Eq. (61) can be computed using [49]:
 n
Z b
 P f (ci )
0
f (z)δ(g(z))dz = i=1 |g (ci )|

a
0

if a < ci < b

,

(62)

otherwise

i)
6= 0. Combining Eqs. (61) and (62) yields the final
where ci are the i roots of g(z), for which holds dg(c
dz
likelihood function for a single measurement as:



σ
(σ m −E)2

if  ≤ Ey0
exp
−

2

2Snoise


π(σ m |x) ∝ exp − (σm −c)2
,
(63)

2S 2
σy0

noise

if

>
c n−1
Hn
E
|1+ E (− E ) |

where c is the value of the theoretical stress when σ(, x) − w = 0 is solved for σ(, x). Choosing the prior
distribution in the form of a modified normal distribution as:


(
(x−x)T Γ−1
x (x−x)
if E ≥ 0 and σy0 ≥ 0 and H ≥ 0 and n ≥ 0
exp −
2
.
(64)
π(x) ∝
0
otherwise
The final form of the posterior distribution reads:
k
 (x − x)T Γ−1 (x − x)  Y
x
π(σim |x),
π(x|σ m ) ∝ exp −
2
i=1

(65)

where π(σim |x) must be calculated according to Eq. (63) for each measurement. To investigate the posterior,
first σ(, x) − w = 0 is numerically solved for σ(, x) and then the adaptive MCMC is employed to find the
statistical characteristics of the posterior.
4. Bayesian inference for uniaxial tensile tests with noise in both stress and strain measurements
All the cases studied in section 3 consider an error in the stress, but not in the strain. Both the measured
stresses and the measured strains may be polluted by statistical errors however. Both statistical errors are
not related to each other, as the devices to measure the forces and strains are independent of each other
(e.g. when a load cell and digital image correlation are used). This section deals with the same four material
models as in the previous sections but considers that the stresses as well as the strains are influenced by
noise.
4.1. Noise distribution
The first calibration experiment is again performed without the use of a specimen. The stress-strain
measurements of this test are shown in Fig. 2(a). The measurements in Fig. 2(a) indicate a normal distribution with a zero mean and a diagonal covariance Γnoise . The noise model and the evolution of the noise
distribution are again investigated by performing a tensile test on a calibration specimen with an exactly
known Young’s modulus. The schematic results are shown in Fig. 2(b). It is clearly visible that the mean
stress and strain values follow a linear, theoretical stress-strain relationship and hence, the covariance does
not change. The calibration results show that the additive noise model can also be used in this section.
14

σ

σ

πnoise (ω)




(a) Fitted noise distribution

(b) The shifted noise distribution

Figure 2: Schematic of the stress-strain measurements (red circles) of the ‘calibration experiments’ for the
case with uncertainty in both stress and strain, including some isolines of the noise distributions (dashed).
The theoretical stress-strain relation (which is exact for the calibration experiments) is presented as the bold
straight line on the diagram on the right.
4.2. Linear elastic
The additive noise model when both the stresses and strains are contaminated by stochastic noise, can
be expressed as follows for the linear elastic model:
(
σ m = E + Ωσ
,
(66)
m =  + Ω
where σ m is the measured stress, m is the measured strain, Ωσ is the stochastic error of the stress measurement and Ω is the stochastic error of the strain measurement. Because the information from both the
measured stress and the measured strain is used here, Bayes’ formula for multiple variables must be employed
[50]:
π(E|σ m , m ) =

π(E)π(m |E)π(σ m |E, m )
.
π(m )π(σ m |m )

(67)

Since the measured strain and the Young’s modulus are statistically independent (measuring a specific strain
does not interfere with the probability that a certain Young’s modulus occurs, i.e. π(m |E) = π(m )), Eq. (67)
can be written as follows:
π(E|σ m , m ) =

π(E)π(σ m |E, m )
,
π(σ m |m )

(68)

or, when one is only interested in relative probabilities, as:
π(E|σ m , m ) ∝ π(E)π(σ m |E, m ).
m

(69)
m

The issue with Eq. (69) is that the the likelihood function, π(σ |E,  ), must again be determined by
integration (over  here [48]), because π(σ m |E, ) can be determined directly, but π(σ m |E, m ) not. To this
end, we write:
Z a
π(σ m |E, m ) =
π(σ m |E, )π(|m )d,
(70)
0

15

where a is defined by the physical upper bound of the tensile tester (i.e. the ratio of the original length of
the specimen and the maximum distance that the clamps can move). Based on Eq. (66) one can express the
conditional probabilities π(σ m |E, ) and π(|m ) as follows:
(
π(σ m |E, ) = πΩσ (σ m − E)
,
(71)
π(|m ) = πΩ (m − )
where πΩσ (ωσ ) and πΩ (ω ) are the noise distributions of the errors in the stress measurements and the strain
measurements, respectively. Based on the calibration results in the previous subsection, the noise distribution
is a two-dimensional normal distribution with no correlation between the noise in the stress measurements
and the noise in the strain measurements (diagonal covariance matrix) and zero mean value as follows:


ω T Γ−1
1
noise ω
p
,
(72)
exp −
πnoise (ω) =
2
2π |Γnoise |

T
where ω = ωσ ω is the noise vector and Γnoise is the covariance matrix of the noise distribution in the
following form:

 2
Sσ 0
,
(73)
0 S2
Combining Eqs. (71) and (72), the following two conditional probabilities are required to construct the
likelihood function (Eq. (70)):



m
2
π(σ m |E, ) = √ 1 exp − (σ −E)
2
2S
2πSσ 
σ
,
(74)
m
2
π(|m ) = √ 1 exp − ( −)
2
2S
2πS




Substitution of Eq. (74) in Eq. (70) then results in the following expression for the likelihood function:
 h m

Z a
1
(m − )2 i
(σ − E)2
π(σ m |E, m ) =
+
exp −
d.
(75)
2πSσ S 0
2Sσ2
2S2
The result of this integral can be expressed analytically as:
√
 p − p2   a − p 
 −p 
p3
2
1
1
1
m
m
exp −
π(σ |E,  ) = √
erf √
− erf √
,
2p3
2p3
2p3
2 2πSσ S

(76)

where erf(·) is the error function [49] and p1 , p2 and p3 are formulated as follows:
p1 =

Eσ m S2 + Sσ2 m
Sσ2 (m )2 + S2 (σ m )2
(Sσ S )2
,
p
=
,
p
=
.
2
3
E 2 S2 + Sσ2
E 2 S2 + Sσ2
E 2 S2 + Sσ2

(77)

Consequently, choosing the same prior as in Eq. (43), the posterior distribution for k measurements reads:

π(E|σ m , m ) ∝ exp

 (E − E)2 
2
2SE

√
√

p3

2 2πSσ S

k Y
k

 (p − p2 )   a − p 
 −p 
2i
1i
1i
1i
exp −
erf √
− erf √
, (78)
2p
2p
2p
3
3
3
i=1

m m
m m
where π(E|σ m , m ) = π(E|(σ1m , m
1 ), (σ2 , 2 ), · · · , (σk , k )).

16

4.3. Linear elastic-perfectly plastic
In case the experimental devices remain the same, but the linear elastic-perfectly plastic model is used,
the additive noise model reads:
(
σ m = σ(, x) + Ωσ
(79)
m =  + Ω 

T
where σ(, x) is given by Eq. (17) and x = E σy0 . For the same noise distribution as in Eq. (72), the
conditional probabilities required to construct the likelihood function (Eq. (70)) are:

 



2
σy0
σy0
!

m

−σy0 h − E
σ −E 1−h − E


1
m
π(σ |E, ) = √2πS exp −
2
2Sσ
.
(80)
σ




m
2

π(|m ) = √ 1 exp − ( −)
2S 2
2πS




Combining Eqs. (70) and (80) the likelihood function can be written as:

 p − p2 h  σy0 − p 
 −p i
1
√
1
2
1
1
√
π(σ |x,  ) =
p3 exp −
erf E√
− erf √
+
2p
2p
2p
2 2πS Sσ
3
3
3
 (σ m − σ )2 h  m − σy0 
 m − a i
y0
E
√
√
erf
S exp −
−
erf
. (81)
2Sσ2
2S
2Sσ
m

m

where p1 , p2 and p3 again given by Eq. (77). Finally, selecting the prior distribution as in Eq. (51), the
posterior distribution for k measurements reads:
m

m

π(x|σ ,  ) ∝ exp

k
 (x − x)T Γ−1 (x − x)  Y
x

2

π(σim |x, m
i ).

(82)

i=1

4.4. Linear elastic-linear hardening

T
For the linear elastic-linear hardening model, the unknown parameters are x = E σy0 H . Assuming
the same experimental equipment as in subsections 4.2 and 4.3 and employing Eqs. (20) and (70), the
likelihood function can be expressed as:

 p − p2 h  σy0 − p 
 −p i
1
√
1
2
1
1
π(σ m |x, m ) = √
p3 exp −
erf E√
− erf √
+
2p3
2p3
2p3
2 2πSσ S
√
√
 β β − β 2 h  β1 a − √β2 
 σy0 β1 −
1
E
β
1 3
2
1
√ exp −
√
√
erf
− erf
2β1
β1
2
2

√β2
β1

i

, (83)

where p1 , p2 and p3 are again given by Eq. (77) and β1 , β2 and β3 are determined as follows:

β1 =

β3 =

HE
H+E

Sσ2

2
1
+ 2 , β2 =
S

HE
(σ m − σy0 ) H+E
+



HE
H+E

2

σy0
E

Sσ2

HE
(σ m − σy0 )2 + 2(σ m − σy0 ) H+E

Sσ2

17

σy0
E

+



HE
H+E

2 

σy0
E

+

m
,
S2

2
m 2

+

( )
.
S2

(84)

Selecting the same posterior distribution as in Eq. (54), the posterior distribution is again of the following
form:
π(x|σ m , m ) ∝ exp

k
 (x − x)T Γ−1 (x − x)  Y
x

2

π(σim |x, m
i ).

(85)

i=1

m m
where π(σim |x, m
i ) is calculated using Eq. (83) for each pair of measurements (σi , i ).

4.5. Linear elastic-nonlinear hardening

T
Like in subsection 3.5 the parameter vector for the current case reads x = E σy0 H n . As the
experimental equipment remains the same, the additive noise model is still valid and hence, the following
expressions for the measured stresses and strains are employed:
(
σ m = σ(x, ) + Ωσ
.
(86)
m =  + Ω 
In contrast to the previous models in this subsection, the theoretical stress (σ(x, ) in the first equation of
Eq. (86)) is an implicit function (Eq. (23)) and as a result the construction of the likelihood function differs.
First, Eq. (86) is rewritten using Eq. (23):



 

n  

σy0
σ(,x)
σ m = E 1 − h  − σy0
+
σ
+
H

−
h

−
+ Ωσ
y0
E
E
E
,
(87)
 m
 =  + Ω
from which it is clear that the measured stress is not only a function of the material parameters (x) and the
theoretical strain () but also of the theoretical stress (σ(x, )).
The following conditional probabilities for the measured stress and measured strain are established:
(
π(σ m |x, , σ(x, )) = πΩσ (σ m − σ(x, ))
.
(88)
π(|m ) = πΩ (m − )
The required likelihood function then reads:
Z a Z +∞
m
m
π(σ |x,  ) =
π(σ m |x, , σ(x, ))π(σ(x, )|x, )π(|m )dσd,
0
m

(89)

0

where π(σ m |x, , σ(x, )) and π(| ) are given by Eq. (88). If the material parameters and the theoretical
strain are given, the conditional probability π(σ(x, )|x, ) reads:
π(σ(x, )|x, ) = δ(σ(x, ) − w),

(90)

where w is again used to abbreviate the right hand side of Eq. (23). Consequently, the resulting likelihood
function reads:
1
π(σ |x,  ) =
2πSσ S
m

m

−

a

+∞

(m − )2
2S2
0
0
 



σ
σ m − E 1 − h  − Ey0 − σy0 + H  −

Z

Z



δ(σ(x, ) − w)exp −

2Sσ2

σ(x,) n
E



h −

σy0 
E

2

dσd. (91)

This can furthermore be expressed as:
π(σ m |x, m ) =

1
(I1 + I2 ),
2πSσ S
18

(92)

where
σy0
E

Z

Z

+∞

I1 =
0

0


2


m
σ
−
E
m
2
( − )
δ(σ(x, ) − w)exp −
−
dσd,
2S2
2Sσ2

(93)

and

Z
I2 =

a
σy0
E

Z

+∞

0




σ(x,) n

m
σ
−
σ
+
H

−
h −
m
2
y0
E
( − )
δ(σ(x, ) − w)exp −
−
2
2S
2Sσ2

σy0 
E

2

dσd. (94)

Using the Dirac delta function properties and Eq. (62), I1 equals to:
√
 −p i
 p − p2 h  σy0 − p 
πp3
1
1
2
1
− erf √
,
I1 = √ exp −
erf E√
2p3
2p3
2p3
2

(95)

where p1 , p2 and p3 are given by Eq. (77). To calculate I2 , one first needs to find the solution for g(z)
in Eq. (62), for which it must be substituted by w. Again using the Dirac delta function properties and
integrating Eq. (94) over the theoretical stress, I2 reads:
Z
I2 =

a
σy0
E



exp −

m

2

( − )
−
2S2

σ m − σy0 + H  −

  2 
w n
E

2Sσ2

1+

1
−

Hn
E (

w n−1
E)

d,

(96)

which can be calculated by numerical approaches (e.g. Simpson’s rule). Finally for the same prior distribution
as in Eq. (64), the posterior distribution for k observations becomes:
π(x|σ m , m ) ∝ exp

k
 (x − x)T Γ−1 (x − x)  Y
x

2

π(σim |x, m
i ).

(97)

i=1

Note that in practice once a sample is drawn by the adaptive MCMC approach, the integral in Eq. (97) and
subsequently the likelihood function can be calculated.
5. Examples
All formulations derived in the previous two sections are investigated below. The effect of the prior
distribution on the posterior distribution is studied, as well as BI’s ability to recover a material parameter
distribution when they are taken from a specific distribution. Also, BI’s ability to recover correlations
between different material parameters is exposed.
5.1. Bayesian inference with noise in stress
This subsection presents several examples for the BI formulations of section 3, in which only a statistical
noise in the stress measurements is considered.
5.1.1. Linear elastic (LE)
Identification of the Young’s modulus. In the first example a specimen with a Young’s modulus of 210 GPa
is considered, which is to be identified. ‘Calibration experiments’ were performed and the noise in the stress
follows the normal distribution of Eq. (40) with Snoise = 0.01 GPa. For only one stress measurement of
σ m = 0.1576 GPa with corresponding strain  = 7.25 × 10−4 , the posterior distribution is calculated using
Eq. (46). Selecting the prior distribution as in Eq. (43) with mean E = 150 GPa and a relatively large
standard deviation SE = 50 GPa, the posterior reads:


(E−µ
)2
π(E|σ m ) ∝ exp − 2S 2post , E ≥ 0,
(98)
post

19

where µpost = 212.6486 GPa and Spost = 13.2964 GPa.
Fig. 3 shows this posterior distribution, as well as the prior distribution and the value predicted by the
least squares method when one measurement is made and also when five measurements are made. Fig. 4
presents the linear elastic responses for the case when one measurement is made and for the case when ten
measurements are made. The figure also shows the stress-strain responses made with Young’s moduli taken
within the 95% credible region (i.e. the region that contains 95% of the posterior) of the posterior.
Two points can be observed based on Fig. 3. First, the strain at which a measurement is made has a
strong influence on the posterior. This can be observed when the posterior of Fig. 3(a) is compared to the
posterior of Fig. 3(b) when only the first measurement is incorporated (the distribution in red, denoted by
π(E|σ1m )). The latter distribution is significantly wider and its MAP point is relatively distant from the
specimen’s Young’s modulus. Hence, a measurement made at a comparatively large strain reduces the width
of the posterior distribution (i.e. reduces the uncertainty).
The second remark is that for an increasing number of measurements, the posterior becomes narrower
and the MAP point moves closer to the specimen’s Young’s modulus.
When the MAP point for a single measurement in Fig. 3(b) (µpost = 207.2821 GPa) is compared with
the result of the least squares method for the same measurement (Els = 210.2216 GPa), one can notice the
effect of the selected prior distribution. One interpretation of this is that the least squares method gives a
more accurate result than BI (although this depends the selected prior), as the result of the least squares
method is closer to the specimen’s Young’s modulus than the MAP point determined using BI. On the
other hand, the result determined using the least squares method is not the actual Young’s modulus of the
specimen (210 GPa), whereas the posterior distribution of BI does include this value. Furthermore, the MAP
point and mean value of BI, comes with an uncertainty in terms of the parameter value itself. This can be
considered an advantage when one wants to include this uncertainty, instead of including one deterministic
value.
The main point is that BI cannot be directly compared to the least squares method, because contrary to
the latter, it considers the Young’s modulus to be stochastic in nature. This means that BI accounts for the
fact that parameters originate from a distribution (for which the user is required to give an a-priori belief,
which is reflected in the prior). In the least squares method the Young’s modulus is a deterministic value,
which can only be determined using the current measurements. Hence, it cannot consider the effect that
other measurements could have been made (due to different noise realisations).

20

0.8

1
:(E)
:(Ej<1m )
:(Ej<1m ; <2m )
:(Ej<1m ; <2m ; <3m )
:(Ej<1m ; <2m ; <3m ; <4m )
:(Ej<1m ; <2m ; <3m ; <4m ; <5m )
Least squares method

0.7
0.8

0.6

0.6

0.4

:

:

0.5
:(Ej< m )

0.3

0.4

:(E)

0.2

Least squares method

0.2

0.1
0

0

100

200
300
E (GPa)

400

0

500

0

(a) One observation

100

200
300
E (GPa)

400

500

(b) Five successive observations

Figure 3: Linear elastic: The prior, the posterior and the value predicted by least squares method for one
measurement (a) and five measurements (b). The distributions are not normalised. The strain at which a
measurement is made has a considerable influence on the posterior. This can be observed when the posterior
of (a) (π(E|σ m ), red line) is compared to the posterior of (b) when only the first measurement is incorporated
(π(E|σ1m ), red line). An increase of the number of measurements leads to narrow posteriors with their MAP
estimates closer to the specimen’s Young’s modulus.

0.3

0.3
95% credible region curves
Measured value
Mean

0.25

0.25
0.2

< (GPa)

< (GPa)

0.2
0.15

0.15

0.1

0.1

0.05

0.05

0

95% credible region curves
Measured value
Mean

0

0.2

0.4

0.6
0

0.8

1

0

1.2
#10!3

(a) One observation

0

0.2

0.4

0.6
0

0.8

1

1.2
#10!3

(b) Ten observations

Figure 4: Linear elastic: The measurements, and the stress-strain curves created using the posterior for (a)
one measurement and (b) ten measurements. Increasing the number of measurements leads to a smaller
uncertainty, since the bandwidth of the responses associated with the 95% credible region is much smaller
in (b) than in (a).
The influence of the prior on the identified Young’s modulus. It is nevertheless interesting to study the effect
of the prior distribution on the MAP point (which is the same as the mean value for the normal posteriors in
this subsection). In Fig. 5 the MAP points are shown as a function of the mean and the standard deviation
21

of the prior. The MAP points are presented for different numbers of measurements. As can be seen, an
increase of the number of measurements leads to a flatter surface, which means that the influence of the
prior distribution decreases.
10 measurements
5 measurements
1 measurement
True value

7post (GPa)

1000
800
600
400
200
0
60
1000

40
500

20
SE (GPa)

0

0

E (GPa)

Figure 5: Linear elastic: The influence of the prior (i.e. the mean value and the standard deviation) on the
resulting MAP point for different numbers of measurements. Increasing the number of the measurements
leads to a flatter surface which indicates a decreasing influence of the prior distribution.
Recovering the heterogeneity of the Young’s modulus. A last important point to show using the linear elastic
model is BI’s ability (or inability for the current formulation) to capture the heterogeneity in the material
parameters. The question here is thus if BI is able to recover the distribution of the Young’s modulus when
multiple specimens are tested and their Young’s moduli are taken from a specific underlying distribution.
To this end, 25 specimens are considered of which the Young’s moduli are taken from a normal distribution
with a mean value of 210 GPa and a standard deviation of 10 GPa (blue curve in Fig. 6). For each specimen
ten measurements are made. The same noise model and noise distribution are applied.
The resulting posterior is represented by the red curve in Fig. 6, which is a (modified) normal distribution
with µpost = 215.3971 GPa and Spost = 0.8561 GPa. The posterior is substantially narrower than the
distribution of the specimens’ Young’s moduli and hence, using the BI formulations of this contribution,
the material heterogeneity cannot be captured. This entails that the width of the posterior distributions
(represented by Spost in this subsection) is only a measure of the uncertainty of the MAP points and the
mean value and not of the material heterogeneity.
To be able to recover the material heterogeneity one needs to consider both the inherent uncertainty in
the material parameters as well as the uncertainty in the measurements. In that case, one does not only
aim to determine the posterior of the material parameter(s), but also the posterior of the variance of the
assumed material parameter distribution. Our future work will focus on this.

22

1
Specimen elastic
modulus PDF
Posterior PDF

0.8

:

0.6

0.4

0.2

0
160

180

200
220
E (GPa)

240

260

Figure 6: Linear elastic: The distribution of the specimens’ Young’s moduli and the resulting posterior. The
PDFs are not normalised. The current formulation is clearly not able to recover the material heterogeneity.
To be able to recover the material heterogeneity one needs to consider both the inherent uncertainty in the
material parameters as well as that in the measurements.
5.1.2. Linear elastic-perfectly plastic (LE-PP)
Identification of the material parameters. In the first example of this subsection one linear elastic-perfectly
plastically behaving specimen is considered with Young’s modulus E = 210 GPa and initial yield stress
σy0 = 0.25 GPa. Twelve measurements are generated by employing the same noise distribution as in the
previous subsection. The prior distribution of Eq. (51) is furthermore selected with the following mean vector
and covariance matrix:




200
2500
0
x=
GPa, Γx =
GPa2 .
(99)
0.29
0
2.7778 × 10−4
Consequently, the posterior for the LE-PP model of subsection 2.2.2 is of the form as presented in Eq. (52),
which is investigated by the MCMC approach given in subsection 2.3.2. Running the chain for 104 samples
whilst burning the first 3000 samples (i.e. the first 3000 samples are not used to determine the mean, the
covariance matrix and the MAP estimate) yields:




29.807
4.1064 × 10−4
208.9859
b post =
b post =
µ
GPa, Γ
GPa2 ,
(100)
0.2578
4.1064 × 10−4 1.5067 × 10−5
and


\ = 208.4475 GPa,
MAP
0.2578

(101)

where the hat sign (b·) denotes the numerical approximation.
Fig. 7(a) shows the samples generated by the adaptive MCMC approach which are used to approximate
the posterior distribution. The domains presented in Fig. 7(b) show which of the measurements are included
in the elastic part and which fall within the plastic part. These discrete domains are a result of the C0 continuity of Eq. (17). In domain ‘a’ (in which no samples are generated by the adaptive MCMC approach),
all the measurements are considered to be in the plastic part. In domain ‘b’ on the other hand, the first
measurement (the one with the smallest strain) is considered to be in the elastic part and the others remain
in the plastic part. Continuing like this, in domain ‘c’ the second measurement is also considered to be in

23

the elastic part. Finally, in domain ‘m’ all measurements are considered to fall within the plastic domain.
Based on Fig. 7(b) the MAP point is clearly located in the domain in which the first six measurements are
considered to be in the elastic part and the remaining in the plastic part.
The 95% credible region is shown together with the posterior distribution in Fig. 8(a). The possible
stress-strain responses inside the credible region are presented in Fig. 8(b).
#10 -4

0.6

0.8

0.4

<y0 (GPa)

0.5

:

#10!3
1

0.6
0.4

..
..
.
...
..
..
..
..
..
..
..
..

0.3
0.2

0
0.3

0.1
0.25

0
180

200
0.2

150

6
5

Domain g

4
3

..
.
..
..
..

2

Domain c
Domain b
Domain a

250

200

220

8
7

Domain f

0.2

<y0 (GPa)

Domain m

1

240

0

E (GPa)

E (GPa)

(a) Samples generated by the adaptive MCMC approach (b) Samples generated by the adaptive MCMC approach
(top view) including different domains

Figure 7: Linear elastic-perfectly plastic: Two different views of the samples generated by the adaptive
MCMC approach to approximate the posterior. The colours represent the value of the posterior, which in
the left image is also shown along the z-axis. In Fig. 7(b) several domains are shown. For each of these
domains, the measurements considered to be in the elastic part and the measurements considered to be in
the plastic part remain the same.

24

#10 -4

0.3

95% credible region

0.3

8
7

0.25

0.28
6

4

0.24

< (GPa)

<y0 (GPa)

0.2
5

0.26

3

0.15
0.1

2

0.22

95% credible region curves
Measured value
Mean

0.05
1

0.2

160

180

200
220
E (GPa)

240

260

0

0

0

0.5

1

1.5
0

(a) The 95% credible region

2

2.5
#10!3

(b) The measurements and the stress-strain curves

Figure 8: Linear elastic-perfectly plastic: The 95% credible region and the posterior distribution (a), the
measurements and the stress-strain curves created using the posterior (b). The responses associated with
the 95% credible region in (b) are plotted using the points inside the 95% credible region marked by the red
ellipse in (a).
The posterior distribution seems be roughly of an an elliptical shape with the primary axes almost
along the E-axis and σy0 -axis. This entails that the correlation between the two material parameters is not
significant. One has to notice though, that the assumed prior is uncorrelated. In other words, the prior
covariance matrix (Γx ) is diagonal. It is therefore interesting to investigate the influence of the off-diagonal
term of the prior covariance matrix on the posterior covariance matrix. In Fig. 9 this influence is graphically
presented for the three terms of the posterior covariance matrix (note that both the prior covariance matrix
and the posterior covariance matrix are symmetric). It seems that an increase of (Γx )12 leads to some
b post )11 and some increasing trend for (Γ
b post )12 . However, it is difficult to assess
decreasing trend for (Γ
whether or not these trends can be considered as meaningful.
The influence of the prior on the correlation between the material parameters. The next example focuses
on the ability of the current formulation to capture a correlation between the Young’s modulus and the
initial yield stress when they are correlated. To this end, ten specimens are considered of which the material
parameters are governed by a normal distribution with the following mean vector and covariance matrix:




210
100
10−4
GPa2 .
(102)
xspc =
GPa, Γspc =
0.25
10−4 1.1111 × 10−4
For each specimen twelve measurements are made. Using the same prior as in the previous example (see
Eq. (99)) and the adaptive MCMC approach for 104 samples whilst burning the first 3000 samples yields:




5.5373
−8.396 × 10−4
211.1077
b
b post =
µ
GPa, Γpost =
GPa2 .
(103)
0.2519
−8.396 × 10−4 1.8174 × 10−6
The MAP estimate is given by:


\ = 210.5923 GPa.
MAP
0.2521

(104)

These results show that the correlation of the posterior is not same as that of the distribution of the
actual material. This corresponds closely with the observation that the formulations in this contribution are
25

b post )11
(a) Effect of (Γx )12 on (Γ

b post )12
(b) Effect of (Γx )12 on (Γ

b post )22
(c) Effect of (Γx )12 on (Γ

Figure 9: Linear elastic-perfectly plastic: Effect of the off-diagonal component of the prior covariance matrix
on the posterior’s covariance matrix. It seems that an increase of (Γx )12 leads to some decreasing trend of
b post )11 and some increasing trend of (Γ
b post )12 . However, it is difficult to assess whether or not a true trend
(Γ
is present in these results. The hat sign (b·) indicates that the values are computed numerically using the
adaptive MCMC approach.

26

b post )11
(a) Effect of (Γx )12 on ((Γ

b post )12
(b) Effect of (Γx )12 on ((Γ

b post )22
(c) Effect of (Γx )12 on (Γ

Figure 10: Linear elastic-perfect plastic: Effect of the off-diagonal components of the prior covariance matrix
on the posterior’s covariance matrix when the measurements are generated from specimens drawn from a
normal distribution given in Eq. (102). No real trends can be observed. The hat sign (b·) indicates that the
values are computed numerically using the adaptive MCMC approach.
not able to capture any of the intrinsic uncertainty of the material parameters. Fig. 10 shows the effect of the
off-diagonal component of the prior covariance matrix (Γx ) on the components of the posterior’s covariance
matrix (Γpost ). Again, no specific trends can be observed.
Application of the perfect plasticity formulation for a material that hardens nonlinearly. In the last example
of this subsection, we will regard how the formulation for the linear elastic-perfectly plastic model behaves
when it is used for a material that hardens nonlinearly. To this end, fifteen measurements are generated
using the linear elastic-nonlinear hardening model with E = 210 GPa, σy0 = 0.25 GPa, H = 2 GPa and
n = 0.5 and the same noise model and noise distribution as in the previous examples. Employing the
prior distribution in the form of Eq. (51) with the mean and covariance matrix of Eq. (99) and running the
adaptive MCMC approach for 104 samples whilst burning the first 3000 samples leads to:




31.2818
−3.0214 × 10−3
206.4357
b
b post =
µ
GPa, Γpost =
GPa2 ,
(105)
0.2941
−3.0214 × 10−3 1.1279 × 10−5
and


\ = 204.8997 GPa.
MAP
0.2944

(106)

The stress-strain curves associated with the 95% credible region in Fig. 11 show that six measurements
belong to the elastic part of the response for the mean.

27

0.35
0.3

< (GPa)

0.25
0.2
0.15
0.1
95% credible region curves
Measured value
Mean
True value curve

0.05
0

0

0.5

1

1.5

2

0

2.5
#10!3

Figure 11: Linear elastic-perfectly plastic: The measurements, the stress-strain curves using µpost including
the 95% credible region and the linear elastic-nonlinear hardening curve used to create the measurements.
5.1.3. Linear elastic-linear hardening (LE-LH)
Identification of the material parameters. This subsection deals with the Bayesian formulation for the linear
elastic-linear hardening model. A specimen with Young’s modulus E = 210 GPa, initial yield stress σy0 =
0.25 GPa and plastic modulus H = 50 GPa is regarded. Twelve measurements are created by employing the
same noise distribution as in the previous subsection. The prior distribution is given by Eq. (54) with the
following properties:




200
2500
0
0
2.7778 × 10−4
0  GPa2 .
x = 0.29 GPa, Γx =  0
(107)
60
0
0
100
The adaptive MCMC algorithm for 104 samples whilst burning the first 3000 samples gives:

b post
µ




36.5642
207.4586
b



0.2533
−1.2746
× 10−2
GPa, Γpost =
=
55.9187
−3.7886

−1.2746 × 10−2
4.0359 × 10−5
−2.6218 × 10−2


−3.7886
−2.6218 × 10−2  GPa2 ,
66.8214

(108)

and


206.9528
\ =  0.2548  GPa.
MAP
55.2838

(109)

Fig. 12 shows the generated samples by the adaptive MCMC approach in the E − σy0 − H space, including
the projections on the E − σy0 , E − H and σy0 − H planes.
The 95% credible region is presented in Fig. 13(a) and the possible stress-strain responses associated with
it are shown in Fig. 13(b).

28

#10 -3

1.5

90

H (GPa)

80
70

1

60
50
40
30

0.5

20
0.3
240

0.25

220
200

<y0 (GPa)

0

180

0.2

E (GPa)

Figure 12: Linear elastic-linear hardening: Samples generated by the adaptive MCMC approach to approximate the posterior distribution and its projection on the three planes.

0.35

80

0.25

60

0.2

< (GPa)

H (GPa)

0.3

40

0.15
0.1

20
0.3
0.25
<y0 (GPa)

0.2

180

220
200
E (GPa)

95% credible region curves
Measured value
Mean

0.05

240

0

0

0.5

1

1.5
0

(a) The 95% credible region

2

2.5
#10!3

(b) The measurements and the stress-strain curves

Figure 13: Linear elastic-linear hardening: The 95% credible region and the posterior distribution (a) and
the measurements and the stress-strain associated with the 95% credible region (b)
5.1.4. Linear elastic-nonlinear hardening (LE-NH)
Identification of the material parameters. For this subsection twelve measurements are generated based on
a specimen with E = 210 GPa, σy0 = 0.25 GPa, H = 2 GPa, n = 0.57 and the same noise distribution as in
the previous subsections. The prior distribution is selected in the form of Eq. (64) with the following mean

29

vector and covariance matrix:




200
2500


 0
0.29
 GPa, Γx = 
x=


 0
2.5
−1
0.57 GPa
0


0
0

0
0
 GPa2 .

0.1111
0
−2
0
0.0025 GPa

0
2.7778 × 10−4
0
0

(110)

Running the adaptive MCMC approach for 104 samples after burning the first 3000 samples yields:


b post
µ

−3.9287 × 10−3
7.7856 × 10−5
−5.2198 × 10−4
1.4961 × 10−4 GPa−1



b post
Γ


209.8035


0.2566
 GPa,
=


2.1602
0.6034 GPa−1

25.6975

−3.9287 × 10−3
=

−9.0948 × 10−2
−6.8572 × 10−3 GPa−1

−9.0948 × 10−2
−5.2198 × 10−4
9.6433 × 10−2
5.7512 × 10−3 GPa−1


−6.8572 × 10−3 GPa−1
−1
1.4961 × 10−4 GPa 
 GPa2 ,
5.7512 × 10−3 GPa−1 
9.6817 × 10−4 GPa−2
(111)

and

209.4125


\ =  0.2551  GPa.
MAP
 2.1266 
0.597 GPa−1


(112)

The possible stress-strain responses using the 95% credible region for the posterior are presented in Fig. 14.
0.4
0.35
0.3

< (GPa)

0.25
0.2
0.15
0.1
95% credible region curves
Measured value
Mean

0.05
0

0

1

2

3
0

4

5
#10!3

Figure 14: Liner elastic-nonlinear hardening: The measurements and the stress-strain curves associated with
the 95% credible region.
5.2. Bayesian inference with noise in stress and strain
All the measurements used in subsection 5.1 are assumed to be exact in the strain. The aim of current
subsection is to show how the previous results change when the strain is also contaminated by noise. To
accurately investigate the influence of this, exactly the same measurements are used as in the previous
subsection.
30

5.2.1. Linear elastic (LE)
Identification of the Young’s modulus. In the first example the same specimen and the same measurements
as in subsection 5.1.1 (E = 210 GPa) are considered. The noise distribution obtained from the ‘calibration
experiments’ is a normal distribution in the form of Eq. (72) with Sσ = 0.01 GPa and S = 0.0001. For
the same measurement pair with σ m = 0.1576 GPa and m = 7.25 × 10−4 , the posterior distribution can be
calculated using Eq. (78). The same prior is selected as in the first example of subsection 5.1.1 (E = 150GPa
and SE = 50 GPa) and the upper limit for the strain is infinity. Running the adaptive MCMC approach
for 104 samples whilst burning the first 3000 samples leads to a posterior with µ
bpost = 202.7767 GPa,
b
\
Spost = 24.1867 GPa and MAP = 197.5282 GPa.
Fig. 15 shows this posterior, together with the posterior of subsection 5.1.1 when only the noise in the
stress is present. It can clearly be seen that the newly established posterior is wider and it does not have the
form of a (modified) normal distribution, as the posterior when only the noise in the stress is considered.
The stress-strain responses created using the 95% credible region are shown in Fig. 16(a). For comparison,
the same responses are shown in Fig. 16(b) when only the stress is polluted by noise. The estimated values
for this case were µpost = 209.1364 GPa, Spost = 9.6642 GPa and MAP = 207.9963 GPa.

0.8
:(E)
:(Ej<m ; 0m )
:(Ej<m )
Least squares method

0.7
0.6

:

0.5
0.4
0.3
0.2
0.1
0

0

100

200

300

400

500

E (GPa)

Figure 15: Linear elastic: The prior and posterior when both the stress and the strain are corrupted by noise
(black dots), the posterior when only the stress is contaminated (red) and the value predicted by the least
squares method (blue dashed). The posterior for the case with noise in the stress and strain is wider than
the posterior for the case with noise in the stress only. It furthermore does not have the form of a (modified)
normal. The distributions are not normalised.

31

0.3

0.3
95% credible region curves
Measured value
Mean

0.25

0.25
0.2

< (GPa)

0.2
0.15

0.15

0.1

0.1

0.05

0.05

0

95% credible region curves
Measured value
Mean

0

0.2

0.4

0.6

0.8

1

0

1.2

0

0.2

0.4

#10!3

(a) Uncertainty in both the stress and the strain
measurements

0.6
0

0.8

1

1.2
#10!3

(b) Uncertainty in the stress only

Figure 16: Linear elastic: The measurements, the stress-strain curves associated with the mean and the 95%
credible region for (a) noise in the stress and strain and (b) noise in the stress only. The uncertainty is larger
for the case with noise in the stress and strain than that for the case with noise in the stress only.

0.3

0.3
95% credible region curves
Measured value
Mean

0.25

0.25
0.2

< (GPa)

< (GPa)

0.2
0.15

0.15

0.1

0.1

0.05

0.05

0

95% credible region curves
Measured value
Mean

0

0.2

0.4

0.6
0

0.8

1

0

1.2
#10!3

(a) Uncertainty in both the stress and the strain measurements

0

0.2

0.4

0.6
0

0.8

1

1.2
#10!3

(b) Uncertainty in the stress only

Figure 17: Linear elastic: The measurements, and the stress-strain curves associated with the mean and the
95% credible region for ten measurements (a) noise in the stress and the strain and (b) noise in the stress
only. The uncertainty is larger for the case with noise in the stress and strain than that for the case with
noise in the stress only.
5.2.2. Linear elastic-perfectly plastic (LE-PP)
Identification of the material parameters. Similar to the example of subsection 5.1.2, in this subsection a
specimen is considered with a Young’s modulus E = 210 GPa and an initial yield stress σy0 = 0.25 GPa.
Furthermore, to be able to compare the current results with those of subsection 5.1.2 the same measurements
32

are considered again, together with the same noise distribution as in the previous subsection. The prior
distribution is furthermore selected as in Eq. (51) with its properties given in Eq. (99). The resulting
posterior is of the form of Eq. (82) with Eq. (81) as the likelihood function. Similar to subsection 5.1.2, the
upper limit for the theoretical strain (a in Eqs. (70) and (81)) is infinite. Employing the adaptive MCMC
approach for 104 samples and burning the first 3000 samples results in:




215.3918
−6.751 × 10−4
204.0458
b
b post =
µ
GPa, Γpost =
GPa2 ,
(113)
0.2572
−6.751 × 10−4 1.4978 × 10−5
and


202.5858
\
MAP =
GPa.
0.2573

(114)

By comparing the posterior’s covariance matrices (Eqs. 100 and 113), one can see that the additional uncerb post )11 . Fig. 18 shows the samples generated
tainty in the strain measurement has a significant effect on (Γ
by the adaptive MCMC approach to approximate the posterior distribution and the associated 95% credible
region.
#10 -4

0.3

95% credible region

#10!3
1

8

0.28

7

0.8
0.6

<y0 (GPa)

:

9

0.4
0.2

6

0.26

5
4

0.24

3

0
0.3

0.22
0.25
<y0 (GPa)

2

250

1

200
0.2

150

0.2

E (GPa)

(a) Samples generated by the adaptive MCMC approach

160

180

200
220
E (GPa)

240

260

0

(b) 95% credible region

Figure 18: Linear elastic-perfectly plastic: The samples generated by the adaptive MCMC approach to
approximate the posterior distribution (a) and the 95% credible region (b).
The stress-strain curves associated with the credible region are shown in Fig. 19(a). The same curves
for the case with noise in the stress measurements only are presented in Fig. 19(b). Comparing Figs. 19(a)
and 19(b), one can see that the additional uncertainty in the strain has a considerable effect on the elastic
response, but not on the plastic response. This is caused by the fact that the plastic part of the response
does not depend on the strain. Practically no influence on the yield stress can be observed, because this
b post )22 of Eq. 100 and
parameter is independent of the strain. This can also be observed by comparing (Γ
that of Eq. 113.

33

0.3

0.25

0.25

0.2

0.2

< (GPa)

< (GPa)

0.3

0.15
0.1

0.1
95% credible region curves
Measured value
Mean

0.05
0

0.15

0

0.5

1

1.5
0

2

95% credible region curves
Measured value
Mean

0.05
0

2.5

0

0.5

#10!3

1

1.5

2

0

(a) Uncertainty in both the stress and the strain
measurements

2.5
#10!3

(b) Uncertainty in the stress measurements only

Figure 19: Linear elastic-perfectly plastic: The measurements, and the stress-strain responses associated
with the mean and the 95% credible region for twelve measurements: (a) noise in the stress and strain
measurements and (b) noise in the stress measurements only. The additional uncertainty in the strain has a
considerable effect on the elastic response, but not on the plastic response. This is caused by the fact that
the plastic part of the response does not depend on the strain.
5.2.3. Linear elastic-linear hardening (LE-LH)
Identification of the material parameters. In this subsection we will regard the results of the BI formulation
for the linear elastic-linear hardening model, when both the stress and the strain measurements are polluted
by statistical noise. The same noise model and noise distribution is employed as in the previous subsections.
The same measurements as in subsection 5.1.3 are furthermore considered, which are created with E =
210 GPa, σy0 = 0.25 GPa and H = 50 GPa. Similar to in the previous subsections, it is assumed that the
theoretical strain ranges from zero to infinity. Selecting the prior distribution in the form of Eq. (54) with
the mean vector and covariance matrix given by Eq. (107), the resulting posterior is in the form of Eq. (85),
where the likelihood function is given by Eq. (83). Employing the adaptive MCMC approach for 104 sample
generations and burning the first 3000 samples yields:

b post
µ




148.2602
204.269
b post = −5.1006 × 10−2
=  0.2553  GPa, Γ
56.08
−14.8463

−5.1006 × 10−2
6.0511 × 10−5
−2.7164 × 10−2


−14.8463
−2.7164 × 10−2  GPa2 ,
76.7817

(115)

and


201.8193
\ =  55.0417  GPa.
MAP
0.256

(116)

Comparing the results above with those in Eqs. (108) and (109) shows that the uncertainty in the strain
has a larger effect on the Young’s modulus and its corresponding components in the posterior’s covariance
matrix than on the initial yield stress and the plastic modulus. Fig. 20 shows the samples generated by the
adaptive MCMC approach and the associated 95% credible region. The associated stress-strain responses are
presented in Fig. 21. The large influence of the additional noise in the strain measurements on the Young’s
modulus compared to its influence on the plastic parameters can clearly be distinguished.
34

#10 -3

1.5

100

80

H (GPa)

H (GPa)

100

1

60
40

80
60
40

0.5

20
0.3

20
0.3

250

0.25
<y0 (GPa)

0.2

150

250

0.25

200

0

200

<y0 (GPa)

E (GPa)

(a) Samples generated by the adaptive MCMC approach

0.2

E (GPa)

150

(b) 95% credible region

0.35

0.35

0.3

0.3

0.25

0.25

0.2

0.2

< (GPa)

< (GPa)

Figure 20: Linear elastic-linear hardening: Samples generated by the adaptive MCMC approach to approximate the posterior distribution, its projection on the three planes (a) and the 95% credible region (b).

0.15
0.1

0.1
95% credible region curves
Measured value
Mean

0.05
0

0.15

0

0.5

1

1.5
0

2

0

2.5
#10

95% credible region curves
Measured value
Mean

0.05

!3

0

0.5

1

1.5
0

(a) Uncertainty in both the stress and the strain
measurements

2

2.5
#10!3

(b) Uncertainty in the stress measurements only

Figure 21: Linear elastic-linear hardening: The measurements, and the stress-strain curves associated with
the mean and the 95% credible region for twelve measurements: (a) noise in the stress and strain measurements and (b) noise in the stress measurements only. For the case with noise in both the stress and strain
measurements, the uncertainty of the Young’s modulus is larger than for the case with noise in the stress
measurements only. Hardly any difference can be distinguished for the plastic parameters however.
5.2.4. Linear elastic-nonlinear hardening (LE-NH)
Identification of the material parameters. The specimen in this subsection is the same as in subsection
5.1.4 and has a Young’s modulus E = 210 GPa, an initial yield stress σy0 = 0.25 GPa, plastic parameters
H = 2 GPa and n = 0.57. The same measurements as in subsection 5.1.4 are considered, as well as the
same noise distribution. The prior distribution is in the form of Eq. (64) with its values given in Eq. (110).
35

Consequently, the posterior distribution is in the form of Eq. (97) with conditional probability π(σ m |x, m )
given by Eq. (91). Assuming again that the theoretical strain ranges from zero to infinity and running the
adaptive MCMC approach for 104 samples whilst burning the first 3000 samples leads to:


b post
µ


b post
Γ

131.9068
 −2.7875 × 10−2
=

0.4345
4.3498 × 10−3 GPa−1


209.2


0.2557
 GPa,
=


2.1918
−1
0.6033 GPa

−2.7875 × 10−2
8.4839 × 10−5
−4.9175 × 10−4
1.7498 × 10−4 GPa−1

0.4345
−4.9175 × 10−4
9.925 × 10−2
6.1873 × 10−3 GPa−1

4.3498 × 10−3
1.7498 × 10−4
6.1877 × 10−3
1.1031 × 10−3


GPa−1
GPa−1 
 GPa2 ,
GPa−1 
GPa−2
(117)

and



208.6431


0.2555
 GPa.
\ =
MAP


2.2619
−1
0.6055 GPa

(118)

0.4

0.4

0.35

0.35

0.3

0.3

0.25

0.25

< (GPa)

< (GPa)

The associated stress-strain curves are shown in Fig. 22 for both the single uncertainty and the double
uncertainty cases. As can be seen, the additional uncertainty in the strain largely influences the Young’s
modulus because the multiplied gain to the strain (E) in the elastic part is substantially larger than the gain
in plastic part.

0.2
0.15

0.15

0.1

0.1
95% credible region curves
Measured value
Mean

0.05
0

0.2

0

1

2

3
0

4

0

5
#10

95% credible region curves
Measured value
Mean

0.05

!3

0

1

2

3
0

(a) Uncertainty in both the stress and the strain measurements

4

5
#10!3

(b) Uncertainty in the stress measurements only

Figure 22: Linear elastic-nonlinear hardening: The measurements, and the stress-strain curves associated
with the mean and the 95% credible region for twelve measurements: (a) noise in the stress and strain
measurements and (b) noise in the stress measurements only. The uncertainty in the strain largely influences
the Young’s modulus because the multiplied gain to the strain (E) in the elastic part is much larger than
the gain in plastic part.

36

6. Conclusions
In this contribution, we proposed a number of Bayesian inference formulations to identify elastoplastic
material parameters of a single spring based on uniaxial tensile results. The main objectives of the paper
were:
(1) to explain and propose BI formulations for the identification of elasto-plastic material parameters,
when only the stress measurements are contaminated by statistical noise,
(2) to propose BI formulations for the identification of elasto-plastic material parameters, when the stress
measurements as well as the strain measurements are contaminated by statistical noise,
(3) expose possible misconceptions of BI.
In order to apply or develop a Bayesian approach, a noise model and noise distribution need to be
formulated, based on calibration results. Our calibration results were artificially created, but this has enabled
us to to make a direct comparison with the values used to create the measurements. We have furthermore
made comparisons with the results of the least squares method.
The examples given in section 5 allow us to draw four conclusions:
(1) The results of BI cannot directly be compared to those of the least squares method. BI assumes that
other measurements could have been made and hence, the actual measurements are only (a limited
number of) realisations of statistical distributions. BI aims to take this awareness into account (via
the prior distribution).
(2) If one however wants to compare the results of both approaches, it is shown in Fig. 3(a) that the
selected prior distribution has a significant effect on the results. Fig. 3(a) also shows that the influence
of the prior decreases significantly if the number of measurements increases.
(3) The standard deviations and correlations of the material parameters established using the ‘standard’ BI
formulations presented in this contribution, do not reflect the heterogeneity of the material parameters.
In other words, they are not representative for the standard deviations and correlations of the material
parameter distributions. Hence, they are only a measure for the uncertainty of the mean values and
maximum-a-posteriori-probability (MAP) points. If one wants to determine the actual distributions
of the material parameters, the ‘standard’ BI formulations presented here need to be extended by
assuming a probability density function (PDF) for the parameter vector, which is also a source of
uncertainty in addition to measurement error. The objective is then to estimate the representative
parameters of the probability density function (e.g. the mean and variance). This is the focus of our
future work.
(4) Comparing the results when only the stress measurements are polluted by statistical noise to those
when also the strain measurements are polluted by statistical noise, has shown that the uncertainty
of the Young’s modulus is substantially influenced by the additional noise in the strain measurements.
The influence of the additional noise in the strain measurements on the uncertainty of the plastic
parameters is limited however. It may therefore be worthwhile to consider noise in the strains (next
to noise in the stresses) when one is interested in elastic parameters. If one is mostly interested in the
plastic parameters however, it seems questionable if the substantially more complex BI formulations
for double error sources are worth the extra effort.
Acknowledgements
The authors would like to acknowledge the financial support from the University of Luxembourg and the
European Research Council Starting Independent Research Grant (ERC Stg grant agreement No. 279578)
entitled “Towards real time multiscale simulation of cutting in nonlinear materials with applications to
surgical simulation and computer guided surgery”. We also acknowledge the helpful discussion with Prof.
J.T. Oden on the philosophy of Bayesian inference.
37

References
[1] C. Gogu, R. Haftka, R. L. Riche, J. Molimard, A. Vautrin, Introduction to the Bayesian approach
applied to elastic constants identification, AIAA journal 48 (5) (2010) 893–903.
[2] D. Higdon, H. Lee, Z. Bi, A Bayesian approach to characterizing uncertainty in inverse problems using
coarse and fine scale information, IEEE Transactions on Signal Processing 50 (2002) 388–399.
[3] J. Wang, N. Zabaras, A Bayesian inference approach to the inverse heat conduction problem, International Journal of Heat and Mass Transfer 47 (17-18) (2004) 3927–3941.
[4] P. Risholm, F. Janoos, I. Norton, A. J. Golby, W. M. Wells, Bayesian characterization of uncertainty
in intra-subject non-rigid registration, Medical image analysis 17 (5) (2013) 538–555.
[5] S. Lan, T. Bui-Thanh, M. Christie, M. Girolami, Emulation of higher-order tensors in manifold Monte
Carlo methods for Bayesian inverse problems, Journal of Computational Physics 308 (2016) 81–101.
[6] J. L. Beck, L. S. Katafygiotis, Updating models and their uncertainties. I: Bayesian statistical framework,
Journal of Engineering Mechanics 124 (4) (1998) 455–461.
[7] C. K. Oh, J. L. Beck, M. Yamada, Bayesian learning using automatic relevance determination prior
with an application to earthquake early warning, Journal of Engineering Mechanics 134 (12) (2008)
1013–1020.
[8] J. Kaipio, E. Somersalo, Statistical and computational inverse problems, Vol. 160, Springer Science &
Business Media, 2006.
[9] J. Isenberg, Progressing from least squares to Bayesian estimation, in: Proceedings of the 1979 ASME
design engineering technical conference, New York, 1979, pp. 1–11.
[10] K. F. Alvin, Finite element model update via Bayesian estimation and minimization of dynamic residuals, AIAA journal 35 (5) (1997) 879–886.
[11] T. Marwala, S. Sibusiso, Finite element model updating using Bayesian framework and modal properties,
Journal of Aircraft 42 (1) (2005) 275–278.
[12] T. C. Lai, K. H. Ip, Parameter estimation of orthotropic plates by Bayesian sensitivity analysis, Composite Structures 34 (1) (1996) 29–42.
[13] F. Daghia, S. de Miranda, F. Ubertini, E. Viola, Estimation of elastic constants of thick laminated
plates within a Bayesian framework, Composite Structures 80 (3) (2007) 461–473.
[14] P. S. Koutsourelakis, A novel Bayesian strategy for the identification of spatially varying material properties and model validation: An application to static elastography, International Journal for Numerical
Methods in Engineering 91 (3) (2012) 249–268.
[15] C. Gogu, W. Yin, R. Haftka, P. Ifju, J. Molimard, R. Le Riche, A. Vautrin, Bayesian identification of
elastic constants in multi-directional laminate from moiré interferometry displacement fields, Experimental Mechanics 53 (4) (2013) 635–648.
[16] M. Muto, J. L. Beck, Bayesian updating and model class selection for hysteretic structural models using
stochastic simulation, Journal of Vibration and Control 14 (1-2) (2008) 7–34.
[17] P. Liu, S. K. Au, Bayesian parameter identification of hysteretic behavior of composite walls, Probabilistic Engineering Mechanics 34 (2013) 101–109.
[18] D. D. Fitzenz, A. Jalobeanu, S. H. Hickman, Integrating laboratory creep compaction data with numerical fault models: A Bayesian framework, Journal of Geophysical Research: Solid Earth 112 (B8),
b08410.
38

[19] T. Most, Identification of the parameters of complex constitutive models: Least squares minimization
vs. Bayesian updating, in: D. Straub (Ed.), Reliability and optimization of structural systems, CRC
press, 2010, pp. 119–130.
[20] B. V. Rosić, A. Kčerová, J. Sýkora, O. Pajonk, A. Litvinenko, H. G. Matthies, Parameter identification
in a probabilistic setting, Engineering Structures 50 (2013) 179–196, Engineering Structures: Modelling
and Computations (special issue IASS-IACM 2012).
[21] W. P. Hernandez, F. C. L. Borges, D. A. Castello, N. Roitman, C. Magluta, Bayesian inference applied
on model calibration of fractional derivative viscoelastic model, in: V. Steffen Jr, D. A. Rade, W. M.
Bessa (Eds.), DINAME 2015 - Proceedings of the XVII International symposium on dynamic problems
of mechanics, Natal, 2015.
[22] J. M. Nichols, W. A. Link, K. D. Murphy, C. C. Olson, A Bayesian approach to identifying structural
nonlinearity using free-decay response: Application to damage detection in composites, Journal of Sound
and Vibration 329 (15) (2010) 2995–3007.
[23] S. Abhinav, C. S. Manohar, Bayesian parameter identification in dynamic state space models using
modified measurement equations, International Journal of Non-Linear Mechanics 71 (2015) 89–103.
[24] S. Madireddy, B. Sista, K. Vemaganti, A Bayesian approach to selecting hyperelastic constitutive models
of soft tissue, Computer Methods in Applied Mechanics and Engineering 291 (2015) 102–122.
[25] J. T. Oden, E. E. Prudencio, A. Hawkins-Daarud, Selection and assesment of phenomenological models
of tumor growth, Mathematical Models and Methods in Applied Sciences 23 (7) (2013) 1309–1338.
[26] J. Chiachı́o, M. Chiachı́o, A. Saxena, S. Sankararaman, G. Rus, K. Goebel, Bayesian model selection
and parameter estimation for fatigue damage progression models in composites, International Journal
of Fatigue 70 (2015) 361–373.
[27] I. Babuška, Z. Sawlan, M. Scavino, B. Szabó, R. Tempone, Bayesian inference and model comparison for
metallic fatigue data, Computer Methods in Applied Mechanics and Engineering 304 (2016) 171–196.
[28] S. Sarkar, D. S. Kosson, S. Mahadevan, J. C. L. Meeussen, H. van der Sloot, J. R. Arnold, K. G. Brown,
Bayesian calibration of thermodynamic parameters for geochemical speciation modeling of cementitious
materials, Cement and Concrete Research 42 (7) (2012) 889–902.
[29] S. L. Cotter, M. Dashti, J. C. Robinson, A. M. Stuart, Bayesian inverse problems for functions and
applications to fluid mechanics, Inverse Problems 25 (11) (2009) 115008.
[30] S. F. Gull, Bayesian data analysis: Straight-line fitting, in: J. Skilling (Ed.), Maximum entropy and
Bayesian methods, Vol. 36 of Fundamental Theories of Physics, Springer Science & Business Media,
Netherlands, 1989, pp. 511–518.
[31] P. Dellaportas, D. A. Stephens, Bayesian analysis of errors-in-variables regression models, Biometrics
51 (3) (1995) 1085–1095.
[32] R. J. Carroll, K. Roeder, L. Wasserman, Flexible parametric measurement error models, Biometrics
55 (1) (1999) 44–54.
[33] R. Scheines, H. Hoijtink, A. Boomsma, Bayesian estimation and testing of structural equation models,
Psychometrika 64 (1) (1999) 37–52.
[34] B. C. Kelly, Some aspects of measurement error in linear regression of astronomical data, The Astrophysical Journal 665 (2) (2007) 1489.
[35] D. W. Hogg, J. Bovy, D. Lang, Data analysis recipes: Fitting a model to data, ArXiv e-printsarXiv:
1008.4686.
39

[36] L. L. Magorou, F. Bos, F. Rouger, Identification of constitutive laws for wood-based panels by means
of an inverse method, Composites Science and Technology 62 (4) (2002) 591–596.
[37] K. Genovese, L. Lamberti, C. Pappalettere, Improved globallocal simulated annealing formulation for
solving non-smooth engineering optimization problems, International Journal of Solids and Structures
42 (1) (2005) 203–237.
[38] A. van der Vaart, Asymptotic statistics, Cambridge University Press, 1998.
[39] T. J. Ulrych, M. D. Sacchi, A. Woodbury, A Bayes tour of inversion: A tutorial, GEOPHYSICS 66 (1)
(2001) 55–69.
[40] J. C. Simo, T. J. R. Hughes, Computational inelasticity, Springer Science & Business Media, New York,
2000.
[41] J. L. Beck, S. K. Au, Bayesian updating of structural models and reliability using Markov chain Monte
Carlo simulation, Journal of Engineering Mechanics 128 (4) (2002) 380–391.
[42] Y. M. Marzouk, H. N. Najm, L. A. Rahn, Stochastic spectral methods for efficient Bayesian solution of
inverse problems, Journal of Computational Physics 224 (2) (2007) 560–586.
[43] J. Kristensen, N. Zabaras, Bayesian uncertainty quantification in the evaluation of alloy properties with
the cluster expansion method, Computer Physics Communications 185 (11) (2014) 2885–2892.
[44] C. Andrieu, N. De Freitas, A. Doucet, M. I. Jordan, An introduction to MCMC for machine learning,
Machine learning 50 (1-2) (2003) 5–43.
[45] S. Brooks, A. Gelman, G. Jones, X. L. Meng, Handbook of Markov chain Monte Carlo, CRC press,
2011.
[46] A. Gelman, G. O. Roberts, W. R. Gilks, Efficient Metropolis jumping rules, in: J. M. Bernardo, J. O.
Berger, A. P. Dawid, A. F. M. Smith (Eds.), Bayesian statistics, Vol. 5 of Oxford Sci. Publ., Oxford
University Press, New York, 1996, pp. 599–607.
[47] H. Haario, E. Saksman, J. Tamminen, Adaptive proposal distribution for random walk Metropolis
algorithm, Computational Statistics 14 (3) (1999) 375–396.
[48] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Rubin, Bayesian data analysis, Chapman & Hall/CRC Texts
in Statistical Science, Chapman & Hall/CRC, 2003.
[49] S. Hassani, Mathematical methods: For students of physics and related fields, Springer Science &
Business Media, New York, 2008.
[50] S. J. D. Prince, Computer vision: Models learning and inference, Cambridge University Press, 2012.

40

