samsam documentation
====================

The samsam package provides two samplers:

- :doc:`_autosummary/samsam.sam` : a Scaled Adaptive Metropolis algorithm (see [1]_, [2]_, [3]_), to robustly obtain samples from a target distribution,
- :doc:`_autosummary/samsam.covis` : a COVariance Importance Sampling algorithm, to efficiently compute the model evidence (or other integrals).

It additionally includes tools (:doc:`_autosummary/samsam.acf`)
to assess the convergence of the sam sampler (ACF, IAT),
and a few commonly used prior distributions (:doc:`_autosummary/samsam.logprior`).

Installation
------------

Using conda
~~~~~~~~~~~

The samsam package can be installed using conda with the following command:

``conda install -c conda-forge samsam``

Using pip
~~~~~~~~~

It can also be installed using pip with:

``pip install samsam``

Example
-------

Let us first define a simple log-probability function:

.. ipython::

   In [1]: import numpy as np
      ...: import matplotlib.pyplot as plt
      ...: from samsam import sam, covis, acf
      ...: from corner import corner
      ...: np.random.seed(0)
      ...:
      ...: def logprob(x):
      ...:     return(-0.5*(np.sum(x**2) + x.size*np.log(2*np.pi)))

Then we run sam to sample this distribution:

.. ipython::

   In [2]: ndim = 10
      ...: nsamples = 100000
      ...: x0 = np.random.normal(0, 100, ndim)
      ...:
      ...: samples, sam_diagnos = sam(x0, logprob, nsamples=nsamples, print_level=0)
      ...: samples = samples[nsamples//4:]

Let us check that sam converged correctly using the ACF/IAT:

.. ipython::

   @savefig acf.png
   In [3]: R = acf.acf(samples)
      ...: tau = np.arange(samples.shape[0])
      ...:
      ...: plt.figure()
      ...: plt.plot(tau[1:], R[1:])
      ...: plt.xscale('log')
      ...: plt.xlim(1, samples.shape[0])
      ...: plt.xlabel('lag')
      ...: plt.ylabel('ACF')
      ...:
      ...: iat = acf.iat(R=R)
      ...: print('IAT:', iat.max())
      ...: print('Effective number of samples:', samples.shape[0]/iat.max())

Now we plot the corner plot of the parameter samples:

.. ipython::

   @savefig corner.png
   In [4]: corner(samples);

Finally we run covis to compute the log-evidence of the model:

.. ipython::

   In [5]: _, _, covis_diagnos = covis(sam_diagnos['mu'], sam_diagnos['cov'], logprob, nsamples=1000, print_level=0)
      ...: # Should be close to 0 since the logprob is correctly normalized
      ...: print('Log-evidence:', covis_diagnos['logevidence'])

API Reference
-------------

.. autosummary::
   :toctree: _autosummary
   :template: autosummary/custom_module.rst
   :recursive:

   samsam
   samsam.acf
   samsam.logprior

References
----------

.. [1] `Haario et al., "An adaptive metropolis algorithm", 2001 <http://dx.doi.org/10.2307/3318737>`_.
.. [2] `Andrieu & Thoms, "A tutorial on adaptive MCMC", 2008, <https://doi.org/10.1007/s11222-008-9110-y>`_.
.. [3] `Delisle et al., "The HARPS search for southern extra-solar planets. XLIII. A compact system of four super-Earth planets orbiting HD 215152", 2018 <http://adsabs.harvard.edu/abs/2018A\%26A...614A.133D>`_.