In this recipe we will demonstrate how to fit a simple model, namely a line, to some data. Although this example is simple, it illustrates what is the proper way of fitting our models to data and infering the parameters of the models.
Let us first import the main packages that we will use:
# show plots inline in the notebook
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, Math
import zeus
In order to create our synthetic data we need to construct a generative probabilistic model.
We start by defining the straight line model and also setting the true values of the model parameters:
# define the model function
def straight_line(x, m, c):
''' A straight line model: y = m*x + c '''
return m*x + c
# set the true values of the model parameters for creating the data
m_true = 3.5 # gradient of the line
c_true = 1.2 # y-intercept of the line
# Set the x-coordinates of the data points
M = 70 # Number of data points
x = np.sort(10.0 * np.random.rand(M)) # their x-coordinates
We are now ready to generate the synthetic data. To this end, we evaluate the model function at the true values of m (slope) and c (y-intercept) and we add some random Gaussian noise of known amplitude sigma.
# create the data - the model plus Gaussian noise
sigma = 3.0 # standard deviation of the noise
data = straight_line(x, m_true, c_true) + sigma * np.random.randn(M)
We can also plot the generative model and the data:
plt.figure(figsize=(9,6))
plt.errorbar(x, data, yerr=sigma, fmt="o", label='data')
plt.plot(x, straight_line(x, m_true, c_true), '-', lw=2, label='model')
plt.legend()
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.show()
The first step to solve a problem is generally to write down the prior and likelihood functions. An important benefit of MCMC is that none of these probability densities need to be normalised.
Here we'll start with the natural logarithm of the prior probability:
def logprior(theta):
''' The natural logarithm of the prior probability. '''
lp = 0.
# unpack the model parameters from the tuple
m, c = theta
# uniform prior on c
cmin = -10. # lower range of prior
cmax = 10. # upper range of prior
# set prior to 1 (log prior to 0) if in the range and zero (-inf) outside the range
lp = 0. if cmin < c < cmax else -np.inf
# Gaussian prior on m
mmu = 3. # mean of the Gaussian prior
msigma = 10. # standard deviation of the Gaussian prior
lp -= 0.5*((m - mmu)/msigma)**2
return lp
We assume that the likelihood is Gaussian (Normal):
def loglike(theta, data, sigma, x):
'''The natural logarithm of the likelihood.'''
# unpack the model parameters
m, c = theta
# evaluate the model
md = straight_line(x, m, c)
# return the log likelihood
return -0.5 * np.sum(((md - data)/sigma)**2)
The log posterior is just the sum of the log prior and the log likelihood probability density functions:
def logpost(theta, data, sigma, x):
'''The natural logarithm of the posterior.'''
return logprior(theta) + loglike(theta, data, sigma, x)
We initialize and run zeus to sample from the posterior distribution. Thin only takes a few lines of code.
ndim = 2 # Number of parameters/dimensions (e.g. m and c)
nwalkers = 10 # Number of walkers to use. It should be at least twice the number of dimensions.
nsteps = 2000 # Number of steps/iterations.
start = 0.01 * np.random.randn(nwalkers, ndim) # Initial positions of the walkers.
sampler = zeus.EnsembleSampler(nwalkers, ndim, logpost, args=[data, sigma, x]) # Initialise the sampler
sampler.run_mcmc(start, nsteps) # Run sampling
sampler.summary # Print summary diagnostics
Lets plot the chains. We can see that the burn-in phase is very brief.
plt.figure(figsize=(16,1.5*ndim))
for n in range(ndim):
plt.subplot2grid((ndim, 1), (n, 0))
plt.plot(sampler.get_chain()[:,:,0], alpha=0.5)
plt.tight_layout()
plt.show()
We discard the first half of the chain elements, thin the samples by a factor of 10, and flatten the resulted chain. We then proceed to plot the marginal posterior distributions:
# flatten the chains, thin them by a factor of 10, and remove the burn-in (first half of the chain)
chain = sampler.get_chain(flat=True, discard=nsteps//2, thin=10)
# plot marginal posterior distributions
fig, axes = zeus.cornerplot(chain, labels=['m', 'c'], truth=[m_true, c_true]);