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

.. only:: html

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

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

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

.. _sphx_glr_auto_examples_inference_plot_bayes_inference.py:


Bayesian inference
==================

.. GENERATED FROM PYTHON SOURCE LINES 8-16

This example shows how to perform Bayesian inference given the uniform prior

:math:`p(\mathbf{x})=p(x_1,x_2)=0.01`

where :math:`x_i \in [-5,5], i=1,2`, and likelihood

:math:`L(\mathbf{x}|\mathbf{d})=\exp \left(-\frac{1}{100}\left(x_1-1\right)^2-\left(x_1^2-x_2\right)^2\right)`.


.. GENERATED FROM PYTHON SOURCE LINES 17-29

.. code-block:: default


    import numpy as np

    ndim = 2
    bounds = np.array([[-5,5],[-5,5]])

    def prior(x):
            return 0.01

    def likelihood(x):
            return np.exp(-(x[0]-1)**2/100 - (x[0]**2-x[1])**2)








.. GENERATED FROM PYTHON SOURCE LINES 30-33

To estimate the posterior using grid estimation, we need to import the
:class:`.GridEstimation` class, create an instance, and call the 
:py:meth:`.GridEstimation.run` method.

.. GENERATED FROM PYTHON SOURCE LINES 34-40

.. code-block:: default


    from psimpy.inference.bayes_inference import GridEstimation

    grid_estimator = GridEstimation(ndim, bounds, prior, likelihood)
    posterior, x_ndim = grid_estimator.run(nbins=[50,40])








.. GENERATED FROM PYTHON SOURCE LINES 41-43

The following figure plots the estimated posterior.


.. GENERATED FROM PYTHON SOURCE LINES 44-61

.. code-block:: default


    import matplotlib.pyplot as plt

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

    # mask insignificant values
    posterior = np.where(posterior < 1e-10, np.nan, posterior)

    contour = ax.contour(x_ndim[0], x_ndim[1], np.transpose(posterior), levels=10)
    plt.colorbar(contour, ax=ax)
    ax.set_xlim(bounds[0,0], bounds[0,1])
    ax.set_ylim(bounds[1,0], bounds[1,1])
    ax.set_title('Grid estimation')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    plt.tight_layout()




.. image-sg:: /auto_examples/inference/images/sphx_glr_plot_bayes_inference_001.png
   :alt: Grid estimation
   :srcset: /auto_examples/inference/images/sphx_glr_plot_bayes_inference_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 62-67

To estimate the posterior using Metropolis Hastings estimation, we need to import
the :class:`.MetropolisHastingsEstimation` class, create an instance, and call the 
:py:meth:`.MetropolisHastingsEstimation.run` method. The 
:py:meth:`.MetropolisHastingsEstimation.run` method has a parameter, ``mh_sampler``,
which takes an instance of :class:`.MetropolisHastings` as argument.

.. GENERATED FROM PYTHON SOURCE LINES 68-92

.. code-block:: default


    from psimpy.inference.bayes_inference import MetropolisHastingsEstimation

    mh_estimator = MetropolisHastingsEstimation(ndim, bounds, prior, likelihood)

    # create a mh_sampler
    from psimpy.sampler.metropolis_hastings import MetropolisHastings
    from scipy.stats import multivariate_normal

    init_state = np.array([-4,-4])
    f_sample = multivariate_normal.rvs
    nburn = 100
    nthin = 10
    seed = 1
    kwgs_f_sample = {'random_state': np.random.default_rng(seed)}

    mh_sampler = MetropolisHastings(ndim=ndim, init_state=init_state,
        f_sample=f_sample, bounds=bounds, nburn=nburn, nthin=nthin, seed=seed,
        kwgs_f_sample=kwgs_f_sample)
    
    nsamples = 5000
    mh_samples, mh_accept = mh_estimator.run(nsamples, mh_sampler)









.. GENERATED FROM PYTHON SOURCE LINES 93-96

The following figure plots the samples drawn from the unnormalized posterior,
which can be used to estimate the posterior and its poperties.


.. GENERATED FROM PYTHON SOURCE LINES 97-110

.. code-block:: default

    import matplotlib.pyplot as plt

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

    ax.scatter(mh_samples[:,0], mh_samples[:,1], s=10, c='r', marker='o', alpha=0.1)
    ax.set_xlim(bounds[0,0], bounds[0,1])
    ax.set_ylim(bounds[1,0], bounds[1,1])
    ax.set_title('MH estimation')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    plt.tight_layout()





.. image-sg:: /auto_examples/inference/images/sphx_glr_plot_bayes_inference_002.png
   :alt: MH estimation
   :srcset: /auto_examples/inference/images/sphx_glr_plot_bayes_inference_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 111-116

.. note:: Besides ``prior`` and ``likelihood``, one can also instantiate
    the :class:`.MetropolisHastingsEstimation` class with

     - ``ln_prior`` and ``ln_likelihood``: Natural logarithm of ``prior`` and ``likelihood``.
     - ``ln_pxl``: Natural logarithm of the product of ``prior`` and ``likelihood``.


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

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


.. _sphx_glr_download_auto_examples_inference_plot_bayes_inference.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_bayes_inference.py <plot_bayes_inference.py>`

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

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


.. only:: html

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

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