.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_calibration/bayesian_calibration/plot_imh_python_distribution.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_imh_python_distribution.py: Sampling from an unnormalized probability density ================================================= We sample from an unnormalized probability density using Metropolis-Hastings algorithms. .. math:: f(x) = \frac{1}{2} (2 + \sin(x)^2) \exp \left[- \left(2 + \cos(3x)^3 + \sin(2x)^3 \right) x \right] \mathbf{1}_{[0, 2 \pi]}(x). This example is drawn from [1]. .. GENERATED FROM PYTHON SOURCE LINES 13-15 Draw the unnormalized probability density ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 15-30 .. code-block:: default import openturns as ot from openturns.viewer import View from numpy import pi ot.RandomGenerator.SetSeed(1) f = ot.SymbolicFunction( "x", "0.5 * (2 + sin(x)^2) * exp( -( 2 + cos(3*x)^3 + sin(2*x)^3) * x )") lower_bound = 0.0 upper_bound = 2.0 * pi graph = f.draw(lower_bound, upper_bound, 100) graph.setTitle("Christian Robert tough density") graph.setXTitle("") graph.setYTitle("") _ = View(graph) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_001.png :alt: Christian Robert tough density :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 31-37 Independent Metropolis-Hastings ------------------------------- Let us use a mixture distribution to approximate the target distribution. This approximation will serve as the instrumental distribution in the independent Metropolis-Hastings algorithm. .. GENERATED FROM PYTHON SOURCE LINES 37-42 .. code-block:: default exp = ot.Exponential(1.0) unif = ot.Normal(5.3, 0.4) instrumentalDistribution = ot.Mixture([exp, unif], [0.9, 0.1]) .. GENERATED FROM PYTHON SOURCE LINES 43-44 Compare the instrumental density to the target density. .. GENERATED FROM PYTHON SOURCE LINES 44-53 .. code-block:: default graph = f.draw(lower_bound, upper_bound, 100) graph.setTitle("Instrumental PDF") graph.setXTitle("") graph.setYTitle("") graph.add(instrumentalDistribution.drawPDF(lower_bound, upper_bound, 100)) graph.setLegendPosition("topright") graph.setLegends(["Unnormalized target density", "Instrumental PDF"]) _ = View(graph) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_002.png :alt: Instrumental PDF :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 54-55 :class:`~MetropolisHastings` and derived classes can work directly with the logarithm of the target density. .. GENERATED FROM PYTHON SOURCE LINES 55-58 .. code-block:: default log_density = ot.ComposedFunction(ot.SymbolicFunction("x", "log(x)"), f) .. GENERATED FROM PYTHON SOURCE LINES 59-60 In this case, it is easier to directly write it as a :class:`~openturns.SymbolicFunction`. .. GENERATED FROM PYTHON SOURCE LINES 60-69 .. code-block:: default log_density = ot.SymbolicFunction( "x", "log(2 + sin(x)^2) - (2 + cos(3*x)^3 + sin(2*x)^3) * x") initialState = ot.Point([3.0]) # not important in this case support = ot.Interval([lower_bound], [upper_bound]) independentMH = ot.IndependentMetropolisHastings( log_density, support, initialState, instrumentalDistribution, [0]) .. GENERATED FROM PYTHON SOURCE LINES 70-71 Get a sample .. GENERATED FROM PYTHON SOURCE LINES 71-76 .. code-block:: default sampleSize = 500 sample = independentMH.getSample(sampleSize) .. GENERATED FROM PYTHON SOURCE LINES 77-78 Plot the PDF of the sample to compare it to the target density .. GENERATED FROM PYTHON SOURCE LINES 78-91 .. code-block:: default kernel = ot.KernelSmoothing() posterior = kernel.build(sample) graph = ot.Graph("Independent Metropolis-Hastings sample: {} points".format( sampleSize), "", "", True, "topright") graph.setBoundingBox(ot.Interval( [lower_bound, 0.0], [upper_bound, f([0.0])[0]])) graph.add(f.draw(lower_bound, upper_bound, 100)) graph.add(posterior.drawPDF(lower_bound, upper_bound, 100)) graph.setColors(ot.Drawable.BuildDefaultPalette(2)) graph.setLegends(["Unnormalized target density", "Sample PDF"]) _ = View(graph) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_003.png :alt: Independent Metropolis-Hastings sample: 500 points :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 92-96 Even with very few sampling points (500), independent Metropolis-Hastings (with a judiciously chosen instrumental distribution) manages to capture the main features of the target density. .. GENERATED FROM PYTHON SOURCE LINES 98-102 Random walk Metropolis-Hastings ------------------------------- Let us use a normal instrumental distribution. .. GENERATED FROM PYTHON SOURCE LINES 102-107 .. code-block:: default instrumentalDistribution = ot.Normal(0.0, 2.5) randomwalkMH = ot.RandomWalkMetropolisHastings( log_density, support, initialState, instrumentalDistribution, [0]) .. GENERATED FROM PYTHON SOURCE LINES 108-109 Get a sample .. GENERATED FROM PYTHON SOURCE LINES 109-111 .. code-block:: default sample = randomwalkMH.getSample(sampleSize) .. GENERATED FROM PYTHON SOURCE LINES 112-113 Plot the PDF of the sample to compare it to the target density .. GENERATED FROM PYTHON SOURCE LINES 113-127 .. code-block:: default kernel = ot.KernelSmoothing() posterior = kernel.build(sample) graph = ot.Graph("Random walk Metropolis-Hastings sample: {} points".format( sampleSize), "", "", True, "topright") graph.setBoundingBox(ot.Interval( [lower_bound, 0.0], [upper_bound, f([0.0])[0]])) graph.add(f.draw(lower_bound, upper_bound, 100)) graph.add(posterior.drawPDF(lower_bound, upper_bound, 100)) graph.setColors(ot.Drawable.BuildDefaultPalette(2)) graph.setLegends(["Unnormalized target density", "Sample PDF"]) _ = View(graph) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_004.png :alt: Random walk Metropolis-Hastings sample: 500 points :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 128-131 References ----------- [1] Marin , J.M. & Robert, C.P. (2007). *Bayesian Core: A Practical Approach to Computational Bayesian Statistics*. Springer-Verlag, New York .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.295 seconds) .. _sphx_glr_download_auto_calibration_bayesian_calibration_plot_imh_python_distribution.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_imh_python_distribution.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_imh_python_distribution.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_