.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_calibration/bayesian_calibration/plot_ackley_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_ackley_distribution.py: Customize your Metropolis-Hastings algorithm ============================================ This simple example shows how you can build your own variant of the Metropolis-Hastings algorithm. .. GENERATED FROM PYTHON SOURCE LINES 10-18 We want to sample from the distribution with support :math:`[-3, 3]^2` whose PDF :math:`f` is proportional to the Ackley function to the tenth power: .. math:: \forall \vect{x} \in [-3, 3]^2 \quad f(\vect{x}) \propto a(\vect{x})^{10}, where :math:`a` is the Ackey function defined in :ref:`use-case-ackley` page. In the following we call it the "Ackley distribution". .. GENERATED FROM PYTHON SOURCE LINES 18-27 .. code-block:: Python import openturns as ot import openturns.experimental as otexp from openturns.viewer import View from openturns.usecases import ackley_function from numpy import exp, format_float_scientific ot.RandomGenerator.SetSeed(100) .. GENERATED FROM PYTHON SOURCE LINES 28-31 Prepare the Metropolis-Hastings algorithm ----------------------------------------- Define the Ackley distribution support and density (up to a constant factor). .. GENERATED FROM PYTHON SOURCE LINES 31-42 .. code-block:: Python am = ackley_function.AckleyModel() ackley = am.model power10 = ot.SymbolicFunction("x", "x^10") ackley_pdf = ot.ComposedFunction(power10, ackley) logarithm = ot.SymbolicFunction("x", "10 * log(x)") ackley_logpdf = ot.ComposedFunction(logarithm, ackley) lb = -3.0 ub = 3.0 support = ot.Interval([lb] * 2, [ub] * 2) .. GENERATED FROM PYTHON SOURCE LINES 43-47 Define the proposal distribution as a :class:`~openturns.Histogram`. Its ticks (on the X axis of the PDF of the histogram) will remain the same, but its frequencies (on the Y axis) will be updated during the course of the Metropolis-Hastings algorithm. .. GENERATED FROM PYTHON SOURCE LINES 47-53 .. code-block:: Python n_bins = 50 myticks = ot.RegularGrid(lb, (ub - lb) / n_bins, n_bins + 1).getValues() frequencies = [1.0] * (myticks.getSize() - 1) proposal = ot.Histogram(myticks, frequencies) .. GENERATED FROM PYTHON SOURCE LINES 54-62 The state of the Markov chain must be converted to an acceptable set of parameters for the :class:`~openturns.Histogram` distribution. This is the job of the *link function*, which we construct with the :class:`~openturns.OpenTURNSPythonFunction` class. It takes a state as input and outputs the parameters (ticks and frequencies) of the proposal distribution. In our case, the ticks will not depend on the inputs, but the frequencies will be outputs of the Ackley function. .. GENERATED FROM PYTHON SOURCE LINES 62-138 .. code-block:: Python parameter_dim = proposal.getParameter().getSize() parameter_desc = proposal.getParameterDescription() class ConditionalAckley(ot.OpenTURNSPythonFunction): """ When executed, this function returns the parameters of a Histogram which approximates the conditional Ackley distribution obtained when one of the 2 coordinates is fixed. To compute the frequencies of the Histogram, this OpenTURNSPythonFunction computes the values of the Ackley function on a regular grid on the line parallel to either the (1, 0) vector (if the second coordinate is fixed) or the (0, 1) vector (if the first coordinate is fixed) containing the point passed as input. The regular grid covers the part of this line which is contained in the support of the Ackley distribution, implicitly defined as the smallest square that contains the cartesian product of the regular grid with itself. For example, if the regular grid covers the interval [-3, 3], then the support is the square [-3, 3] x [-3, 3]. Parameters ---------- marginal : int The marginal whose value is *not* fixed. If 0, then the line of the regular grid is parallel to the (1, 0) vector. If 1, then the line of the regular grid is parallel to the (0, 1) vector. ticks : RegularGrid Ticks of the Histogram distribution. """ def __init__(self, marginal, ticks): super().__init__(2, parameter_dim) self.setInputDescription(["X0", "X1"]) # input: X0 and X1 coordinates self.setOutputDescription(parameter_desc) # output: parameters of the Histogram self._marginal = marginal # parameter which does not vary after initialization offset = (ticks[1] - ticks[0]) / 2 self._marginal_inputs = ot.Sample.BuildFromPoint(ticks)[0:-1] + offset # _marginal_inputs contains the varying coordinate of the points in the regular grid self._size = self._marginal_inputs.getSize() self._ticks = ticks def _exec(self, X): """ Execute the function on a point X = (X0, X1). Parameters ---------- X : list of 2 floats Point through which the line containing the regular grid passes. Returns ------- parameters : :class:`~openturns.Point` Parameters of the :class:`~openturns.Histogram`. """ inputs = ot.Sample(self._size, X) # sample of inputs for the Ackley function # All input points are initialized at the point X passed as argument. # Replace the varying coordinate with the values of the regular grid. inputs[:, self._marginal] = self._marginal_inputs # Compute the Ackley function on these inputs. outputs = exp(ackley_logpdf(inputs).asPoint()) # The outputs are the unnormalized frequencies of the Histogram # proposal distribution, but the Histogram.setParameter() method # expects a full set of parameters. # The easiest way to provide it is to construct a new Histogram object # with the adequate frequencies and call its getParameter() method. return ot.Histogram(self._ticks, outputs).getParameter() .. GENERATED FROM PYTHON SOURCE LINES 139-149 The 2 components of the state of the Markov chain will be updated one after the other, not simultaneously. We define 2 :class:`~openturns.experimental.UserDefinedMetropolisHastings` algorithms encapsulated within a :class:`~openturns.Gibbs` algorithm, so we need 2 link functions, each corresponding to one of the marginals of the Ackley distribution. Note that thanks to the :class:`~openturns.OpenTURNSPythonFunction` class, we were able to only code one template to be used by two different functions instead of directly coding two functions with the :class:`~openturns.PythonFunction` class. .. GENERATED FROM PYTHON SOURCE LINES 149-153 .. code-block:: Python link_function_0 = ot.Function(ConditionalAckley(0, myticks)) link_function_1 = ot.Function(ConditionalAckley(1, myticks)) .. GENERATED FROM PYTHON SOURCE LINES 154-162 Let us illustrate the first of these functions. We can start by evaluating it at :math:`(0.5, 1.5)`. Let :math:`(X_0, X_1)` be a bidimensional random variable following the Ackley distribution. The output is the set of parameters for a :class:`~openturns.Histogram` distribution that approximates the distribution of :math:`X_0` conditional on :math:`X_1 = 1.5`. Let us update the Histogram we created before with this set of parameters. .. GENERATED FROM PYTHON SOURCE LINES 162-167 .. code-block:: Python par = link_function_0([0.5, 1.5]) proposal.setParameter(par) print(par) .. rst-class:: sphx-glr-script-out .. code-block:: none [-3,0.12,0.574885,0.12,0.590611,0.12,0.616354,0.12,0.575301,0.12,0.474915,0.12,0.345447,0.12,0.217746,0.12,0.127608,0.12,0.0892587,0.12,0.086088,0.12,0.093021,0.12,0.0909145,0.12,0.0769473,0.12,0.0565498,0.12,0.0353399,0.12,0.0193349,0.12,0.0118315,0.12,0.01083,0.12,0.0129736,0.12,0.0148585,0.12,0.0146287,0.12,0.0124841,0.12,0.0091527,0.12,0.00577608,0.12,0.00381056,0.12,0.00381056,0.12,0.00577608,0.12,0.0091527,0.12,0.0124841,0.12,0.0146287,0.12,0.0148585,0.12,0.0129736,0.12,0.01083,0.12,0.0118315,0.12,0.0193349,0.12,0.0353399,0.12,0.0565498,0.12,0.0769473,0.12,0.0909145,0.12,0.093021,0.12,0.086088,0.12,0.0892587,0.12,0.127608,0.12,0.217746,0.12,0.345447,0.12,0.474915,0.12,0.575301,0.12,0.616354,0.12,0.590611,0.12,0.574885]#101 .. GENERATED FROM PYTHON SOURCE LINES 168-171 The `link_function_0` function computes histogram frequencies proportional to the values of the unnormalized Ackley distribution PDF along the :math:`X_1 = 1.5` line. .. GENERATED FROM PYTHON SOURCE LINES 171-188 .. code-block:: Python title = "Ackley PDF (up to a constant factor) and $X_1 = 1.5$ cross-section" graph = ot.Graph(title, "$X_0$", "$X_1$", True) line = ot.Curve([[-3, 1.5], [3, 1.5]], "black", "dashed", 2) graph.add(line) graph.setLegendPosition("upper right") contour = ackley_pdf.draw([lb] * 2, [ub] * 2) reversed_colors = [color for color in reversed(contour.getColors())] contour.setColors(reversed_colors) legends = contour.getLegends() legends = [format_float_scientific(float(v), precision=1) for v in legends[:-1]] contour.setLegends(legends) graph.add(contour) graph.setLegendCorner([1.0, 1.0]) graph.setLegendPosition("upper left") view = View(graph) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_ackley_distribution_001.png :alt: Ackley PDF (up to a constant factor) and $X_1 = 1.5$ cross-section :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_ackley_distribution_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 189-191 Let us compare the histogram to the unnormalized PDF (we still need to rescale it to make it appear properly on the graph). .. GENERATED FROM PYTHON SOURCE LINES 191-202 .. code-block:: Python scaling = ot.SymbolicFunction("x", "1.3e-10 * x") scaled_ackley_pdf = ot.ComposedFunction(scaling, ackley_pdf) graph = proposal.drawPDF() graph.setXTitle("") graph.add(scaled_ackley_pdf.draw(0, 0, [0.0, 1.5], -3.0, 3.0, 100)) graph.setLegends(["Histogram PDF", "Rescaled unnormalized PDF"]) graph.setLegendPosition("top") graph.setTitle("Conditional distribution of $X_0$ given $X_1 = 1.5$") _ = View(graph) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_ackley_distribution_002.png :alt: Conditional distribution of $X_0$ given $X_1 = 1.5$ :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_ackley_distribution_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 203-205 Let us now do the same think with `link_function_1`. This time, the relevant cross-section is along the line :math:`X_1 = 0.5`. .. GENERATED FROM PYTHON SOURCE LINES 205-237 .. code-block:: Python par = link_function_1([0.5, 1.5]) proposal.setParameter(par) print(par) title = "Ackley PDF (up to a constant factor) and $X_0 = 0.5$ cross-section" graph = ot.Graph(title, "$X_0$", "$X_1$", True) line = ot.Curve([[0.5, -3], [0.5, 3]], "black", "dashed", 2) graph.add(line) graph.setLegendPosition("upper right") contour = ackley_pdf.draw([lb] * 2, [ub] * 2) reversed_colors = [color for color in reversed(contour.getColors())] contour.setColors(reversed_colors) legends = contour.getLegends() legends = [format_float_scientific(float(v), precision=1) for v in legends[:-1]] contour.setLegends(legends) graph.add(contour) graph.setLegendCorner([1.0, 1.0]) graph.setLegendPosition("upper left") view = View(graph) scaling = ot.SymbolicFunction("x", "3.1e-10 * x") scaled_ackley_pdf = ot.ComposedFunction(scaling, ackley_pdf) graph = proposal.drawPDF() graph.setXTitle("") graph.add(scaled_ackley_pdf.draw(1, 0, [0.5, 0.0], -3.0, 3.0, 100)) graph.setLegends(["Histogram PDF", "Rescaled unnormalized PDF"]) graph.setLegendPosition("top") graph.setTitle("Conditional distribution of $X_1$ given $X_0 = 0.5$") _ = View(graph) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_ackley_distribution_003.png :alt: Ackley PDF (up to a constant factor) and $X_0 = 0.5$ cross-section :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_ackley_distribution_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_ackley_distribution_004.png :alt: Conditional distribution of $X_1$ given $X_0 = 0.5$ :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_ackley_distribution_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none [-3,0.12,0.714193,0.12,0.709153,0.12,0.714195,0.12,0.635879,0.12,0.49337,0.12,0.331163,0.12,0.187923,0.12,0.0968651,0.12,0.0598784,0.12,0.0526439,0.12,0.0522747,0.12,0.0457978,0.12,0.0334477,0.12,0.0202285,0.12,0.00973413,0.12,0.00378105,0.12,0.00159728,0.12,0.0011016,0.12,0.001075,0.12,0.00097469,0.12,0.00070403,0.12,0.00040715,0.12,0.0001854,0.12,6.70262e-05,0.12,2.77594e-05,0.12,2.77594e-05,0.12,6.70262e-05,0.12,0.0001854,0.12,0.00040715,0.12,0.00070403,0.12,0.00097469,0.12,0.001075,0.12,0.0011016,0.12,0.00159728,0.12,0.00378105,0.12,0.00973413,0.12,0.0202285,0.12,0.0334477,0.12,0.0457978,0.12,0.0522747,0.12,0.0526439,0.12,0.0598784,0.12,0.0968651,0.12,0.187923,0.12,0.331163,0.12,0.49337,0.12,0.635879,0.12,0.714195,0.12,0.709153,0.12,0.714193]#101 .. GENERATED FROM PYTHON SOURCE LINES 238-239 We now choose the initial state of the Markov chain. .. GENERATED FROM PYTHON SOURCE LINES 239-242 .. code-block:: Python initialState = [0.1] * 2 .. GENERATED FROM PYTHON SOURCE LINES 243-247 Sample from the Ackley distribution ----------------------------------- We can finally define the two Metropolis-Hastings algorithms and the Gibbs algorithm which encapsulates them. .. GENERATED FROM PYTHON SOURCE LINES 247-258 .. code-block:: Python gmh_0 = otexp.UserDefinedMetropolisHastings( ackley_logpdf, support, initialState, proposal, link_function_0, [0] ) gmh_1 = otexp.UserDefinedMetropolisHastings( ackley_logpdf, support, initialState, proposal, link_function_1, [1] ) gibbs = ot.Gibbs([gmh_0, gmh_1]) sample = gibbs.getSample(100) .. GENERATED FROM PYTHON SOURCE LINES 259-260 Print the acceptance rates of the two Metropolis-Hastings samplers .. GENERATED FROM PYTHON SOURCE LINES 260-266 .. code-block:: Python mhlist = gibbs.getMetropolisHastingsCollection() rate_gmh_0 = mhlist[0].getAcceptanceRate() rate_gmh_1 = mhlist[1].getAcceptanceRate() print("Acceptance rates: {} and {}".format(rate_gmh_0, rate_gmh_1)) .. rst-class:: sphx-glr-script-out .. code-block:: none Acceptance rates: 0.99 and 0.91 .. GENERATED FROM PYTHON SOURCE LINES 267-269 Draw the contour plot of the Ackley distribution PDF (up to a multiplicative constant) and the MCMC sample: we can see the sample is accurate. .. GENERATED FROM PYTHON SOURCE LINES 269-282 .. code-block:: Python title = "Ackley PDF (up to a constant factor) and Metropolis-Hastings sample" graph = ot.Graph(title, "$X_0$", "$X_1$", True) graph.add(ackley_pdf.draw([lb] * 2, [ub] * 2)) reversed_colors = [color for color in reversed(graph.getColors())] graph.setColors(reversed_colors) graph.add(ot.Cloud(sample, "black", "plus")) legends = graph.getLegends() legends = [format_float_scientific(float(v), precision=1) for v in legends[:-1]] graph.setLegends(legends) graph.setLegendCorner([1.0, 1.0]) graph.setLegendPosition("upper left") view = View(graph) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_ackley_distribution_005.png :alt: Ackley PDF (up to a constant factor) and Metropolis-Hastings sample :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_ackley_distribution_005.png :class: sphx-glr-single-img .. _sphx_glr_download_auto_calibration_bayesian_calibration_plot_ackley_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_ackley_distribution.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_ackley_distribution.py `