.. 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 unscaled probability density ============================================= .. GENERATED FROM PYTHON SOURCE LINES 7-27 In this example, we show different ways to sample a distribution which PDF is known up to a normalization constant: - Case 1: computing the normalization factor and using a :class:`~openturns.PythonDistribution`, - Case 2: using the Ratio of Uniforms algorithm with the class :class:`~openturns.experimental.RatioOfUniforms` - Case 3: using some Metropolis-Hastings algorithms with the class :class:`~openturns.IndependentMetropolisHastings` and :class:`~openturns.RandomWalkMetropolisHastings`. Consider a distribution whose probability density function :math:`p` is known up to the normalization factor :math:`c`: let :math:`f` be a function such that :math:`p = cf` with :math:`c \in \Rset^+_*`. We illustrate the case with: .. 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). First, we draw the unscaled probability density function. .. GENERATED FROM PYTHON SOURCE LINES 27-50 .. code-block:: Python import openturns as ot import openturns.experimental as otexp import openturns.viewer as otv from math import pi from time import monotonic 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 range_PDF = ot.Interval(lower_bound, upper_bound) graph = f.draw(lower_bound, upper_bound, 512) graph.setTitle( r"Christian Robert function: $f(x) = 0.5(2 + sin^2 x) e^{ -x( 2 + cos^33x + sin^3 2x)}, \quad x \in [0, 2\pi]$" ) graph.setXTitle("x") graph.setYTitle(r"$f(x)$") view = otv.View(graph) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_001.svg :alt: Christian Robert function: $f(x) = 0.5(2 + sin^2 x) e^{ -x( 2 + cos^33x + sin^3 2x)}, \quad x \in [0, 2\pi]$ :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_001.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 51-56 Case 1: Computation of the normalization factor ----------------------------------------------- The best thing to do is to compute the normalization factor thanks to the integration algorithms of the library. We show how to compute the normalization factor using a :class:`~openturns.GaussKronrod` quadrature formula. .. GENERATED FROM PYTHON SOURCE LINES 56-60 .. code-block:: Python denom_fact = ot.GaussKronrod().integrate(f, range_PDF)[0] norm_fact = 1.0 / denom_fact print("normalization factor = ", norm_fact) .. rst-class:: sphx-glr-script-out .. code-block:: none normalization factor = 1.65618702904904 .. GENERATED FROM PYTHON SOURCE LINES 61-67 Thus, we can define the exact PDF expression: .. math:: p(x) = \dfrac{f(x)}{\int_0^{2\pi} f(u)\, du} .. GENERATED FROM PYTHON SOURCE LINES 67-69 .. code-block:: Python exact_PDF = ot.LinearCombinationFunction([f], [norm_fact]) .. GENERATED FROM PYTHON SOURCE LINES 70-77 Then we define a :class:`~openturns.PythonDistribution` which: - *computePDF* is computed from the exact expression, - *computeCDF* is computed using an integration algorithm on the *computePDF*. Doing that way, we use the generic sampler of the distribution, based on the CDF inversion method, implemented in the *getSample* method. .. GENERATED FROM PYTHON SOURCE LINES 77-112 .. code-block:: Python class NewDistribution(ot.PythonDistribution): def __init__(self): super(NewDistribution, self).__init__(1) self.PDF_ = exact_PDF self.logPDF_ = ot.ComposedFunction( ot.SymbolicFunction("x", "log(x)"), self.PDF_ ) self.rangeLow_ = 0.0 self.rangeUp_ = 2.0 * pi def getRange(self): return ot.Interval(self.rangeLow_, self.rangeUp_) def computeLogPDF(self, X): if X[0] < self.rangeLow_ or X[0] >= self.rangeUp_: return -ot.SpecFunc.Infinity return self.logPDF_(X)[0] def computePDF(self, X): if X[0] < self.rangeLow_ or X[0] >= self.rangeUp_: return 0.0 return self.PDF_(X)[0] def computeCDF(self, X): if X[0] < self.rangeLow_: return 0.0 if X[0] >= self.rangeUp_: return 1.0 return ot.GaussLegendre([64]).integrate( self.PDF_, ot.Interval(self.rangeLow_, X[0]) )[0] .. GENERATED FROM PYTHON SOURCE LINES 113-114 We create the new distribution that uses the generic sampler: .. GENERATED FROM PYTHON SOURCE LINES 114-116 .. code-block:: Python newDist_generic = ot.Distribution(NewDistribution()) .. GENERATED FROM PYTHON SOURCE LINES 117-118 We can draw the exact PDF and the PDF estimated from a sample generated by the Ratio of Uniforms algorithm. .. GENERATED FROM PYTHON SOURCE LINES 120-139 .. code-block:: Python size = 500 sample = newDist_generic.getSample(size) ks_algo = ot.KernelSmoothing() ks_algo.setBoundaryCorrection(True) ks_algo.setLowerBound(lower_bound) ks_algo.setUpperBound(upper_bound) ks_pdf_GS = ks_algo.build(sample) g = ks_pdf_GS.drawPDF(-0.5, upper_bound * 1.1, 1001) draw_exact = exact_PDF.draw(lower_bound, upper_bound, 1001).getDrawable(0) draw_exact.setLineWidth(2) g.add(draw_exact) g.setLegends(["GS size = " + str(size), "exact pdf"]) g.setLegendPosition("topright") g.setTitle("Generic sampler (GS)") g.setXTitle("x") view = otv.View(g) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_002.svg :alt: Generic sampler (GS) :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_002.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 140-141 We can change the sampler of that new distribution by implementing the method *getSample* and *getRealization* as follows: .. GENERATED FROM PYTHON SOURCE LINES 141-183 .. code-block:: Python class NewDistribution_RoU(ot.PythonDistribution): def __init__(self): super(NewDistribution_RoU, self).__init__(1) self.PDF_ = exact_PDF self.logPDF_ = ot.ComposedFunction( ot.SymbolicFunction("x", "log(x)"), self.PDF_ ) self.rangeLow_ = 0.0 self.rangeUp_ = 2.0 * pi self.sampler_ = otexp.RatioOfUniforms(self.logPDF_, self.getRange()) def getRange(self): return ot.Interval(self.rangeLow_, self.rangeUp_) def computeLogPDF(self, X): if X[0] < self.rangeLow_ or X[0] >= self.rangeUp_: return -ot.SpecFunc.Infinity return self.logPDF_(X)[0] def computePDF(self, X): if X[0] < self.rangeLow_ or X[0] >= self.rangeUp_: return 0.0 return self.PDF_(X)[0] def computeCDF(self, X): if X[0] < self.rangeLow_: return 0.0 if X[0] >= self.rangeUp_: return 1.0 return ot.GaussLegendre([32]).integrate( self.PDF_, ot.Interval(self.rangeLow_, X[0]) )[0] def getRealization(self): return self.sampler_.getRealization() def getSample(self, n): return self.sampler_.getSample(n) .. GENERATED FROM PYTHON SOURCE LINES 184-185 We create the new distribution that uses the Ratio of Uniforms sampler: .. GENERATED FROM PYTHON SOURCE LINES 185-187 .. code-block:: Python newDist_RoU = ot.Distribution(NewDistribution_RoU()) .. GENERATED FROM PYTHON SOURCE LINES 188-190 We compare the sampling speed of the distribution with the Ratio of Uniforms algorithm to the generic sampling speed. The Ratio of Uniforms algorithm proves to be much quicker than the generic method. .. GENERATED FROM PYTHON SOURCE LINES 190-203 .. code-block:: Python sizeRoU = 10000 sizeGeneric = 100 t0 = monotonic() sample_RoU = newDist_RoU.getSample(sizeRoU) t1 = monotonic() sample_generic = newDist_generic.getSample(sizeGeneric) t2 = monotonic() speedRoU = sizeRoU / (t1 - t0) speedGeneric = sizeGeneric / (t2 - t1) print(f"Is Ratio of Uniforms faster ? {speedRoU > 10 * speedGeneric}") .. rst-class:: sphx-glr-script-out .. code-block:: none Is Ratio of Uniforms faster ? True .. GENERATED FROM PYTHON SOURCE LINES 204-210 Case 2: Direct use the Ratio of Uniforms algorithm -------------------------------------------------- In that case, we want to use the :class:`~openturns.experimental.RatioOfUniforms` algorithm to sample :math:`p`. We need to compute the function :math:`\log f` and its range. We create the function :math:`\log f` and the :class:`~openturns.experimental.RatioOfUniforms`: .. GENERATED FROM PYTHON SOURCE LINES 210-213 .. code-block:: Python log_UnscaledPDF = ot.ComposedFunction(ot.SymbolicFunction("x", "log(x)"), f) ratio_algo = otexp.RatioOfUniforms(log_UnscaledPDF, range_PDF) .. GENERATED FROM PYTHON SOURCE LINES 214-215 We can draw the exact PDF and the PDF estimated from a sample generated by the Ratio of Uniforms algorithm. .. GENERATED FROM PYTHON SOURCE LINES 217-234 .. code-block:: Python size = 100000 sample = ratio_algo.getSample(size) ks_algo = ot.KernelSmoothing() ks_algo.setBoundaryCorrection(True) ks_algo.setLowerBound(lower_bound) ks_algo.setUpperBound(upper_bound) ks_pdf = ks_algo.build(sample) g = ks_pdf.drawPDF(-0.5, upper_bound * 1.1, 1001) g.add(draw_exact) g.setLegends(["RoU size = " + str(size), "exact PDF"]) g.setLegendPosition("topright") g.setTitle("Ratio of Uniforms (RoU) ") g.setXTitle("x") view = otv.View(g) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_003.svg :alt: Ratio of Uniforms (RoU) :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_003.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 235-238 In order to facilitate the comparison with the use of Metropolis-Hastings based algorithms, we generate a sample of the same size and we draw the estimated PDF. .. GENERATED FROM PYTHON SOURCE LINES 238-256 .. code-block:: Python size2 = 500 sample = ratio_algo.getSample(size2) ks_algo = ot.KernelSmoothing() ks_algo.setBoundaryCorrection(True) ks_algo.setLowerBound(lower_bound) ks_algo.setUpperBound(upper_bound) ks_pdf_RoU = ks_algo.build(sample) draw_RoU = ks_pdf_RoU.drawPDF(-0.5, upper_bound * 1.1, 1001).getDrawable(0) draw_RoU.setLineWidth(2) g.add(draw_RoU) g.setLegends( ["RoU size = " + str(size), "exact PDF", "RoU (size = " + str(size2) + ")"] ) g.setTitle("Ratio of Uniforms") view = otv.View(g) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_004.svg :alt: Ratio of Uniforms :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_004.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 257-258 By default, the parameter :math:`r=1`. We can get the associated acceptance ratio: .. GENERATED FROM PYTHON SOURCE LINES 258-260 .. code-block:: Python print("Acceptance ratio = ", ratio_algo.getAcceptanceRatio()) .. rst-class:: sphx-glr-script-out .. code-block:: none Acceptance ratio = 0.10015724687759783 .. GENERATED FROM PYTHON SOURCE LINES 261-263 We can estimate the normalization factor with the Ratio of Uniforms algorithm (see the documenttaion of :class:`~openturns.experimental.RatioOfUniforms`): .. GENERATED FROM PYTHON SOURCE LINES 263-265 .. code-block:: Python print("Normalization factor estimate = ", ratio_algo.getC()) .. rst-class:: sphx-glr-script-out .. code-block:: none Normalization factor estimate = 1.6477858299397723 .. GENERATED FROM PYTHON SOURCE LINES 266-267 We can change the :math:`r` parameter and check the associated acceptance ratio: .. GENERATED FROM PYTHON SOURCE LINES 267-271 .. code-block:: Python r_new = 0.5 ratio_algo.setR(r_new) print("New acceptance ratio = ", ratio_algo.getAcceptanceRatio()) .. rst-class:: sphx-glr-script-out .. code-block:: none New acceptance ratio = 0.11111358030178449 .. GENERATED FROM PYTHON SOURCE LINES 272-278 Case 3(a): Use of 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 278-283 .. code-block:: Python exp = ot.Exponential(1.0) unif = ot.Normal(5.3, 0.4) instrumental_dist = ot.Mixture([exp, unif], [0.9, 0.1]) .. GENERATED FROM PYTHON SOURCE LINES 284-285 Compare the instrumental density to the exact PDF. .. GENERATED FROM PYTHON SOURCE LINES 285-293 .. code-block:: Python graph = instrumental_dist.drawPDF(lower_bound, upper_bound, 512) graph.add(draw_exact) graph.setLegends(["Instrumental PDF", "exact PDF"]) graph.setLegendPosition("upper right") graph.setTitle("Independent Metropolis-Hastings: Instrumental PDF") graph.setXTitle("x") _ = otv.View(graph) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_005.svg :alt: Independent Metropolis-Hastings: Instrumental PDF :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_005.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 294-296 :class:`~openturns.MetropolisHastings` and derived classes can work directly with the logarithm of the unscaled target density. .. GENERATED FROM PYTHON SOURCE LINES 296-304 .. code-block:: Python log_density = ot.ComposedFunction(ot.SymbolicFunction("x", "log(x)"), f) initial_State = ot.Point([3.0]) # not important in this case independent_IMH = ot.IndependentMetropolisHastings( log_density, range_PDF, initial_State, instrumental_dist, [0] ) .. GENERATED FROM PYTHON SOURCE LINES 305-306 Get a sample .. GENERATED FROM PYTHON SOURCE LINES 306-310 .. code-block:: Python sample_Size = 500 sample_IMH = independent_IMH.getSample(sample_Size) .. GENERATED FROM PYTHON SOURCE LINES 311-312 Plot the PDF of the sample to compare it to the target density .. GENERATED FROM PYTHON SOURCE LINES 312-326 .. code-block:: Python ks_algo = ot.KernelSmoothing() ks_algo.setBoundaryCorrection(True) ks_algo.setLowerBound(0.0) ks_algo.setUpperBound(2.0 * pi) posterior_IMH = ks_algo.build(sample_IMH) g_IMH = posterior_IMH.drawPDF(-0.5, upper_bound * 1.1, 1001) g_IMH.add(draw_exact) g_IMH.setLegends(["IMH size = {}".format(sample_Size), "exact PDF"]) g_IMH.setTitle("Independent Metropolis-Hastings (IMH)") g_IMH.setXTitle("x") g_IMH.setLegendPosition("topright") view = otv.View(g_IMH) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_006.svg :alt: Independent Metropolis-Hastings (IMH) :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_006.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 327-331 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 333-336 Case 3(b): Use of Random walk Metropolis-Hastings ------------------------------------------------- Let us use a normal instrumental distribution. .. GENERATED FROM PYTHON SOURCE LINES 336-342 .. code-block:: Python instrumental_dist = ot.Normal(0.0, 2.5) randomwalk_MH = ot.RandomWalkMetropolisHastings( log_density, range_PDF, initial_State, instrumental_dist, [0] ) .. GENERATED FROM PYTHON SOURCE LINES 343-344 Get a sample .. GENERATED FROM PYTHON SOURCE LINES 344-346 .. code-block:: Python sample_RWMH = randomwalk_MH.getSample(sample_Size) .. GENERATED FROM PYTHON SOURCE LINES 347-348 Plot the PDF of the sample to compare it to the target density .. GENERATED FROM PYTHON SOURCE LINES 348-363 .. code-block:: Python ks_algo = ot.KernelSmoothing() ks_algo.setBoundaryCorrection(True) ks_algo.setLowerBound(0.0) ks_algo.setUpperBound(2.0 * pi) posterior_RWMH = ks_algo.build(sample_RWMH) g_RWMH = posterior_RWMH.drawPDF(-0.5, upper_bound * 1.1, 1001) g_RWMH.add(draw_exact) g_RWMH.setLegends(["RWMH size = {}".format(sample_Size), "exact PDF"]) g_RWMH.setTitle("Random walk Metropolis-Hastings (RWMH)") g_RWMH.setXTitle("x") g_RWMH.setLegendPosition("topright") view = otv.View(g_RWMH) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_007.svg :alt: Random walk Metropolis-Hastings (RWMH) :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_007.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 364-368 Final comparison ---------------- We plot on the same graph the estimated PDF with all the previous algorithms with the same sample size. sphinx_gallery_thumbnail_number = 8 .. GENERATED FROM PYTHON SOURCE LINES 368-381 .. code-block:: Python g = ks_pdf_GS.drawPDF(-0.5, upper_bound * 1.1, 1001) g.add(posterior_RWMH.drawPDF(-0.5, upper_bound * 1.1, 1001)) g.add(posterior_IMH.drawPDF(-0.5, upper_bound * 1.1, 1001)) g.add(ks_pdf_RoU.drawPDF(-0.5, upper_bound * 1.1, 1001)) draw_exact.setLineStyle("dashed") draw_exact.setColor("black") g.add(draw_exact) g.setLegends(["GS", "RWMH", "IMH", "RoU", "exact PDF"]) g.setTitle("Comparison of samplers (size = 500)") g.setXTitle("x") view = otv.View(g) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_008.svg :alt: Comparison of samplers (size = 500) :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_imh_python_distribution_008.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 382-384 .. code-block:: Python view.ShowAll() .. GENERATED FROM PYTHON SOURCE LINES 385-388 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 ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_imh_python_distribution.zip `