Bayesian calibration of a computer code¶
In this example we are going to compute the parameters of a computer model thanks to Bayesian estimation.
Let us denote the observation
sample,
the model prediction,
the density function of
observation
conditional on model prediction
, and
the calibration parameters
we wish to estimate.
The posterior distribution is given by Bayes theorem:
(where
means “proportional to”, regarded as a function of
) and is approximated here by the empirical
distribution of the sample
generated by
the Metropolis-Hastings algorithm. This means that any quantity
characteristic of the posterior distribution (mean, variance, quantile,
…) is approximated by its empirical counterpart.
Our model (i.e. the compute code to calibrate) is a standard normal linear regression, where
where and
we use a normal prior on
:
The following objects need to be defined in order to perform Bayesian calibration:
- The conditional density
must be defined as a probability distribution
- The computer model must be implemented thanks to the
ParametricFunction class. This takes a value of
as input, and outputs the vector of model predictions
, as defined above (the vector of covariates
is treated as a known constant). When doing that, we have to keep in mind that
will be used as the vector of parameters corresponding to the distribution specified for
. For instance, if
is normal, this means that
must be a vector containing the mean and variance of
- The prior density
encoding the set of possible values for the calibration parameters, each value being weighted by its a priori probability, reflecting the beliefs about the possible values of
before consideration of the experimental data. Again, this is implemented as a probability distribution
- The Metropolis-Hastings algorithm that samples from the posterior
distribution of the calibration parameters requires a vector
initial values for the calibration parameters, as well as the proposal laws used to update each parameter sequentially.
In [1]:
from __future__ import print_function
import openturns as ot
import math as m
In [2]:
# Number of covariates
covNum = 1
# Dimension of the observations
obsDim = 1
# Dimension of the vector of parameters to calibrate
paramDim = 3
# The number of obesrvations
obsSize = 10
- Define the inputs
In [3]:
x = ot.Sample(obsSize, covNum)
for i in range(obsSize):
x[i, 0] = -2 + 5. * i / 9.
x
Out[3]:
v0 | |
---|---|
0 | -2.0 |
1 | -1.4444444444444444 |
2 | -0.8888888888888888 |
3 | -0.33333333333333326 |
4 | 0.22222222222222232 |
5 | 0.7777777777777777 |
6 | 1.3333333333333335 |
7 | 1.8888888888888888 |
8 | 2.4444444444444446 |
9 | 3.0 |
- Define the vector of observations
In [4]:
y_obs = ot.Sample(obsSize, obsDim)
y_obs[0, 0] = -9.50794871493506
y_obs[1, 0] = -3.83296694500105
y_obs[2, 0] = -2.44545713047953
y_obs[3, 0] = 0.0803625289211318
y_obs[4, 0] = 1.01898069723583
y_obs[5, 0] = 0.661725805623086
y_obs[6, 0] = -1.57581204592385
y_obs[7, 0] = -2.95308465670895
y_obs[8, 0] = -8.8878164296758
y_obs[9, 0] = -13.0812290405651
y_obs
Out[4]:
v0 | |
---|---|
0 | -9.50794871493506 |
1 | -3.83296694500105 |
2 | -2.44545713047953 |
3 | 0.0803625289211318 |
4 | 1.01898069723583 |
5 | 0.661725805623086 |
6 | -1.57581204592385 |
7 | -2.95308465670895 |
8 | -8.8878164296758 |
9 | -13.0812290405651 |
- Define the vector of parameters matching each observation
In [5]:
p = ot.Sample(obsSize, paramDim)
for i in range(obsSize):
for j in range(paramDim):
p[i, j] = (-2 + 5. * i / 9.) ** j
p
Out[5]:
v0 | v1 | v2 | |
---|---|---|---|
0 | 1.0 | -2.0 | 4.0 |
1 | 1.0 | -1.4444444444444444 | 2.0864197530864197 |
2 | 1.0 | -0.8888888888888888 | 0.7901234567901234 |
3 | 1.0 | -0.33333333333333326 | 0.11111111111111106 |
4 | 1.0 | 0.22222222222222232 | 0.04938271604938276 |
5 | 1.0 | 0.7777777777777777 | 0.6049382716049381 |
6 | 1.0 | 1.3333333333333335 | 1.7777777777777781 |
7 | 1.0 | 1.8888888888888888 | 3.567901234567901 |
8 | 1.0 | 2.4444444444444446 | 5.975308641975309 |
9 | 1.0 | 3.0 | 9.0 |
- Define the parametric model
that associates each observation
and values of the parameters
to the parameters of the distribution of the corresponding observation: here
In [6]:
fullModel = ot.SymbolicFunction(
['p1', 'p2', 'p3', 'x1', 'x2', 'x3'], ['p1*x1+p2*x2+p3*x3', '1.0'])# assume sigma=1 is known
# f:(x1, x2, x3)->(z) with parameters (p1, p2, p3)
model = ot.ParametricFunction(fullModel, list(range(paramDim)), [0.0]*paramDim)
- Define the distribution of observations
conditional on model predictions
Note that its parameter dimension is the one of ,
so the model must be adjusted accordingly
In [7]:
conditional = ot.Normal()
conditional
Out[7]:
Normal(mu = 0, sigma = 1)
- Define the the prior distribution
of the parameter
In [8]:
sigma0 = ot.Point(paramDim, 10.0) # standard deviations
C0 = ot.CorrelationMatrix(paramDim) # variance matrix
for i in range(paramDim):
C0[i, i] = sigma0[i] * sigma0[i]
m0 = ot.Point(paramDim, 0.0) # mean
prior = ot.Normal(m0, C0)
prior.setDescription(['p1', 'p2', 'p3'])
prior
Out[8]:
Normal(mu = [0,0,0], sigma = [10,10,10], R = [[ 1 0 0 ]
[ 0 1 0 ]
[ 0 0 1 ]])
In [9]:
# Additional inputs for the random walk Metropolis-Hastings
# sampler constructor:
# 1) Initial state: the prior mean
theta0 = prior.getMean()
# 2) Proposal distribution: uniform
proposal = [ot.Uniform(-1., 1.)] * paramDim
# Creation of the Random Walk Metropolis-Hastings (RWMH) sampler
RWMHsampler = ot.RandomWalkMetropolisHastings(
prior, conditional, model, p, y_obs, theta0, proposal)
In [10]:
# Tuning of the RWMH algorithm:
# 1) Strategy of calibration for the random walk (trivial example: default)
strategy = ot.CalibrationStrategyCollection(paramDim)
RWMHsampler.setCalibrationStrategyPerComponent(strategy)
# 2) Other parameters
RWMHsampler.setVerbose(True)
RWMHsampler.setThinning(1)
RWMHsampler.setBurnIn(2000)
# Ready to generate a sample from the posterior distribution
# of the parameters theta
sampleSize = 10000
sample = RWMHsampler.getSample(sampleSize)
# Look at the acceptance rate
# (basic checking of the efficiency of the tuning;
# value close to 0.2 usually recommended)
print('acceptance rate=', end=' ')
print(RWMHsampler.getAcceptanceRate())
acceptance rate= [0.455417,0.292583,0.13425]
In [11]:
# Build the distribution of the posterior by kernel smoothing
sample = RWMHsampler.getSample(1000)
kernel = ot.KernelSmoothing()
posterior = kernel.build(sample)
In [12]:
# Display prior vs posterior for p1, the first parameter of the model
p_index = 1
graph = posterior.getMarginal(0).drawPDF()
prior_drawable = prior.getMarginal(0).drawPDF().getDrawable(0)
prior_drawable.setColor('blue')
graph.add(prior_drawable)
graph.setLegends(['posterior', 'prior'])
graph.setTitle('p%d parameter calibration' % (p_index+1))
graph
Out[12]: