.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/sampling/mcmc/mcmc_metropolis_hastings.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 or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_sampling_mcmc_mcmc_metropolis_hastings.py: Sampling Rosenbrock distribution using Metropolis-Hastings ============================================================ In this example, the Metropolis-Hastings is employed to generate samples from a Rosenbrock distribution. The method illustrates various aspects of the UQpy :class:`.MCMC` class: - various ways of defining the target pdf to sample from, - definition of input parameters required by the algorithm (proposal_type and proposal_scale for :class:`.MetropolisHastings`), - running several chains in parallel, - call diagnostics functions. .. GENERATED FROM PYTHON SOURCE LINES 16-18 Import the necessary libraries. Here we import standard libraries such as numpy and matplotlib, but also need to import the :class:`.MCMC` class from UQpy. .. GENERATED FROM PYTHON SOURCE LINES 21-28 .. code-block:: default from UQpy.sampling import MetropolisHastings import numpy as np import matplotlib.pyplot as plt import time .. GENERATED FROM PYTHON SOURCE LINES 29-41 Explore various ways of defining the target pdf ----------------------------------------------- Define the Rosenbrock probability density function up to a scale factor. Here the pdf is defined directly in the python script - define the Rosenbrock probability density function up to a scale factor, this function only takes as input parameter the point x where to compute the pdf, - define a pdf function that also takes as argument a set of parameters params, - define a function that computes the log pdf up to a constant. Alternatively, the pdf can be defined in an external file that defines a distribution and its :code:`pdf` or :code:`log_pdf` methods (Rosenbrock.py) .. GENERATED FROM PYTHON SOURCE LINES 44-67 .. code-block:: default def rosenbrock_no_params(x): return np.exp(-(100 * (x[:, 1] - x[:, 0] ** 2) ** 2 + (1 - x[:, 0]) ** 2) / 20) def log_rosenbrock_with_param(x, p): return -(100 * (x[:, 1] - x[:, 0] ** 2) ** 2 + (1 - x[:, 0]) ** 2) / p x = MetropolisHastings(dimension=2, pdf_target=rosenbrock_no_params, burn_length=500, jump=50, n_chains=1, nsamples=500) print(x.samples.shape) plt.figure() plt.plot(x.samples[:, 0], x.samples[:, 1], 'o', alpha=0.5) plt.show() plt.figure() x = MetropolisHastings(dimension=2, log_pdf_target=log_rosenbrock_with_param, burn_length=500, jump=50, n_chains=1, args_target=(20,), nsamples=500) plt.plot(x.samples[:, 0], x.samples[:, 1], 'o') plt.show() .. GENERATED FROM PYTHON SOURCE LINES 68-69 In the following, a custom Rosenbrock distribution is defined and its :code:`log_pdf` method is used. .. GENERATED FROM PYTHON SOURCE LINES 72-92 .. code-block:: default from UQpy.distributions import DistributionND class Rosenbrock(DistributionND): def __init__(self, p=20.): super().__init__(p=p) def pdf(self, x): return np.exp(-(100 * (x[:, 1] - x[:, 0] ** 2) ** 2 + (1 - x[:, 0]) ** 2) / self.parameters['p']) def log_pdf(self, x): return -(100 * (x[:, 1] - x[:, 0] ** 2) ** 2 + (1 - x[:, 0]) ** 2) / self.parameters['p'] log_pdf_target = Rosenbrock(p=20.).log_pdf x = MetropolisHastings(dimension=2, pdf_target=log_pdf_target, burn_length=500, jump=50, n_chains=1, nsamples=500) .. GENERATED FROM PYTHON SOURCE LINES 93-94 In the following, we show that if :code:`burn_length` is set to :math:`0`, the first sample is always the seed. .. GENERATED FROM PYTHON SOURCE LINES 97-111 .. code-block:: default seed = [[0., -1.], [1., 1.]] x = MetropolisHastings(dimension=2, log_pdf_target=log_pdf_target, burn_length=0, jump=50, seed=seed, concatenate_chains=False, nsamples=500) print(x.samples.shape) plt.plot(x.samples[:, 0, 0], x.samples[:, 0, 1], 'o', alpha=0.5) plt.plot(x.samples[:, 1, 0], x.samples[:, 1, 1], 'o', alpha=0.5) plt.show() print(seed) print(x.samples[0, :, :]) .. GENERATED FROM PYTHON SOURCE LINES 112-115 The algorithm-specific parameters for MetropolisHastings are proposal and proposal_is_symmetric ------------------------------------------------------------------------------------------------- The default proposal is standard normal (symmetric). .. GENERATED FROM PYTHON SOURCE LINES 118-142 .. code-block:: default # Define a few proposals to try out from UQpy.distributions import JointIndependent, Normal, Uniform proposals = [JointIndependent([Normal(), Normal()]), JointIndependent([Uniform(loc=-0.5, scale=1.5), Uniform(loc=-0.5, scale=1.5)]), Normal()] proposals_is_symmetric = [True, False, False] fig, ax = plt.subplots(ncols=3, figsize=(16, 4)) for i, (proposal, symm) in enumerate(zip(proposals, proposals_is_symmetric)): print(i) try: x = MetropolisHastings(dimension=2, burn_length=500, jump=100, log_pdf_target=log_pdf_target, proposal=proposal, proposal_is_symmetric=symm, n_chains=1, nsamples=1000) ax[i].plot(x.samples[:, 0], x.samples[:, 1], 'o') except ValueError as e: print(e) print('This last call fails because the proposal is in dimension 1, while the target distribution is' ' in dimension 2') plt.show() .. GENERATED FROM PYTHON SOURCE LINES 143-147 Run several chains in parallel ------------------------------------------------------------------------------- The user can provide the total number of samples :code:`nsamples`, or the number of samples per chain :code:`nsamples_per_chain`. .. GENERATED FROM PYTHON SOURCE LINES 150-169 .. code-block:: default x = MetropolisHastings(dimension=2, log_pdf_target=log_pdf_target, jump=1000, burn_length=500, seed=[[0., 0.], [1., 1.]], concatenate_chains=False, nsamples=100) plt.plot(x.samples[:, 0, 0], x.samples[:, 0, 1], 'o', label='chain 1', alpha=0.5) plt.plot(x.samples[:, 1, 0], x.samples[:, 1, 1], 'o', label='chain 2', alpha=0.5) print(x.samples.shape) plt.legend() plt.show() x = MetropolisHastings(dimension=2, log_pdf_target=log_pdf_target, jump=1000, burn_length=500, seed=[[0., 0.], [1., 1.]], concatenate_chains=False, nsamples_per_chain=100) plt.plot(x.samples[:, 0, 0], x.samples[:, 0, 1], 'o', label='chain 1', alpha=0.5) plt.plot(x.samples[:, 1, 0], x.samples[:, 1, 1], 'o', label='chain 2', alpha=0.5) print(x.samples.shape) plt.legend() plt.show() .. GENERATED FROM PYTHON SOURCE LINES 170-172 Initialize without nsamples... then call run ------------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 175-193 .. code-block:: default t = time.time() x = MetropolisHastings(dimension=2, log_pdf_target=log_pdf_target, jump=1000, burn_length=500, seed=[[0., 0.], [1., 1.]], concatenate_chains=False) print('Elapsed time for initialization: {} s'.format(time.time() - t)) t = time.time() x.run(nsamples=100) print('Elapsed time for running MCMC: {} s'.format(time.time() - t)) print('nburn, jump at first run: {}, {}'.format(x.burn_length, x.jump)) print('total nb of samples: {}'.format(x.samples.shape[0])) plt.plot(x.samples[:, 0, 0], x.samples[:, 0, 1], 'o', label='chain 1') plt.plot(x.samples[:, 1, 0], x.samples[:, 1, 1], 'o', label='chain 2') plt.legend() plt.show() .. GENERATED FROM PYTHON SOURCE LINES 194-197 Run another example with a bivariate distributon with copula dependence - use random_state to always have the same outputs ---------------------------------------------------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 200-213 .. code-block:: default from UQpy.distributions import Normal, Gumbel, JointCopula, JointIndependent dist_true = JointCopula(marginals=[Normal(), Normal()], copula=Gumbel(theta=2.)) proposal = JointIndependent(marginals=[Normal(scale=0.2), Normal(scale=0.2)]) sampler = MetropolisHastings(dimension=2, log_pdf_target=dist_true.log_pdf, proposal=proposal, seed=[0., 0.], random_state=123, nsamples=500) print(sampler.samples.shape) print(np.round(sampler.samples[-5:], 4)) plt.plot(sampler.samples[:, 0], sampler.samples[:, 1], linestyle='none', marker='+') .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.000 seconds) .. _sphx_glr_download_auto_examples_sampling_mcmc_mcmc_metropolis_hastings.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/SURGroup/UQpy/master?urlpath=lab/tree/notebooks/auto_examples/sampling/mcmc/mcmc_metropolis_hastings.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: mcmc_metropolis_hastings.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: mcmc_metropolis_hastings.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_