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

.. only:: html

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

        :ref:`Go to the end <sphx_glr_download_dev__auto_dev_JointDistribution.py>`
        to download the full example code.

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

.. _sphx_glr_dev__auto_dev_JointDistribution.py:


Joint Distribution technical details
====================================

    This shows technical aspects of the JointDistribution class.
    These are mostly relevant for developers and advanced users.
    See :doc:`Joint distribution tutorial <../../user/_auto_tutorials/JointDistribution>`
    for a more user-friendly introduction.

.. GENERATED FROM PYTHON SOURCE LINES 12-14

Setup
-----

.. GENERATED FROM PYTHON SOURCE LINES 14-27

.. code-block:: Python


    import sys; sys.path.append("../..")
    import numpy as np
    from cuqi.distribution import Gaussian, Gamma, JointDistribution
    from cuqi.testproblem import Deconvolution1D

    # Model and data
    A, y_obs, _ = Deconvolution1D().get_components()

    # Model dimensions
    n = A.domain_dim
    m = A.range_dim








.. GENERATED FROM PYTHON SOURCE LINES 28-46

Defining hierarchical Bayesian models
-------------------------------------

The joint distribution is general enough to allow us to define hierarchical
Bayesian models.

For example, consider the following extension of the Bayesian model earlier:

.. math::

    \begin{align}
        d &\sim \mathrm{Gamma}(1, 10^{-4}) \\
        l &\sim \mathrm{Gamma}(1,10^{-4}) \\
        \mathbf{x} &\sim \mathcal{N}(\mathbf{0}, d^{-1} \mathbf{I}_n) \\
        \mathbf{y} &\sim \mathcal{N}(\mathbf{A} \mathbf{x}, s^{-1} \mathbf{I}_m)
    \end{align}

We can write this model in CUQIpy as follows:

.. GENERATED FROM PYTHON SOURCE LINES 46-53

.. code-block:: Python


    # Define distribution
    d = Gamma(1, 1e-4)
    l = Gamma(1, 1e-4)
    x = Gaussian(np.zeros(n), lambda d: 1/d)
    y = Gaussian(lambda x: A@x, lambda l: 1/l, geometry=m)








.. GENERATED FROM PYTHON SOURCE LINES 54-55

Define joint distribution p(d,l,x,y)

.. GENERATED FROM PYTHON SOURCE LINES 55-58

.. code-block:: Python

    joint = JointDistribution(d, l, x, y)
    print(joint)





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

 .. code-block:: none

    JointDistribution(
        Equation: 
            p(d,l,x,y) = p(d)p(l)p(x|d)p(y|x,l)
        Densities: 
            d ~ CUQI Gamma.
            l ~ CUQI Gamma.
            x ~ CUQI Gaussian. Conditioning variables ['d'].
            y ~ CUQI Gaussian. Conditioning variables ['x', 'l'].
    )




.. GENERATED FROM PYTHON SOURCE LINES 59-60

Define posterior

.. GENERATED FROM PYTHON SOURCE LINES 60-63

.. code-block:: Python

    posterior = joint(y=y_obs)
    print(posterior)





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

 .. code-block:: none

    JointDistribution(
        Equation: 
            p(d,l,x|y) ∝ p(d)p(l)p(x|d)L(x,l|y)
        Densities: 
            d ~ CUQI Gamma.
            l ~ CUQI Gamma.
            x ~ CUQI Gaussian. Conditioning variables ['d'].
            y ~ CUQI Gaussian Likelihood function. Parameters ['x', 'l'].
    )




.. GENERATED FROM PYTHON SOURCE LINES 64-84

Enabling Gibbs sampling
-----------------------

One of the main purposes of the joint distribution is to be able to sample
using a Gibbs scheme. For this to work we need to be able to define the
distribution of each variable conditioned on the other variables.

That is, we need to define

.. math::

      C_\mathbf{x}&=p(\mathbf{x} \mid \hat{\mathbf{y}}, \hat{d}, \hat{l})\\
      C_d&=p(d \mid \hat{\mathbf{y}}, \hat{\mathbf{x}}, \hat{l})\\
      C_l&=p(l \mid \hat{\mathbf{y}}, \hat{\mathbf{x}}, \hat{d})

Assuming we have some fixed values for :math:`\mathbf{x}`, :math:`d` and :math:`l`,
which we have denoted with the hat symbol. These simply indicate any fixed value.

Then we can simply condition the joint distribution on these values and it will
handle the rest.

.. GENERATED FROM PYTHON SOURCE LINES 84-100

.. code-block:: Python


    # Assume we want to condition on these values
    yh = y_obs
    xh = np.ones(n)
    dh = 1
    lh = 1

    # The conditionals can be computed as follows:
    Cx = joint(y=yh, d=dh, l=lh)
    Cd = joint(y=yh, x=xh, l=lh)
    Cl = joint(y=yh, x=xh, d=dh)

    # We can try inspecting one of these conditional distributions.
    # Notice how conditional distributions changed into a posterior distribution.
    print(Cd)





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

 .. code-block:: none

    Posterior(
        Equation:
             p(d|x) ∝ L(d|x)p(d)
        Densities:
            x ~ CUQI Gaussian Likelihood function. Parameters ['d'].
            d ~ CUQI Gamma.
     )




.. GENERATED FROM PYTHON SOURCE LINES 101-108

Going into the internals of the joint distribution
--------------------------------------------------

A useful "hack" is to have the joint distribution act as if it were a single
parameter density. This is achieved by calling the `._as_stacked()` method.
This returns a new "stacked" joint distribution that the samplers/solvers
can use as if it were any other Density.

.. GENERATED FROM PYTHON SOURCE LINES 108-113

.. code-block:: Python


    posterior_stacked = posterior._as_stacked()

    print(posterior_stacked)





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

 .. code-block:: none

    _StackedJointDistribution(
        Equation: 
            p(d,l,x|y) ∝ p(d)p(l)p(x|d)L(x,l|y)
        Densities: 
            d ~ CUQI Gamma.
            l ~ CUQI Gamma.
            x ~ CUQI Gaussian. Conditioning variables ['d'].
            y ~ CUQI Gaussian Likelihood function. Parameters ['x', 'l'].
    )




.. GENERATED FROM PYTHON SOURCE LINES 114-116

Here the dimension is the sum of the dimensions of the individual variables.
The geometry is assumed as a default geometry matching the dimension.

.. GENERATED FROM PYTHON SOURCE LINES 116-119

.. code-block:: Python

    print(f"Stacked Posterior dimension: {posterior_stacked.dim}")
    print(f"Stacked Posterior geometry: {posterior_stacked.geometry}")





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

 .. code-block:: none

    Stacked Posterior dimension: 130
    Stacked Posterior geometry: _DefaultGeometry1D(130,)




.. GENERATED FROM PYTHON SOURCE LINES 120-122

Finally you can evaluate the log density of the stacked using a single vector
of parameters.

.. GENERATED FROM PYTHON SOURCE LINES 122-126

.. code-block:: Python


    logd = posterior_stacked.logd(np.ones(posterior_stacked.dim))

    print(logd)




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

 .. code-block:: none

    [-370.76796028]





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

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


.. _sphx_glr_download_dev__auto_dev_JointDistribution.py:

.. only:: html

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

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

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

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

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

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: JointDistribution.zip <JointDistribution.zip>`


.. only:: html

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

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