Crash Course

The best way to get started with nautilus is to apply it to a problem. Here, we give a minimal example of using the code to estimate parameter posteriors and the Bayesian evidence \(\mathcal{Z}\). In this example, our “model” has three parameters: \(a\), \(b\) and \(c\). Furthermore, let’s assume flat priors for \(a\) and \(b\) in the range \([-5, +5]\) and a Gaussian prior for \(c\) with \(0\) mean and \(2\) scatter. The priors can be specified as follows.

from scipy.stats import norm
from nautilus import Prior

prior = Prior()
prior.add_parameter('a', dist=(-5, +5))
prior.add_parameter('b', dist=(-5, +5))
prior.add_parameter('c', dist=norm(loc=0, scale=2.0))

The next step is to define the likelihood. For simplicity, let’s assume that all three parameters have zero mean and unit variance but that \(a\) and \(c\) have a correlation coefficient of \(0.9\).

import numpy as np
from scipy.stats import multivariate_normal

def likelihood(param_dict):
    x = np.array([param_dict['a'], param_dict['b'], param_dict['c']])
    return multivariate_normal.logpdf(
        x, mean=np.zeros(3), cov=[[1, 0, 0.90], [0, 1, 0], [0.90, 0, 1]])

Now, we are ready to run the sampler.

from nautilus import Sampler

sampler = Sampler(prior, likelihood, n_live=1000)
sampler.run(verbose=True)

With verbose=True, the sample will give updates on the status. The output will look something like this.

Starting the nautilus sampler...
Please report issues at github.com/johannesulf/nautilus.
Status    | Bounds | Ellipses | Networks | Calls    | f_live | N_eff | log Z
Computing | 7      | 1        | 4        | 15600    | 0.3454 | 8039  | +0.01

The status is either computing the likehihoods (Computing), sampling new points from boundaries (Sampling), or creating new bounds (Bounding). The subsequent entries show the total number of bounds, the number of ellipsoids in the latest bound, the number of neural networks in the latest bound, the total number of likelihood evaluations, the fractional evidence in the live set, the effective sample size and the natural log of the evidence. Once the algorithm is finished, we can plot the posterior using the handy corner Python package.

import corner
import matplotlib.pyplot as plt

points, log_w, log_l = sampler.posterior()
ndim = points.shape[1]
fig, axes = plt.subplots(ndim, ndim, figsize=(3.5, 3.5))
corner.corner(points, weights=np.exp(log_w), bins=20, labels=prior.keys,
              plot_datapoints=False, plot_density=False,
              fill_contours=True, levels=(0.68, 0.95),
              range=np.ones(ndim) * 0.999, fig=fig)

The Bayesian evidence \(\log \mathcal{Z}\) was also estimated during the run.

print('log Z: {:.2f}'.format(sampler.log_z))

Output:

log Z: -6.34