The standard Bayesian approach to model fitting involves sampling the posterior, usually via a variant of Markov Chain Monte Carlo (MCMC). Though there are many very sophisticated MCMC samplers out there, the most simple algorithm (Metropolis-Hastings) is rather straightforward to code.
Here we'll walk through creating our own Metropolis-Hastings sampler from scratch.
(The following is based on material put together by Adrian Price-Whelan)
As usual, we start with some imports:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# use seaborn plotting defaults
# If this causes an error, you can comment it out.
import seaborn as sns
sns.set()
And load some data:
from fig_code import linear_data_sample
x, y, dy = linear_data_sample()
Walk through all the following steps, filling-in the code along the way.
First plot the data to see what we're looking at (Use a plt.errorbar()
plot with the provided data)
plt.errorbar(x, y, dy, fmt='o');
We're going to fit a line to the data, as we've done through the lecture:
def model(theta, x):
# the `theta` argument is a list of parameter values, e.g., theta = [m, b] for a line
return theta[0] + theta[1] * x
We'll start with the assumption that the data are independent and identically distributed so that the likelihood is simply a product of Gaussians (one big Gaussian). We'll also assume that the uncertainties reported are correct, and that there are no uncertainties on the x
data. We need to define a function that will evaluate the (ln)likelihood of the data, given a particular choice of your model parameters. A good way to structure this function is as follows:
def ln_likelihood(theta, x, y, dy):
# we will pass the parameters (theta) to the model function
# the other arguments are the data
return -0.5 * np.sum(np.log(2 * np.pi * dy ** 2) + ((y - model(theta, x)) / dy) ** 2)
What about priors? Remember your prior only depends on the model parameters, but be careful about what kind of prior you are specifying for each parameter. Do we need to properly normalize the probabilities?
def ln_prior(theta):
# flat prior: log(1) = 0
return 0
Now we can define a function that evaluates the (ln)posterior probability, which is just the sum of the ln prior and ln likelihood:
def ln_posterior(theta, x, y, dy):
return ln_prior(theta) + ln_likelihood(theta, x, y, dy)
Now write a function to actually run a Metropolis-Hastings MCMC sampler. Ford (2005) includes a great step-by-step walkthrough of the Metropolis-Hastings algorithm, and we'll base our code on that. Fill-in the steps mentioned in the comments below:
def run_mcmc(ln_posterior, nsteps, ndim, theta0, stepsize, args=()):
"""
Run a Markov Chain Monte Carlo
Parameters
----------
ln_posterior: callable
our function to compute the posterior
nsteps: int
the number of steps in the chain
theta0: list
the starting guess for theta
stepsize: float
a parameter controlling the size of the random step
e.g. it could be the width of the Gaussian distribution
args: tuple (optional)
additional arguments passed to ln_posterior
"""
# Create the array of size (nsteps, ndims) to hold the chain
# Initialize the first row of this with theta0
chain = np.zeros((nsteps, ndim))
chain[0] = theta0
# Create the array of size nsteps to hold the log-likelihoods for each point
# Initialize the first entry of this with the log likelihood at theta0
log_likes = np.zeros(nsteps)
log_likes[0] = ln_posterior(chain[0], *args)
# Loop for nsteps
for i in range(1, nsteps):
# Randomly draw a new theta from the proposal distribution.
# for example, you can do a normally-distributed step by utilizing
# the np.random.randn() function
theta_new = chain[i - 1] + stepsize * np.random.randn(ndim)
# Calculate the probability for the new state
log_like_new = ln_likelihood(theta_new, *args)
# Compare it to the probability of the old state
# Using the acceptance probability function
# (remember that you've computed the log probability, not the probability!)
log_p_accept = log_like_new - log_likes[i - 1]
# Chose a random number r between 0 and 1 to compare with p_accept
r = np.random.rand()
# If p_accept>1 or p_accept>r, accept the step
# Else, do not accept the step
if log_p_accept > np.log(r):
chain[i] = theta_new
log_likes[i] = log_like_new
else:
chain[i] = chain[i - 1]
log_likes[i] = log_likes[i - 1]
return chain
Now run the MCMC code on the data provided.
chain = run_mcmc(ln_posterior, 10000, 2, [0, 1], 0.1, (x, y, dy))
Plot the position of the walker as a function of step number for each of the parameters. Are the chains converged?
plt.plot(chain[:, 0])
plt.plot(chain[:, 1]);
Make histograms of the samples for each parameter. Should you include all of the samples?
plt.hist(chain[1000:, 0], alpha=0.5)
plt.hist(chain[1000:, 1], alpha=0.5);
It's also sometimes useful to view a two-dimensional histogram:
plt.hist2d(chain[1000:, 0], chain[1000:, 1], bins=20)
plt.xlabel('intercept')
plt.ylabel('slope')
plt.grid(False);
Report to us your constraints on the model parameters -- you have some freedom in interpreting what this means...
theta_best = np.mean(chain[1000:], 0)
print("slope, intercept =", theta_best)
slope, intercept = [-2.74153022 3.1508756 ]