.. 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 :ref:`Go to the end ` 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-32 .. code-block:: Python 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 33-39 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 39-44 .. code-block:: Python 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 45-46 Compare the instrumental density to the target density. .. GENERATED FROM PYTHON SOURCE LINES 46-55 .. code-block:: Python 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("upper right") 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 56-57 :class:`~MetropolisHastings` and derived classes can work directly with the logarithm of the target density. .. GENERATED FROM PYTHON SOURCE LINES 57-60 .. code-block:: Python log_density = ot.ComposedFunction(ot.SymbolicFunction("x", "log(x)"), f) .. GENERATED FROM PYTHON SOURCE LINES 61-62 In this case, it is easier to directly write it as a :class:`~openturns.SymbolicFunction`. .. GENERATED FROM PYTHON SOURCE LINES 62-73 .. code-block:: Python 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 74-75 Get a sample .. GENERATED FROM PYTHON SOURCE LINES 75-80 .. code-block:: Python sampleSize = 500 sample = independentMH.getSample(sampleSize) .. GENERATED FROM PYTHON SOURCE LINES 81-82 Plot the PDF of the sample to compare it to the target density .. GENERATED FROM PYTHON SOURCE LINES 82-99 .. code-block:: Python kernel = ot.KernelSmoothing() posterior = kernel.build(sample) graph = ot.Graph( "Independent Metropolis-Hastings sample: {} points".format(sampleSize), "", "", True, "upper right", ) 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 100-104 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 106-110 Random walk Metropolis-Hastings ------------------------------- Let us use a normal instrumental distribution. .. GENERATED FROM PYTHON SOURCE LINES 110-116 .. code-block:: Python instrumentalDistribution = ot.Normal(0.0, 2.5) randomwalkMH = ot.RandomWalkMetropolisHastings( log_density, support, initialState, instrumentalDistribution, [0] ) .. GENERATED FROM PYTHON SOURCE LINES 117-118 Get a sample .. GENERATED FROM PYTHON SOURCE LINES 118-120 .. code-block:: Python sample = randomwalkMH.getSample(sampleSize) .. GENERATED FROM PYTHON SOURCE LINES 121-122 Plot the PDF of the sample to compare it to the target density .. GENERATED FROM PYTHON SOURCE LINES 122-140 .. code-block:: Python kernel = ot.KernelSmoothing() posterior = kernel.build(sample) graph = ot.Graph( "Random walk Metropolis-Hastings sample: {} points".format(sampleSize), "", "", True, "upper right", ) 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 141-144 References ----------- [1] Marin , J.M. & Robert, C.P. (2007). *Bayesian Core: A Practical Approach to Computational Bayesian Statistics*. Springer-Verlag, New York .. _sphx_glr_download_auto_calibration_bayesian_calibration_plot_imh_python_distribution.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_imh_python_distribution.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_imh_python_distribution.py `