Note

Go to the end to download the full example code.

Bayesian calibration of the flooding model¶

Abstract¶

The goal of this example is to present the Bayesian calibration of the flooding model. We use the RandomWalkMetropolisHastings and Gibbs classes and simulate a sample of the posterior distribution using The Metropolis-Hastings Algorithm.

Parameters to calibrate¶

The vector of parameters to calibrate is:

$\vect{\theta} = (K_s,Z_v,Z_m).$

The reference values are:

$K_s = 30, \qquad Z_v = 50, \qquad Z_m = 55.$

Observations¶

We consider the probabilistic model:

$H = \model(Q,K_s,Z_v,Z_m) + \epsilon$

where:

$\epsilon \sim \mathcal{N}(0,\sigma^2)$

Thus, the error of the water height follows a normal distribution with a zero mean and a standard variation $\sigma$ .

We have some observations $\{(Q_i,h_i)\}_{i=1,...,\sampleSize}$ . Each observation is a couple made of the flowrate and the corresponding river height. We assume that the observations are associated to independent errors. We consider a sample size equal to:

$\sampleSize = 10.$

Variables¶

$Q$ : Input. Observed.
$K_s$ , $Z_v$ , $Z_m$ : Input. Calibrated.
$h$ : Output. Observed.

Analysis¶

In the description of the flooding model, we see that only one parameter can be identified. Hence, calibrating this model requires some regularization. In this example, we use Bayesian methods as a way to regularize the model.

Generate the observations¶

from matplotlib import pyplot as plt
from openturns.usecases import flood_model
import openturns.viewer as otv
import numpy as np
import openturns as ot

A basic implementation of the probabilistic model is available in the usecases module :

fm = flood_model.FloodModel()

We define the model $\model$ which has 4 inputs and one output $h$ .

The non linear least squares algorithm does not take into account bounds in the parameters. Therefore, we ensure that the output is computed whatever the inputs. The model fails into two situations:

if $K_s<0$ ,
if $Z_v-Z_m<0$ .

In these cases, we return an infinite number.

def functionFlooding(X):
    L = 5.0e3
    B = 300.0
    Q, K_s, Z_v, Z_m = X
    alpha = (Z_m - Z_v) / L
    if alpha < 0.0 or K_s <= 0.0:
        H = np.inf
    else:
        H = (Q / (K_s * B * np.sqrt(alpha))) ** (3.0 / 5.0)
    return [H]

g = ot.PythonFunction(4, 1, functionFlooding)
g = ot.MemoizeFunction(g)
g.setOutputDescription(["H (m)"])

We load the input distribution for the flow $Q$ .

Q = fm.Q

Set the parameters to be calibrated.

K_s = ot.Dirac(30.0)
Z_v = ot.Dirac(50.0)
Z_m = ot.Dirac(55.0)
K_s.setDescription(["Ks (m^(1/3)/s)"])
Z_v.setDescription(["Zv (m)"])
Z_m.setDescription(["Zm (m)"])

We create the joint input distribution of $(Q, K_s, Z_v, Z_m)$ .

inputRandomVector = ot.JointDistribution([Q, K_s, Z_v, Z_m])

We import some noisy observations of the flow rate (column 0) and the height (column 1) already stored of the flooding model in the field data.

Q_H_obs = fm.data
Qobs = Q_H_obs.getMarginal(0)
Hobs = Q_H_obs.getMarginal(1)

graph = ot.Graph("Observations", "Q (m3/s)", "h (m)", True)
cloud = ot.Cloud(Q_H_obs)
# cloud = ot.Cloud(Qobs, Hobs)
graph.add(cloud)
view = otv.View(graph)

Setting the calibration parameters¶

We assume that the output of the model $H$ is random and follows a normal distribution with the standard deviation $\sigma = 0.5$ and centered on $\mu = \model(Q, K_s, Z_v, Z_m)$ .

Then, we define the parametric model $\vect z = f_Q(\vect\theta)$ that associates each observation $Q$ and value of the parameters $\vect \theta = (K_s, Z_v, Z_m)$ to the parameters $\vect z=(\mu, \sigma)$ . We want to get the posterior distribution of the random vector $\vect \Theta$ that fits the best to the observations $(Q_i, h_i)_i$ .

def fullModelPy(X):
    Q, K_s, Z_v, Z_m = X
    mu = g(X)[0]
    sigma = 0.5  # (m^2) The standard deviation of the observation error.
    return [mu, sigma]


fullModel = ot.PythonFunction(4, 2, fullModelPy)
linkFunction = ot.ParametricFunction(fullModel, [0], [np.nan])
print(linkFunction)

ParametricEvaluation(class=PythonEvaluation name=OpenTURNSPythonFunction, parameters positions=[0], parameters=[x0 : nan], input positions=[1,2,3])

We define the prior distribution of $\vect\Theta$ . We assume that it is a normal distribution. In the Bayesian framework, this is called the prior normal distribution. In the data assimilation framework, this is called the background. We assume that the prior distribution has independent components, that the mean marginal values are equal to the reference values and that the marginal standard deviations are known.

KsInitial = 20.0
ZvInitial = 49.0
ZmInitial = 51.0
parameterPriorMean = [KsInitial, ZvInitial, ZmInitial]
paramDim = len(parameterPriorMean)

sigmaKs = 5.0
sigmaZv = 1.0
sigmaZm = 1.0

parameterPriorCovariance = ot.CovarianceMatrix(paramDim)
parameterPriorCovariance[0, 0] = sigmaKs**2
parameterPriorCovariance[1, 1] = sigmaZv**2
parameterPriorCovariance[2, 2] = sigmaZm**2

Define the prior distribution of $\vect\Theta$ .

prior = ot.Normal(parameterPriorMean, parameterPriorCovariance)
prior.setDescription(["Ks", "Zv", "Zm"])

Define the distribution of the output $H | \vect{z}$ conditional on model predictions. Note that its parameter dimension is the one of $\vect{z}$ , so the model must be adjusted accordingly. In other words, the input argument of the setParameter method of the conditional distribution must be equal to the dimension of the output of the model. Hence, we do not have to set the actual parameters: only the type of distribution is used.

conditional = ot.Normal()

The proposed steps for $K_s$ $Z_v$ and $Z_m$ will all follow uniform distributions, but with different supports.

proposal = [ot.Uniform(-5.0, 5.0), ot.Uniform(-1.0, 1.0), ot.Uniform(-1.0, 1.0)]

Build a Gibbs sampler¶

initialState = parameterPriorMean
mh_coll = [
    ot.RandomWalkMetropolisHastings(prior, initialState, proposal[i], [i])
    for i in range(paramDim)
]
for mh in mh_coll:
    mh.setLikelihood(conditional, Hobs, linkFunction, Qobs)
sampler = ot.Gibbs(mh_coll)

Generate a sample from the posterior distribution of $\vect \Theta$ .

sampleSize = 1000
sample = sampler.getSample(sampleSize)

Look at the acceptance rates of the random walk Metropolis-Hastings samplers.

[mh.getAcceptanceRate() for mh in sampler.getMetropolisHastingsCollection()]

[0.365, 0.361, 0.367]

Build the posterior distribution of $\vect \Theta$ by kernel smoothing.

kernel = ot.KernelSmoothing()
posterior = kernel.build(sample)

Display prior vs posterior for each parameter.

def plot_bayesian_prior_vs_posterior_pdf(prior, posterior):
    """
    Plot the prior and posterior distribution of a Bayesian calibration

    Parameters
    ----------
    prior : ot.Distribution(dimension)
        The prior.
    posterior : ot.Distribution(dimension)
        The posterior.

    Return
    ------
    grid : ot.GridLayout(1, dimension)
        The prior and posterior PDF for each marginal.
    """
    paramDim = prior.getDimension()
    grid = ot.GridLayout(1, paramDim)
    parameterNames = prior.getDescription()
    for parameter_index in range(paramDim):
        graph = ot.Graph("", parameterNames[parameter_index], "PDF", True)
        # Prior
        curve = prior.getMarginal(parameter_index).drawPDF().getDrawable(0)
        curve.setLineStyle(
            ot.ResourceMap.GetAsString("CalibrationResult-PriorLineStyle")
        )
        curve.setLegend("Prior")
        graph.add(curve)
        # Posterior
        curve = posterior.getMarginal(parameter_index).drawPDF().getDrawable(0)
        curve.setLineStyle(
            ot.ResourceMap.GetAsString("CalibrationResult-PosteriorLineStyle")
        )
        curve.setLegend("Posterior")
        graph.add(curve)
        #
        if parameter_index < paramDim - 1:
            graph.setLegends([""])
        if parameter_index > 0:
            graph.setYTitle("")
        graph.setLegendPosition("upper right")
        grid.setGraph(0, parameter_index, graph)
    grid.setTitle("Bayesian calibration")
    return grid

sphinx_gallery_thumbnail_number = 2

grid = plot_bayesian_prior_vs_posterior_pdf(prior, posterior)
otv.View(
    grid,
    figure_kw={"figsize": (8.0, 3.0)},
    legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"},
)
plt.subplots_adjust(right=0.8, bottom=0.2, wspace=0.3)

otv.View.ShowAll()

OpenTURNS

An Open source initiative for the Treatment of Uncertainties, Risks'N Statistics

Table of Contents

Previous topic

Next topic

This Page