Open In Colab

Basic usage

In this tutorial we demonstrate basic usage of harmonic, using emcee as the sampler.

Import packages

To begin with we need to import packages needed.

[12]:
import numpy as np
import sys
sys.path.insert(0, '../../..')
import harmonic as hm
sys.path.append("../../../examples")
import utils
import emcee
import scipy.special as sp
import time
import matplotlib.pyplot as plt
import gc

Define Bayesian posterior function

Now we will need to define the log-posterior function of interest.

As a working example for this basic tutorial we consider a log-likelihood given a standard 2-dimensional Gaussian

\[f(x) = -\frac{1}{2}x^{T}\Sigma^{-1}x\]

where for simplicity we have taken the mean \(\mu=0\) and dropped scaling factors, and assume a trivial uniform prior over an infinite interval. Under such conditions the log-posterior is given by

[13]:
def ln_posterior(x, inv_cov):
    """Compute log_e of posterior of n dimensional multivariate Gaussian.

    Args:

        x: Position at which to evaluate posterior.

    Returns:

        double: Value of posterior at x.

    """

    return -np.dot(x,np.dot(inv_cov,x))/2.0

Compute samples using emcee

We then sample from the posterior using an MCMC algorithm. While any MCMC approach can be used we sample using the emcee package.

First we will need to define and initialise some variables.

[14]:
# Define parameters for emcee sampling
ndim = 2                    # Dimensions
nchains = 200                # total number of chains to compute
samples_per_chain = 5000    # number of samples per chain
nburn = 2000                # number of samples to discard as burn in

# Construct a trivial inverse covariance (identity matrix)
inv_cov = np.zeros((ndim,ndim))
diag_cov = np.ones(ndim)
np.fill_diagonal(inv_cov, diag_cov)

Now we need to run the sampler.

[15]:
# initialise random seed
np.random.seed(1)

# Set initial random position and state
pos = np.random.rand(ndim * nchains).reshape((nchains, ndim))
rstate = np.random.get_state()

# Instantiate and execute sampler
sampler = emcee.EnsembleSampler(nchains, ndim, ln_posterior, args=[inv_cov])
(pos, prob, state) = sampler.run_mcmc(pos, samples_per_chain, rstate0=rstate)

# Collect samples into contiguous numpy arrays (discarding burn in)
samples = np.ascontiguousarray(sampler.chain[:,nburn:,:])
lnprob = np.ascontiguousarray(sampler.lnprobability[:,nburn:])

Compute evidence using harmonic

The harmonic package requires only posterior samples. There are no constraints on the type of sampling algorithm used.

Once we have posterior samples to hand, they can be post-processed using harmonic to compute the Bayesian evidence.

Collating samples into harmonic

We first configure the chains into a harmonic-friendly shape, which we do as follows.

[16]:
# Instantiate harmonic's chains class
chains = hm.Chains(ndim)
chains.add_chains_3d(samples, lnprob)

Since we will subsequently learn the target distribution \(\varphi\) we split the samples into training and inference sets (we often use the common machine learning terminology test for the inference data-set).

[17]:
# Split the chains into the ones which will be used to train and peform inference
chains_train, chains_infer = hm.utils.split_data(chains, training_proportion=0.5)

Train the machine learning model

For simplicity we will manually select hyper-parameters and train a machine learning model using chains_train. See other tutorials that use cross-validation to select hyper-parameters and different models.

Simply select the model we wish to adopt and fit the model.

[18]:
domains = [np.array([1E-1,1E1])]  # hyper-sphere bounding domain

# Select model
model = hm.model.HyperSphere(ndim, domains)

# Train model
fit_success = model.fit(chains_train.samples, chains_train.ln_posterior)

Compute the Bayesian evidence

Finally we simply compute the learnt harmonic mean estimator as follows.

[19]:
# Instantiate harmonic's evidence class
ev = hm.Evidence(chains_infer.nchains, model)

# Pass the evidence class the inference chains and compute the log of the evidence!
ev.add_chains(chains_infer)
evidence, evidence_std = ev.compute_evidence()

Results

Let’s check the evidence value computed and also plot the posterior.

Numerical integration

As this is a standard 2-dimensional Gaussian the evidence is analytic and given by

[20]:
def ln_analytic_evidence(ndim, cov):
    """Compute analytic evidence for nD Gaussian.

    Args:

        ndim: Dimension of Gaussian.

        cov: Covariance matrix.

    Returns:

        double: Analytic evidence.

    """

    ln_norm_lik = 0.5*ndim*np.log(2*np.pi) + 0.5*np.log(np.linalg.det(cov))
    return ln_norm_lik

evidence_analytic = np.exp(ln_analytic_evidence(ndim, inv_cov))

Let’s compare the value computed by harmonic and by numerical integration.

[21]:
print('evidence (harmonic) = {} +/- {}'.format(evidence, evidence_std))
print('evidence (analytic) = {}'.format(evidence_analytic))
print('nsigma = {}'.format(np.abs(evidence - evidence_analytic) / evidence_std))
evidence (harmonic) = 6.284739837844395 +/- 0.01510410316980982
evidence (analytic) = 6.283185307179585
nsigma = 0.10292108358453812

As expected, the evidence computed by harmonic is close to that computed by numerical integration.

Posterior triangle plot

Out of interest let’s also plot slices of the posterior using these samples to see what we’re working with!

[22]:
utils.plot_getdist(samples.reshape((-1, ndim)))
../_images/tutorials_basic_usage_27_0.png