.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_calibration/bayesian_calibration/plot_bayesian_calibration_flooding.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_calibration_bayesian_calibration_plot_bayesian_calibration_flooding.py: Bayesian calibration of the flooding model ========================================== .. GENERATED FROM PYTHON SOURCE LINES 7-16 Abstract -------- The goal of this example is to present the Bayesian calibration of the :ref:`flooding model`. We use the :class:`~openturns.RandomWalkMetropolisHastings` and :class:`~openturns.Gibbs` classes and simulate a sample of the posterior distribution using :ref:`metropolis_hastings`. .. GENERATED FROM PYTHON SOURCE LINES 19-80 Parameters to calibrate ----------------------- The vector of parameters to calibrate is: .. math:: \vect{\theta} = (K_s,Z_v,Z_m). The variables to calibrate are :math:`(K_s,Z_v,Z_m)` and are set to the following values: .. math:: K_s = 30, \qquad Z_v = 50, \qquad Z_m = 55. Observations ------------ In this section, we describe the statistical model associated with the :math:`n` observations. The errors of the water heights are associated with a gaussian distribution with a zero mean and a standard variation equal to: .. math:: \sigma=0.1. Therefore, the observed water heights are: .. math:: H_i = G(Q_i,K_s,Z_v,Z_m) + \epsilon_i for :math:`i=1,...,n` where .. math:: \epsilon \sim \mathcal{N}(0,\sigma^2) and we make the hypothesis that the observation errors are independent. We consider a sample size equal to: .. math:: n=20. The observations are the couples :math:`\{(Q_i,H_i)\}_{i=1,...,n}`, i.e. each observation is a couple made of the flowrate and the corresponding river height. Variables --------- - :math:`Q` : Input. Observed. - :math:`K_s`, :math:`Z_v`, :math:`Z_m` : Input. Calibrated. - :math:`H`: Output. Observed. Analysis -------- In the description of the :ref:`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. .. GENERATED FROM PYTHON SOURCE LINES 82-84 Generate the observations ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 86-92 .. code-block:: Python from matplotlib import pyplot as plt from openturns.usecases import flood_model import openturns.viewer as viewer import numpy as np import openturns as ot .. GENERATED FROM PYTHON SOURCE LINES 93-94 A basic implementation of the probabilistic model is available in the usecases module : .. GENERATED FROM PYTHON SOURCE LINES 94-96 .. code-block:: Python fm = flood_model.FloodModel() .. GENERATED FROM PYTHON SOURCE LINES 97-107 We define the model :math:`g` which has 4 inputs and one output H. The nonlinear least squares does not take into account for bounds in the parameters. Therefore, we ensure that the output is computed whatever the inputs. The model fails into two situations: * if :math:`K_s<0`, * if :math:`Z_v-Z_m<0`. In these cases, we return an infinite number. .. GENERATED FROM PYTHON SOURCE LINES 109-123 .. code-block:: Python 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] .. GENERATED FROM PYTHON SOURCE LINES 124-128 .. code-block:: Python g = ot.PythonFunction(4, 1, functionFlooding) g = ot.MemoizeFunction(g) g.setOutputDescription(["H (m)"]) .. GENERATED FROM PYTHON SOURCE LINES 129-130 We load the input distribution for :math:`Q`. .. GENERATED FROM PYTHON SOURCE LINES 130-132 .. code-block:: Python Q = fm.Q .. GENERATED FROM PYTHON SOURCE LINES 133-134 Set the parameters to be calibrated. .. GENERATED FROM PYTHON SOURCE LINES 136-143 .. code-block:: Python 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)"]) .. GENERATED FROM PYTHON SOURCE LINES 144-145 We create the joint input distribution. .. GENERATED FROM PYTHON SOURCE LINES 147-149 .. code-block:: Python inputRandomVector = ot.JointDistribution([Q, K_s, Z_v, Z_m]) .. GENERATED FROM PYTHON SOURCE LINES 150-151 Create a Monte-Carlo sample of the output :math:`H`. .. GENERATED FROM PYTHON SOURCE LINES 153-157 .. code-block:: Python nbobs = 20 inputSample = inputRandomVector.getSample(nbobs) outputH = g(inputSample) .. GENERATED FROM PYTHON SOURCE LINES 158-159 Generate the observation noise and add it to the output of the model. .. GENERATED FROM PYTHON SOURCE LINES 161-167 .. code-block:: Python sigmaObservationNoiseH = 0.1 # (m) noiseH = ot.Normal(0.0, sigmaObservationNoiseH) ot.RandomGenerator.SetSeed(0) sampleNoiseH = noiseH.getSample(nbobs) Hobs = outputH + sampleNoiseH .. GENERATED FROM PYTHON SOURCE LINES 168-169 Plot the Y observations versus the X observations. .. GENERATED FROM PYTHON SOURCE LINES 171-173 .. code-block:: Python Qobs = inputSample[:, 0] .. GENERATED FROM PYTHON SOURCE LINES 174-180 .. code-block:: Python graph = ot.Graph("Observations", "Q (m3/s)", "H (m)", True) cloud = ot.Cloud(Qobs, Hobs) graph.add(cloud) view = viewer.View(graph) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_bayesian_calibration_flooding_001.svg :alt: Observations :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_bayesian_calibration_flooding_001.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 181-183 Setting the calibration parameters ---------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 185-190 Define the parametric model :math:`\vect z = f_Q(\vect\theta)` that associates each observation :math:`Q` and value of the parameters :math:`\vect \theta = (K_s, Z_v, Z_m)` to the parameters of the distribution of the corresponding observation: here :math:`\vect z=(\mu, \sigma)` with :math:`\mu = G(Q, K_s, Z_v, Z_m)` and :math:`\sigma = 0.5`. .. GENERATED FROM PYTHON SOURCE LINES 193-204 .. code-block:: Python 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) .. rst-class:: sphx-glr-script-out .. code-block:: none ParametricEvaluation(class=PythonEvaluation name=OpenTURNSPythonFunction, parameters positions=[0], parameters=[x0 : nan], input positions=[1,2,3]) .. GENERATED FROM PYTHON SOURCE LINES 205-207 Define the value of the reference values of the :math:`\vect\theta` parameter. In the Bayesian framework, this is called the mean of the *prior* Gaussian distribution. In the data assimilation framework, this is called the *background*. .. GENERATED FROM PYTHON SOURCE LINES 209-215 .. code-block:: Python KsInitial = 20.0 ZvInitial = 49.0 ZmInitial = 51.0 parameterPriorMean = [KsInitial, ZvInitial, ZmInitial] paramDim = len(parameterPriorMean) .. GENERATED FROM PYTHON SOURCE LINES 216-217 Define the covariance matrix of the parameters :math:`\vect\theta` to calibrate. .. GENERATED FROM PYTHON SOURCE LINES 219-223 .. code-block:: Python sigmaKs = 5.0 sigmaZv = 1.0 sigmaZm = 1.0 .. GENERATED FROM PYTHON SOURCE LINES 224-229 .. code-block:: Python parameterPriorCovariance = ot.CovarianceMatrix(paramDim) parameterPriorCovariance[0, 0] = sigmaKs**2 parameterPriorCovariance[1, 1] = sigmaZv**2 parameterPriorCovariance[2, 2] = sigmaZm**2 .. GENERATED FROM PYTHON SOURCE LINES 230-231 Define the prior distribution :math:`\pi(\vect\theta)` of the parameter :math:`\vect\theta` .. GENERATED FROM PYTHON SOURCE LINES 233-236 .. code-block:: Python prior = ot.Normal(parameterPriorMean, parameterPriorCovariance) prior.setDescription(["Ks", "Zv", "Zm"]) .. GENERATED FROM PYTHON SOURCE LINES 237-242 Define the distribution of observations :math:`\vect{y} | \vect{z}` conditional on model predictions. Note that its parameter dimension is the one of :math:`\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. .. GENERATED FROM PYTHON SOURCE LINES 244-246 .. code-block:: Python conditional = ot.Normal() .. GENERATED FROM PYTHON SOURCE LINES 247-250 The proposed steps for :math:`K_s` :math:`Z_v` and :math:`Z_m` will all follow uniform distributions, but with different supports. .. GENERATED FROM PYTHON SOURCE LINES 252-254 .. code-block:: Python proposal = [ot.Uniform(-5.0, 5.0), ot.Uniform(-1.0, 1.0), ot.Uniform(-1.0, 1.0)] .. GENERATED FROM PYTHON SOURCE LINES 255-257 Build a Gibbs sampler --------------------- .. GENERATED FROM PYTHON SOURCE LINES 257-267 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 268-269 Generate a sample from the posterior distribution of :math:`\vect \theta`. .. GENERATED FROM PYTHON SOURCE LINES 271-274 .. code-block:: Python sampleSize = 1000 sample = sampler.getSample(sampleSize) .. GENERATED FROM PYTHON SOURCE LINES 275-276 Look at the acceptance rates of the random walk Metropolis-Hastings samplers. .. GENERATED FROM PYTHON SOURCE LINES 276-279 .. code-block:: Python [mh.getAcceptanceRate() for mh in sampler.getMetropolisHastingsCollection()] .. rst-class:: sphx-glr-script-out .. code-block:: none [0.361, 0.39, 0.389] .. GENERATED FROM PYTHON SOURCE LINES 280-281 Build the distribution of the posterior by kernel smoothing. .. GENERATED FROM PYTHON SOURCE LINES 281-285 .. code-block:: Python kernel = ot.KernelSmoothing() posterior = kernel.build(sample) .. GENERATED FROM PYTHON SOURCE LINES 286-287 Display prior vs posterior for each parameter. .. GENERATED FROM PYTHON SOURCE LINES 287-335 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 336-337 sphinx_gallery_thumbnail_number = 2 .. GENERATED FROM PYTHON SOURCE LINES 337-345 .. code-block:: Python grid = plot_bayesian_prior_vs_posterior_pdf(prior, posterior) viewer.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) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_bayesian_calibration_flooding_002.svg :alt: Bayesian calibration :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_bayesian_calibration_flooding_002.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 346-347 .. code-block:: Python viewer.View.ShowAll() .. _sphx_glr_download_auto_calibration_bayesian_calibration_plot_bayesian_calibration_flooding.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_bayesian_calibration_flooding.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_bayesian_calibration_flooding.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_bayesian_calibration_flooding.zip `