.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/sampler/plot_metropolis_hastings.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_auto_examples_sampler_plot_metropolis_hastings.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_sampler_plot_metropolis_hastings.py:


Metropolis Hastings sampling
============================

.. GENERATED FROM PYTHON SOURCE LINES 8-22

This example shows how to draw samples using Metropolis Hastings sampling.
The target probability distribution is

:math:`p(\mathbf{x})=p(x_1, x_2) \propto \exp \left(-\frac{1}{100}\left(x_1-1\right)^2-\left(x_1^2-x_2\right)^2\right)`

where :math:`x_1 \in [-5,5]` and :math:`x_2 \in [-5,5]`.

It should be noted that the right hand side of the equation is an unnormalized
probability density function since its integral is not equal to :math:`1`.
This can happen, for example, when the normalization constant is unknown or difficult
to compute.

We can define the target probability distribution in Python as follows:


.. GENERATED FROM PYTHON SOURCE LINES 23-32

.. code-block:: default


    import numpy as np

    def target_dist(x):
        if (x[0]>=-5 and x[0]<=5) and (x[1]>=-5 and x[1]<=5):
            return np.exp(-0.01*(x[0]-1)**2 - (x[0]**2-x[1])**2)
        else:
            return 0








.. GENERATED FROM PYTHON SOURCE LINES 33-34

The figure below shows how the target distribution looks like.

.. GENERATED FROM PYTHON SOURCE LINES 35-60

.. code-block:: default


    import matplotlib.pyplot as plt
    import itertools

    x1_values = np.linspace(-5,5.1,100)
    x2_values = np.linspace(-5,5.1,100)

    target_values = np.zeros((100, 100))
    for i, j in itertools.product(range(100), range(100)):
        x1 = x1_values[i]
        x2 = x2_values[j]
        target_values[i,j] = target_dist(np.array([x1, x2]))

    fig , ax = plt.subplots(figsize=(4,4))

    ax.contourf(x1_values, x2_values, np.transpose(target_values), levels=10, cmap='Blues')

    ax.set_title('(unnormalized) target distribution')
    ax.set_xlabel(r'$x_1$')
    ax.set_ylabel(r'$x_2$')
    ax.set_xlim(-5,5)
    ax.set_ylim(-5,5)
    plt.tight_layout()





.. image-sg:: /auto_examples/sampler/images/sphx_glr_plot_metropolis_hastings_001.png
   :alt: (unnormalized) target distribution
   :srcset: /auto_examples/sampler/images/sphx_glr_plot_metropolis_hastings_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 61-67

To perform Metropolis Hastings sampling, we need to choose a proposal distribution
which can be used to determine a new state ``x'`` given a current state ``x`` at
each iteration. This is defined by the parameter ``f_sample`` which should be a
function. A usual choice is to choose the new state from a Gaussian distribution
centered at the current state.


.. GENERATED FROM PYTHON SOURCE LINES 68-75

.. code-block:: default


    from scipy.stats import multivariate_normal

    f_sample = multivariate_normal.rvs
    # make the samples reproducible
    kwgs_f_sample = {'random_state': np.random.default_rng(seed=1)}








.. GENERATED FROM PYTHON SOURCE LINES 76-77

Then, we create an instance of :class:`.MetropolisHastings` class.

.. GENERATED FROM PYTHON SOURCE LINES 78-85

.. code-block:: default


    from psimpy.sampler.metropolis_hastings import MetropolisHastings

    mh_sampler = MetropolisHastings(ndim=2, init_state=np.array([-4,-4]),
        f_sample=f_sample, target=target_dist, nburn=1000, nthin=5,
        seed=1, kwgs_f_sample=kwgs_f_sample)








.. GENERATED FROM PYTHON SOURCE LINES 86-88

Next, we call the :py:meth:`.MetropolisHastings.sample` method to draw required 
number of samples.

.. GENERATED FROM PYTHON SOURCE LINES 89-95

.. code-block:: default


    mh_samples, mh_accept = mh_sampler.sample(nsamples=5000)

    print("Acceptance ratio: ", np.mean(mh_accept))






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Acceptance ratio:  0.3712




.. GENERATED FROM PYTHON SOURCE LINES 96-97

The following figure shows how the samples look like.

.. GENERATED FROM PYTHON SOURCE LINES 98-117

.. code-block:: default


    import matplotlib.pyplot as plt
    import itertools

    # sphinx_gallery_thumbnail_number = 2
    fig , ax = plt.subplots(figsize=(4,4))

    ax.contourf(x1_values, x2_values, np.transpose(target_values), levels=10, cmap='Blues')
    ax.scatter(mh_samples[:,0], mh_samples[:,1], s=5, c='r', marker='o',alpha=0.05)

    ax.set_title('Metropolis Hastings samples')
    ax.set_xlabel(r'$x_1$')
    ax.set_ylabel(r'$x_2$')
    ax.set_xlim(-5,5)
    ax.set_ylim(-5,5)
    plt.tight_layout()






.. image-sg:: /auto_examples/sampler/images/sphx_glr_plot_metropolis_hastings_002.png
   :alt: Metropolis Hastings samples
   :srcset: /auto_examples/sampler/images/sphx_glr_plot_metropolis_hastings_002.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  4.704 seconds)


.. _sphx_glr_download_auto_examples_sampler_plot_metropolis_hastings.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example


    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_metropolis_hastings.py <plot_metropolis_hastings.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_metropolis_hastings.ipynb <plot_metropolis_hastings.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_