# 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 :

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 -2.0 -1.4444444444444444 -0.8888888888888888 -0.33333333333333326 0.22222222222222232 0.7777777777777777 1.3333333333333335 1.8888888888888888 2.4444444444444446 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 -9.50794871493506 -3.83296694500105 -2.44545713047953 0.0803625289211318 1.01898069723583 0.661725805623086 -1.57581204592385 -2.95308465670895 -8.8878164296758 -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 1.0 -2.0 4.0 1.0 -1.4444444444444444 2.0864197530864197 1.0 -0.8888888888888888 0.7901234567901234 1.0 -0.33333333333333326 0.11111111111111106 1.0 0.22222222222222232 0.04938271604938276 1.0 0.7777777777777777 0.6049382716049381 1.0 1.3333333333333335 1.7777777777777781 1.0 1.8888888888888888 3.567901234567901 1.0 2.4444444444444446 5.975308641975309 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]:

# It is possible to create a random vector based on the RWMHsampler
myRandomVector = ot.PosteriorRandomVector(RWMHsampler)
myRandomVector.getSample(10)

Out[11]:

v0 v1 v2 0.43934747504736293 0.6058776149007925 -1.827413438026229 0.38285157077930076 1.0952621102581648 -1.827413438026229 0.38285157077930076 0.9687945505321682 -1.827413438026229 0.38285157077930076 0.9687945505321682 -1.827413438026229 0.6339531419468756 0.9687945505321682 -1.827413438026229 0.6339531419468756 0.9687945505321682 -1.827413438026229 0.6339531419468756 0.9687945505321682 -1.827413438026229 0.6339531419468756 0.9687945505321682 -1.827413438026229 0.5463809953171936 0.9687945505321682 -1.827413438026229 0.5463809953171936 0.9687945505321682 -1.827413438026229
In [12]:

# Build the distribution of the posterior by kernel smoothing
sample = myRandomVector.getSample(1000)
kernel = ot.KernelSmoothing()
posterior = kernel.build(sample)

In [13]:

# 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')

Out[13]: