.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/sampling/mcmc/mcmc_algorithm_comparison.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_algorithm_comparison.py: Comparison of various MCMC algorithms ============================================================ .. GENERATED FROM PYTHON SOURCE LINES 9-15 This script illustrates performance of various MCMC algorithms currently integrated in UQpy: - Metropolis Hastings (MH) - Modified Metropolis Hastings (MMH) - Affine Invariant with Stretch moves (Stretch) - Adaptive Metropolis with delayed rejection (DRAM) .. GENERATED FROM PYTHON SOURCE LINES 20-21 Import the necessary libraries. .. GENERATED FROM PYTHON SOURCE LINES 24-33 .. code-block:: default import time import matplotlib.pyplot as plt import numpy as np from UQpy.distributions import Normal, Gumbel, JointCopula, JointIndependent, Uniform from UQpy.sampling import MetropolisHastings, Stretch, ModifiedMetropolisHastings, DREAM, DRAM .. GENERATED FROM PYTHON SOURCE LINES 34-38 Affine invariant with Stretch moves ----------------------------------- This algorithm requires as seed a few samples near the region of interest. Here :class:`.MetropolisHastings` is first run to obtain few samples, used as seed within the :class:`.Stretch` algorithm. .. GENERATED FROM PYTHON SOURCE LINES 40-60 .. code-block:: default def log_Rosenbrock(x): return (-(100 * (x[:, 1] - x[:, 0] ** 2) ** 2 + (1 - x[:, 0]) ** 2) / 20) x = MetropolisHastings(dimension=2, burn_length=0, jump=10, n_chains=1, log_pdf_target=log_Rosenbrock, nsamples=5000) print(x.samples.shape) plt.figure() plt.plot(x.samples[:, 0], x.samples[:, 1], 'o') plt.show() x = Stretch(burn_length=0, jump=10, log_pdf_target=log_Rosenbrock, seed=x.samples[:10].tolist(), scale=2., nsamples=5000) print(x.samples.shape) plt.figure() plt.plot(x.samples[:, 0], x.samples[:, 1], 'o') plt.show() .. GENERATED FROM PYTHON SOURCE LINES 61-63 DREAM algorithm: compare with :class:`.MetropolisHastings` (inputs parameters are set as their default values) --------------------------------------------------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 65-83 .. code-block:: default # Define a function to sample seed uniformly distributed in the 2d space ([-20, 20], [-4, 4]) prior_sample = lambda nsamples: np.array([[-2, -2]]) + np.array([[4, 4]]) * JointIndependent( [Uniform(), Uniform()]).rvs(nsamples=nsamples) fig, ax = plt.subplots(ncols=2, figsize=(12, 4)) seed = prior_sample(nsamples=7) x = MetropolisHastings(dimension=2, burn_length=500, jump=50, seed=seed.tolist(), log_pdf_target=log_Rosenbrock, nsamples=1000) ax[0].plot(x.samples[:, 0], x.samples[:, 1], 'o') x = DREAM(dimension=2, burn_length=500, jump=50, seed=seed.tolist(), log_pdf_target=log_Rosenbrock, nsamples=1000) ax[1].plot(x.samples[:, 0], x.samples[:, 1], 'o') plt.show() .. GENERATED FROM PYTHON SOURCE LINES 84-86 DRAM algorithm ----------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 88-122 .. code-block:: default fig, ax = plt.subplots(ncols=2, figsize=(12, 4)) t = time.time() seed = prior_sample(nsamples=1) x = MetropolisHastings(dimension=2, burn_length=500, jump=10, seed=seed.tolist(), log_pdf_target=log_Rosenbrock, nsamples=1000) ax[0].plot(x.samples[:, 0], x.samples[:, 1], 'o') ax[0].set_xlabel('x1') ax[0].set_ylabel('x2') ax[0].set_title('algorithm: MH') print('time to run MH' + ': {}s'.format(time.time() - t)) x = DRAM(dimension=2, burn_length=500, jump=10, seed=seed.tolist(), log_pdf_target=log_Rosenbrock, save_covariance=True, nsamples=1000) ax[1].plot(x.samples[:, 0], x.samples[:, 1], 'o') ax[1].set_xlabel('x1') ax[1].set_ylabel('x2') ax[1].set_title('algorithm: DRAM') print('time to run DRAM' + ': {}s'.format(time.time() - t)) plt.show() # look at the covariance adaptivity fig, ax = plt.subplots() adaptive_covariance = np.array(x.adaptive_covariance) for i in range(2): ax.plot(np.sqrt(adaptive_covariance[:, 0, i, i]), label='dimension {}'.format(i)) ax.set_title('Adaptive proposal std. dev. in both dimensions') ax.legend() plt.show() .. GENERATED FROM PYTHON SOURCE LINES 123-126 MMH: target pdf is given as a joint pdf ---------------------------------------- The target pdf should be a 1 dimensional distribution or set of 1D distributions. .. GENERATED FROM PYTHON SOURCE LINES 128-141 .. code-block:: default from UQpy.distributions import Normal proposal = [Normal(), Normal()] proposal_is_symmetric = [False, False] x = ModifiedMetropolisHastings(dimension=2, burn_length=500, jump=50, log_pdf_target=log_Rosenbrock, proposal=proposal, proposal_is_symmetric=proposal_is_symmetric, n_chains=1, nsamples=500) fig, ax = plt.subplots() ax.plot(x.samples[:, 0], x.samples[:, 1], linestyle='none', marker='.') .. GENERATED FROM PYTHON SOURCE LINES 142-144 MMH: target pdf is given as a couple of independent marginals ---------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 146-161 .. code-block:: default log_pdf_target = [Normal(loc=0., scale=5.).log_pdf, Normal(loc=0., scale=20.).log_pdf] proposal = [Normal(), Normal()] proposal_is_symmetric = [True, True] x = ModifiedMetropolisHastings(dimension=2, burn_length=100, jump=10, log_pdf_target=log_pdf_target, proposal=proposal, proposal_is_symmetric=proposal_is_symmetric, n_chains=1, nsamples=1000) fig, ax = plt.subplots() ax.plot(x.samples[:, 0], x.samples[:, 1], linestyle='none', marker='.') plt.show() print(x.samples.shape) .. GENERATED FROM PYTHON SOURCE LINES 162-164 Use random_state to provide repeated results ------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 166-178 .. code-block:: default seed = Uniform().rvs(nsamples=2 * 10).reshape((10, 2)) dist_true = JointCopula(marginals=[Normal(), Normal()], copula=Gumbel(theta=2.)) proposal = JointIndependent(marginals=[Normal(scale=0.2), Normal(scale=0.2)]) for _ in range(3): sampler = ModifiedMetropolisHastings(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_algorithm_comparison.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_algorithm_comparison.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: mcmc_algorithm_comparison.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: mcmc_algorithm_comparison.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_