.. 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 Click :ref:`here ` 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 6-10 Abstract -------- The goal of this example is to present the statistical hypotheses of the bayesian calibration of the :ref:`flooding model`. .. GENERATED FROM PYTHON SOURCE LINES 13-63 Parameters to calibrate ----------------------- The vector of parameters to calibrate is: .. math:: \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. Analysis -------- In this model, the variables :math:`Z_m` and :math:`Z_v` are not identifiables, since only the difference :math:`Z_m-Z_v` matters. Hence, calibrating this model requires some regularization. .. GENERATED FROM PYTHON SOURCE LINES 65-67 Generate the observations ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 69-74 .. code-block:: default import numpy as np import openturns as ot ot.Log.Show(ot.Log.NONE) import openturns.viewer as viewer .. GENERATED FROM PYTHON SOURCE LINES 75-76 A basic implementation of the probabilistic model is available in the usecases module : .. GENERATED FROM PYTHON SOURCE LINES 76-79 .. code-block:: default from openturns.usecases import flood_model as flood_model fm = flood_model.FloodModel() .. GENERATED FROM PYTHON SOURCE LINES 80-88 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 90-102 .. code-block:: default 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 103-107 .. code-block:: default g = ot.PythonFunction(4, 1, functionFlooding) g = ot.MemoizeFunction(g) g.setOutputDescription(["H (m)"]) .. GENERATED FROM PYTHON SOURCE LINES 108-109 We load the input distribution for :math:`Q`. .. GENERATED FROM PYTHON SOURCE LINES 109-111 .. code-block:: default Q = fm.Q .. GENERATED FROM PYTHON SOURCE LINES 112-113 Set the parameters to be calibrated. .. GENERATED FROM PYTHON SOURCE LINES 115-122 .. code-block:: default 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 123-124 We create the joint input distribution. .. GENERATED FROM PYTHON SOURCE LINES 126-128 .. code-block:: default inputRandomVector = ot.ComposedDistribution([Q, K_s, Z_v, Z_m]) .. GENERATED FROM PYTHON SOURCE LINES 129-130 Create a Monte-Carlo sample of the output H. .. GENERATED FROM PYTHON SOURCE LINES 132-136 .. code-block:: default nbobs = 20 inputSample = inputRandomVector.getSample(nbobs) outputH = g(inputSample) .. GENERATED FROM PYTHON SOURCE LINES 137-138 Generate the observation noise and add it to the output of the model. .. GENERATED FROM PYTHON SOURCE LINES 140-146 .. code-block:: default sigmaObservationNoiseH = 0.1 # (m) noiseH = ot.Normal(0.,sigmaObservationNoiseH) ot.RandomGenerator.SetSeed(0) sampleNoiseH = noiseH.getSample(nbobs) Hobs = outputH + sampleNoiseH .. GENERATED FROM PYTHON SOURCE LINES 147-148 Plot the Y observations versus the X observations. .. GENERATED FROM PYTHON SOURCE LINES 150-152 .. code-block:: default Qobs = inputSample[:,0] .. GENERATED FROM PYTHON SOURCE LINES 153-159 .. code-block:: default graph = ot.Graph("Observations","Q (m3/s)","H (m)",True) cloud = ot.Cloud(Qobs,Hobs) graph.add(cloud) view = viewer.View(graph) .. image:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_bayesian_calibration_flooding_001.png :alt: Observations :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 160-162 Setting the calibration parameters ---------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 164-165 Define the parametric model :math:`z = f(x,\theta)` that associates each observation :math:`x_i` and values of the parameters :math:`\theta_i` to the parameters of the distribution of the corresponding observation: here :math:`z=(\mu, \sigma)` .. GENERATED FROM PYTHON SOURCE LINES 167-177 .. code-block:: default def fullModelPy(X): Q, K_s, Z_v, Z_m = X H = g(X)[0] sigmaH = 0.5 # (m^2) The standard deviation of the observation error. return [H,sigmaH] fullModel = ot.PythonFunction(4, 2, fullModelPy) model = ot.ParametricFunction(fullModel, [0], Qobs[0]) model .. raw:: html

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



.. GENERATED FROM PYTHON SOURCE LINES 178-179 Define the value of the reference values of the :math:`\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 181-187 .. code-block:: default KsInitial = 20. ZvInitial = 49. ZmInitial = 51. parameterPriorMean = [KsInitial, ZvInitial, ZmInitial] paramDim = len(parameterPriorMean) .. GENERATED FROM PYTHON SOURCE LINES 188-189 Define the covariance matrix of the parameters :math:`\theta` to calibrate. .. GENERATED FROM PYTHON SOURCE LINES 191-195 .. code-block:: default sigmaKs = 5. sigmaZv = 1. sigmaZm = 1. .. GENERATED FROM PYTHON SOURCE LINES 196-202 .. code-block:: default parameterPriorCovariance = ot.CovarianceMatrix(paramDim) parameterPriorCovariance[0,0] = sigmaKs**2 parameterPriorCovariance[1,1] = sigmaZv**2 parameterPriorCovariance[2,2] = sigmaZm**2 parameterPriorCovariance .. raw:: html

[[ 25 0 0 ]
[ 0 1 0 ]
[ 0 0 1 ]]



.. GENERATED FROM PYTHON SOURCE LINES 203-204 Define the the prior distribution :math:`\pi(\underline{\theta})` of the parameter :math:`\underline{\theta}` .. GENERATED FROM PYTHON SOURCE LINES 206-210 .. code-block:: default prior = ot.Normal(parameterPriorMean,parameterPriorCovariance) prior.setDescription(['Ks', 'Zv', 'Zm']) prior .. raw:: html

Normal(mu = [20,49,51], sigma = [5,1,1], R = [[ 1 0 0 ]
[ 0 1 0 ]
[ 0 0 1 ]])



.. GENERATED FROM PYTHON SOURCE LINES 211-214 Define the distribution of observations :math:`\underline{y} | \underline{z}` conditional on model predictions. Note that its parameter dimension is the one of :math:`\underline{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 216-219 .. code-block:: default conditional = ot.Normal() conditional .. raw:: html

Normal(mu = 0, sigma = 1)



.. GENERATED FROM PYTHON SOURCE LINES 220-221 Proposal distribution: uniform. .. GENERATED FROM PYTHON SOURCE LINES 223-226 .. code-block:: default proposal = [ot.Uniform(-5., 5.),ot.Uniform(-1., 1.),ot.Uniform(-1., 1.)] proposal .. rst-class:: sphx-glr-script-out Out: .. code-block:: none [class=Uniform name=Uniform dimension=1 a=-5 b=5, class=Uniform name=Uniform dimension=1 a=-1 b=1, class=Uniform name=Uniform dimension=1 a=-1 b=1] .. GENERATED FROM PYTHON SOURCE LINES 227-231 Test the MCMC sampler --------------------- The MCMC sampler essentially computes the log-likelihood of the parameters. .. GENERATED FROM PYTHON SOURCE LINES 233-235 .. code-block:: default mymcmc = ot.MCMC(prior, conditional, model, Qobs, Hobs, parameterPriorMean) .. GENERATED FROM PYTHON SOURCE LINES 236-238 .. code-block:: default mymcmc.computeLogLikelihood(parameterPriorMean) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none -146.3520105560049 .. GENERATED FROM PYTHON SOURCE LINES 239-241 Test the Metropolis-Hastings sampler ------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 243-244 - Creation of the Random Walk Metropolis-Hastings (RWMH) sampler. .. GENERATED FROM PYTHON SOURCE LINES 246-248 .. code-block:: default initialState = parameterPriorMean .. GENERATED FROM PYTHON SOURCE LINES 249-252 .. code-block:: default RWMHsampler = ot.RandomWalkMetropolisHastings( prior, conditional, model, Qobs, Hobs, initialState, proposal) .. GENERATED FROM PYTHON SOURCE LINES 253-254 Tuning of the RWMH algorithm. .. GENERATED FROM PYTHON SOURCE LINES 256-257 Strategy of calibration for the random walk (trivial example: default). .. GENERATED FROM PYTHON SOURCE LINES 259-262 .. code-block:: default strategy = ot.CalibrationStrategyCollection(paramDim) RWMHsampler.setCalibrationStrategyPerComponent(strategy) .. GENERATED FROM PYTHON SOURCE LINES 263-264 Other parameters. .. GENERATED FROM PYTHON SOURCE LINES 266-270 .. code-block:: default RWMHsampler.setVerbose(True) RWMHsampler.setThinning(1) RWMHsampler.setBurnIn(200) .. GENERATED FROM PYTHON SOURCE LINES 271-272 Generate a sample from the posterior distribution of the parameters theta. .. GENERATED FROM PYTHON SOURCE LINES 274-277 .. code-block:: default sampleSize = 1000 sample = RWMHsampler.getSample(sampleSize) .. GENERATED FROM PYTHON SOURCE LINES 278-279 Look at the acceptance rate (basic checking of the efficiency of the tuning; value close to 0.2 usually recommended). .. GENERATED FROM PYTHON SOURCE LINES 281-283 .. code-block:: default RWMHsampler.getAcceptanceRate() .. raw:: html

[0.513333,0.623333,0.595833]



.. GENERATED FROM PYTHON SOURCE LINES 284-285 Build the distribution of the posterior by kernel smoothing. .. GENERATED FROM PYTHON SOURCE LINES 287-290 .. code-block:: default kernel = ot.KernelSmoothing() posterior = kernel.build(sample) .. GENERATED FROM PYTHON SOURCE LINES 291-292 Display prior vs posterior for each parameter. .. GENERATED FROM PYTHON SOURCE LINES 294-309 .. code-block:: default from openturns.viewer import View import pylab as pl fig = pl.figure(figsize=(12, 4)) for parameter_index in range(paramDim): graph = posterior.getMarginal(parameter_index).drawPDF() priorGraph = prior.getMarginal(parameter_index).drawPDF() priorGraph.setColors(['blue']) graph.add(priorGraph) graph.setLegends(['Posterior', 'Prior']) ax = fig.add_subplot(1, paramDim, parameter_index+1) _ = ot.viewer.View(graph, figure=fig, axes=[ax]) _ = fig.suptitle("Bayesian calibration") .. image:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_bayesian_calibration_flooding_002.png :alt: Bayesian calibration :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.898 seconds) .. _sphx_glr_download_auto_calibration_bayesian_calibration_plot_bayesian_calibration_flooding.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. 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-jupyter :download:`Download Jupyter notebook: plot_bayesian_calibration_flooding.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_