.. 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.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.py: Bayesian calibration of a computer code ======================================= .. GENERATED FROM PYTHON SOURCE LINES 7-97 In this example we compute the parameters of a computer model thanks to Bayesian estimation. We use the :class:`~openturns.RandomWalkMetropolisHastings` and :class:`~openturns.Gibbs` classes and simulate a sample of the posterior distribution using :ref:`metropolis_hastings`. Let us denote :math:`(y_1, \dots, y_n)` the observation sample, :math:`(\vect z_1, \ldots, \vect z_n) = (f(x_1|\vect\theta), \ldots, f(x_n|\vect\theta))` the model prediction, :math:`p(y |\vect z)` the density function of observation :math:`y` conditional on model prediction :math:`\vect z`, and :math:`\vect\theta \in \mathbb{R}^p` the calibration parameters we wish to estimate. The posterior distribution is given by Bayes theorem: .. math:: \pi(\vect\theta | \vect y) \quad \propto \quad L\left(\vect y | \vect\theta\right) \times \pi(\vect\theta) where :math:`\propto` means "proportional to", regarded as a function of :math:`\vect\theta`. The posterior distribution is approximated here by the empirical distribution of the sample :math:`\vect\theta^1, \ldots, \vect\theta^N` 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 .. math:: y_i = \theta_1 + x_i \theta_2 + x_i^2 \theta_3 + \varepsilon_i where :math:`\varepsilon_i \stackrel{i.i.d.}{\sim} \mathcal N(0, 1)`. The "true" value of :math:`\theta` is: .. math:: \vect \theta_{true} = (-4.5,4.8,2.2)^T. We use a normal prior on :math:`\vect\theta`: .. math:: \pi(\vect\theta) = \mathcal N(\vect{\mu}_\vect{\theta}, \mat{\Sigma}_\vect{\theta}) where .. math:: \vect{\mu}_\vect{\theta} = \begin{pmatrix} -3 \\ 4 \\ 1 \end{pmatrix} is the mean of the prior and .. math:: \mat{\Sigma}_\vect{\theta} = \begin{pmatrix} \sigma_{\theta_1}^2 & 0 & 0 \\ 0 & \sigma_{\theta_2}^2 & 0 \\ 0 & 0 & \sigma_{\theta_3}^2 \end{pmatrix} is the prior covariance matrix with .. math:: \sigma_{\theta_1} = 2, \qquad \sigma_{\theta_2} = 1, \qquad \sigma_{\theta_3} = 1.5. The following objects need to be defined in order to perform Bayesian calibration: - The conditional density :math:`p(y|\vect z)` must be defined as a probability distribution. - The computer model must be implemented thanks to the :class:`~openturns.ParametricFunction` class. This takes a value of :math:`\vect\theta` as input, and outputs the vector of model predictions :math:`\vect z`, as defined above (the vector of covariates :math:`\vect x = (x_1, \ldots, x_n)` is treated as a known constant). When doing that, we have to keep in mind that :math:`\vect z` will be used as the vector of parameters corresponding to the distribution specified for :math:`p(y |\vect z)`. For instance, if :math:`p(y|\vect z)` is normal, this means that :math:`\vect z` must be a vector containing the mean and standard deviation of :math:`y`. - The prior density :math:`\pi(\vect\theta)` 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 :math:`\vect\theta` before consideration of the experimental data. Again, this is implemented as a probability distribution. - Metropolis-Hastings algorithm(s), possibly used in tandem with a Gibbs algorithm in order to sample from the posterior distribution of the calibration parameters. This example uses the :class:`~openturns.ParametricFunction` class. Please read its documentation and :doc:`/auto_functional_modeling/vectorial_functions/plot_parametric_function` for a detailed example. .. GENERATED FROM PYTHON SOURCE LINES 99-105 .. code-block:: Python import openturns as ot import openturns.viewer as viewer from matplotlib import pylab as plt ot.Log.Show(ot.Log.NONE) .. GENERATED FROM PYTHON SOURCE LINES 106-107 Dimension of the vector of parameters to calibrate .. GENERATED FROM PYTHON SOURCE LINES 107-111 .. code-block:: Python paramDim = 3 # The number of obesrvations obsSize = 10 .. GENERATED FROM PYTHON SOURCE LINES 112-113 Define the observed inputs :math:`x_i`. .. GENERATED FROM PYTHON SOURCE LINES 115-121 .. code-block:: Python xmin = -2.0 xmax = 3.0 step = (xmax - xmin) / (obsSize - 1) rg = ot.RegularGrid(xmin, step, obsSize) x_obs = rg.getVertices() .. GENERATED FROM PYTHON SOURCE LINES 122-128 Define the parametric model :math:`\vect z = f(x,\vect\theta)` that associates each observation :math:`x` and value of :math:`\vect \theta` to the parameters of the distribution of the corresponding observation :math:`y`: here :math:`\vect z=(\mu, \sigma)` where :math:`\mu`, the first output of the model, is the mean and :math:`\sigma`, the second output of the model, is the standard deviation. .. GENERATED FROM PYTHON SOURCE LINES 130-134 .. code-block:: Python fullModel = ot.SymbolicFunction( ["x", "theta1", "theta2", "theta3"], ["theta1+theta2*x+theta3*x^2", "1.0"] ) .. GENERATED FROM PYTHON SOURCE LINES 135-142 To differentiate between the two classes of inputs (:math:`x` and :math:`\vect\theta`), we define a :class:`~openturns.ParametricFunction` from `fullModel` and make the first input (the observations :math:`x`) its *parameter*: :math:`f_x(\vect \theta) := f(x, \vect \theta)`. We set :math:`x = 1` as a placeholder, but :math:`x` will actually take the values :math:`x_i` of the observations when we sample :math:`\vect\theta`. .. GENERATED FROM PYTHON SOURCE LINES 142-146 .. code-block:: Python linkFunction = ot.ParametricFunction(fullModel, [0], [1.0]) print(linkFunction) .. rst-class:: sphx-glr-script-out .. code-block:: none ParametricEvaluation([x,theta1,theta2,theta3]->[theta1+theta2*x+theta3*x^2,1.0], parameters positions=[0], parameters=[x : 1], input positions=[1,2,3]) .. GENERATED FROM PYTHON SOURCE LINES 147-148 Define the observation noise :math:`\varepsilon {\sim} \mathcal N(0, 1)` and create a sample from it. .. GENERATED FROM PYTHON SOURCE LINES 150-155 .. code-block:: Python ot.RandomGenerator.SetSeed(0) noiseStandardDeviation = 1.0 noise = ot.Normal(0, noiseStandardDeviation) noiseSample = noise.getSample(obsSize) .. GENERATED FROM PYTHON SOURCE LINES 156-158 Define the vector of observations :math:`y_i`, here sampled using the "true" value of :math:`\vect \theta`: :math:`\vect \theta_{true}`. .. GENERATED FROM PYTHON SOURCE LINES 160-162 .. code-block:: Python thetaTrue = [-4.5, 4.8, 2.2] .. GENERATED FROM PYTHON SOURCE LINES 163-168 .. code-block:: Python y_obs = ot.Sample(obsSize, 1) for i in range(obsSize): linkFunction.setParameter(x_obs[i]) y_obs[i, 0] = linkFunction(thetaTrue)[0] + noiseSample[i, 0] .. GENERATED FROM PYTHON SOURCE LINES 169-170 Draw the model predictions vs the observations. .. GENERATED FROM PYTHON SOURCE LINES 172-180 .. code-block:: Python functionnalModel = ot.ParametricFunction(fullModel, [1, 2, 3], thetaTrue) graphModel = functionnalModel.getMarginal(0).draw(xmin, xmax) observations = ot.Cloud(x_obs, y_obs) graphModel.add(observations) graphModel.setLegends(["Model", "Observations"]) graphModel.setLegendPosition("upper left") view = viewer.View(graphModel) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_bayesian_calibration_001.png :alt: y0 as a function of x :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_bayesian_calibration_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 181-184 Define the distribution of observations :math:`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. .. GENERATED FROM PYTHON SOURCE LINES 186-188 .. code-block:: Python conditional = ot.Normal() .. GENERATED FROM PYTHON SOURCE LINES 189-190 Define the mean :math:`\mu_\theta`, the covariance matrix :math:`\Sigma_\theta`, then the prior distribution :math:`\pi(\vect\theta)` of the parameter :math:`\vect\theta`. .. GENERATED FROM PYTHON SOURCE LINES 192-194 .. code-block:: Python thetaPriorMean = [-3.0, 4.0, 1.0] .. GENERATED FROM PYTHON SOURCE LINES 195-203 .. code-block:: Python sigma0 = [2.0, 1.0, 1.5] # standard deviations thetaPriorCovarianceMatrix = ot.CovarianceMatrix(paramDim) for i in range(paramDim): thetaPriorCovarianceMatrix[i, i] = sigma0[i] ** 2 prior = ot.Normal(thetaPriorMean, thetaPriorCovarianceMatrix) prior.setDescription(["theta1", "theta2", "theta3"]) .. GENERATED FROM PYTHON SOURCE LINES 204-207 The proposed steps for :math:`\theta_1`, :math:`\theta_2` and :math:`\theta_3` will all follow a uniform distribution. .. GENERATED FROM PYTHON SOURCE LINES 207-210 .. code-block:: Python proposal = ot.Uniform(-1.0, 1.0) .. GENERATED FROM PYTHON SOURCE LINES 211-213 Test the Metropolis-Hastings sampler ------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 215-217 Creation of a single component random walk Metropolis-Hastings (RWMH) sampler. This involves a combination of the RWMH and the Gibbs algorithms. .. GENERATED FROM PYTHON SOURCE LINES 219-221 .. code-block:: Python initialState = thetaPriorMean .. GENERATED FROM PYTHON SOURCE LINES 222-225 We create a :class:`~openturns.RandomWalkMetropolisHastings` sampler for each component. Each sampler must be aware of the joint prior distribution. We also use the same proposal distribution, but this is not mandatory. .. GENERATED FROM PYTHON SOURCE LINES 225-231 .. code-block:: Python mh_coll = [ ot.RandomWalkMetropolisHastings(prior, initialState, proposal, [i]) for i in range(paramDim) ] .. GENERATED FROM PYTHON SOURCE LINES 232-234 Each sampler must be made aware of the likelihood. Otherwise we would sample from the prior! .. GENERATED FROM PYTHON SOURCE LINES 234-238 .. code-block:: Python for mh in mh_coll: mh.setLikelihood(conditional, y_obs, linkFunction, x_obs) .. GENERATED FROM PYTHON SOURCE LINES 239-240 Finally, the :class:`~openturns.Gibbs` algorithm is constructed from all Metropolis-Hastings samplers. .. GENERATED FROM PYTHON SOURCE LINES 240-243 .. code-block:: Python sampler = ot.Gibbs(mh_coll) .. GENERATED FROM PYTHON SOURCE LINES 244-245 Generate a sample from the posterior distribution of the parameters :math:`\vect \theta`. .. GENERATED FROM PYTHON SOURCE LINES 247-250 .. code-block:: Python sampleSize = 10000 sample = sampler.getSample(sampleSize) .. GENERATED FROM PYTHON SOURCE LINES 251-254 Look at the acceptance rate (basic check of the sampling efficiency: values close to :math:`0.2` are usually recommended for Normal posterior distributions). .. GENERATED FROM PYTHON SOURCE LINES 256-258 .. code-block:: Python [mh.getAcceptanceRate() for mh in sampler.getMetropolisHastingsCollection()] .. rst-class:: sphx-glr-script-out .. code-block:: none [0.282, 0.2944, 0.3035] .. GENERATED FROM PYTHON SOURCE LINES 259-260 Build the distribution of the posterior by kernel smoothing. .. GENERATED FROM PYTHON SOURCE LINES 262-265 .. code-block:: Python kernel = ot.KernelSmoothing() posterior = kernel.build(sample) .. GENERATED FROM PYTHON SOURCE LINES 266-267 Display prior vs posterior for each parameter. .. GENERATED FROM PYTHON SOURCE LINES 267-315 .. 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 316-317 sphinx_gallery_thumbnail_number = 2 .. GENERATED FROM PYTHON SOURCE LINES 317-326 .. 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) plt.show() .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_bayesian_calibration_002.png :alt: Bayesian calibration :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_bayesian_calibration_002.png :class: sphx-glr-single-img .. _sphx_glr_download_auto_calibration_bayesian_calibration_plot_bayesian_calibration.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.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_bayesian_calibration.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_bayesian_calibration.zip `