.. _sphx_glr_examples_plot_shells.py:


===============
Gaussian Shells
===============

Toy likelihood model for stress testing multiple-ellipsoid method.

The problem is:

.. math:: 
   \mathcal{L}(\theta) = \mathrm{circ}(\theta; c_1, r_1, w_1) +
                         \mathrm{circ}(\theta; c_2, r_2, w_2)

where
 
.. math::
   \mathrm{circ}(\theta; c, r, w) = \frac{1}{\sqrt{2 \pi w^2}}
       \exp \left[ - \frac{(|\theta - c| - r)^2}{2 w^2} \right]




.. code-block:: python


    import math
    import time
    from collections import OrderedDict

    import numpy as np
    from numpy.random import RandomState
    from matplotlib import pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D

    import nestle

    rstate = RandomState(0)







In the following block, we define the problem. We use r = 2 and w = 0.1,
meaning that the  gaussian is quite narrow compared to the size of the
sphere.



.. code-block:: python


    r = 2.
    w = 0.1

    const = math.log(1. / math.sqrt(2. * math.pi * w**2))

    def logcirc(theta, c):
        d = np.sqrt(np.sum((theta - c)**2, axis=-1))  # |theta - c|
        return const - (d - r)**2 / (2. * w**2)

    def loglike(theta, c1, c2):
        return np.logaddexp(logcirc(theta, c1), logcirc(theta, c2))

    def prior_transform(x):
        """Defines a flat prior between -6 and 6 in all dimensions."""
        return 12. * x - 6.








Visualize
---------

It helps to visualize the surface in two dimensions. Here, we plot the
likelihood evaluated on a fine grid and the sample points from nested
sampling.



.. code-block:: python


    # likelihood surface in 2-d
    xx, yy = np.meshgrid(np.linspace(-6., 6., 200),
                         np.linspace(-6., 6., 200))
    c1 = np.array([-3.5, 0.])
    c2 = np.array([3.5, 0.])
    Z = np.exp(loglike(np.dstack((xx, yy)), c1, c2))

    # nested sampling result
    c1 = np.array([-3.5, 0.])
    c2 = np.array([3.5, 0.])
    f = lambda theta: loglike(theta, c1, c2)
    res = nestle.sample(f, prior_transform, 2, method='multi', npoints=1000,
                        rstate=rstate)

    fig = plt.figure(figsize=(14., 6.))
    ax = fig.add_subplot(121, projection='3d')
    ax.plot_surface(xx, yy, Z, rstride=1, cstride=1, linewidth=0, cmap='coolwarm')
    ax.set_xlim(-6., 6.)
    ax.set_ylim(-6., 6.)
    ax.set_zlim(0., 4.)
    ax.set_zlabel('L')
    ax.set_title('Likelihood evaluated on fine grid')

    ax = fig.add_subplot(122, projection='3d')
    ax.scatter(res.samples[:,0], res.samples[:, 1], np.exp(res.logl),
               marker='.', c=np.exp(res.logl), linewidths=(0.,), cmap='coolwarm')
    ax.set_xlim(-6., 6.)
    ax.set_ylim(-6., 6.)
    ax.set_zlim(0., 4.)
    ax.set_zlabel('L')
    ax.set_title('Nested sampling points');





.. image:: /examples/images/sphx_glr_plot_shells_001.png
    :align: center




Scaling with dimension
----------------------

Here, we demonstrate how the algorithm scales with dimension and compare
the total evidence to the analytic answer.



.. code-block:: python


    npoints = 1000

    def run(ndim):
        """Convenience function for running in any dimension"""

        c1 = np.zeros(ndim)
        c1[0] = -3.5
        c2 = np.zeros(ndim)
        c2[0] = 3.5
        f = lambda theta: loglike(theta, c1, c2)
    
        return nestle.sample(f, prior_transform, ndim, method='multi',
                             npoints=npoints, rstate=rstate)

    # Run over dimensions and save time for each run.
    results = OrderedDict()
    for ndim in [2, 5, 10, 20]:
        t0 = time.time()
        results[ndim] = run(ndim)
        results[ndim].time = time.time() - t0

    analytic_logz = {2: -1.75,
                     5: -5.67,
                     10: -14.59,
                     20: -36.09}

    print("D  analytic    logz  logzerr   nlike  eff(%)   time")
    for ndim, res in results.items():
        eff = 100. * res.niter/(res.ncall - npoints)
        print("{:2d}   {:6.2f}  {:6.2f}     {:4.2f}  {:6d}   {:5.2f}  {:6.2f}"
              .format(ndim, analytic_logz[ndim], res.logz, res.logzerr,
                      res.ncall, eff, res.time))




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

 Out::

    D  analytic    logz  logzerr   nlike  eff(%)   time
     2    -1.75   -1.80     0.05   11912   37.55    1.64
     5    -5.67   -5.71     0.08   23413   35.74    5.11
    10   -14.59  -14.52     0.12   50193   34.19   10.97
    20   -36.09  -35.74     0.19  193214   19.78   40.58


**Total running time of the script:** ( 1 minutes  1.536 seconds)



.. container:: sphx-glr-footer


  .. container:: sphx-glr-download

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



  .. container:: sphx-glr-download

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

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

    `Generated by Sphinx-Gallery <http://sphinx-gallery.readthedocs.io>`_