In [ ]:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import diymcmc

We've written a convenience function to grab the test data:

In [ ]:

x,y,sigma_y = diymcmc.get_data1()

Start by plotting the data -- just make an errorbar() plot with the provided data.

In [ ]:

# MAKE PLOT

What model should we use for the data? Why? Create a function to evaluate this model for a given choice of parameters and x data:

In [ ]:

def model(pars, x):
    # the `pars` argument is a list of parameter values, e.g., pars = [m, b] for a line
    pass

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:

In [ ]:

def ln_likelihood(pars, x, y, y_unc):
    # we will pass the parameters (pars) to the model function
    # the other arguments are the data
    pass 

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?

In [ ]:

def ln_prior(pars):
    pass

Now we can define a function that evaluates the (ln)posterior probability, which is just the sum of the ln prior and ln likelihood:

In [ ]:

def ln_posterior(pars, x, y, y_unc):
    return ln_prior(pars) + ln_likelihood(pars, x, y, y_unc)

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

In [ ]:

def run_mcmc(ln_posterior, nsteps, ndim, p0, stepsize, args=()):
    
    # Set up the chain, and initialize it
    # (need an array holding the positions at each step
    
    # Set up an array to hold the probabilities at each step
    
    # Calculate the probability for the first step
    
    # Loop for nsteps
    for i in np.linspace(1,nsteps-1,nsteps-1):
        # Randomly choose new model parameters for the trial state
        
        # Calculate the probability for the new state
        
        # Compare it to the probability of the old state
        # Using the acceptance probability function
        
        # Chose a random number u between 0 and 1 to compare with p_accept
        
        # If p_accept>1 or p_accept>u, accept the step
            # Save the position to the chain
            
            # Save the probability to that array
            
        # Else, do not accept the step
            # Set the position and probability are equal to the last value
            
    # Return the chain and probabilities

Now run the MCMC code on the data provided.

In [ ]:

# run it!

Plot the position of the walker as a function of step number for each of the parameters. Are the chains converged?

In [ ]:

Make histograms of the samples for each parameter. Should you include all of the samples?

In [ ]:

Report to us your constraints on the model parameters -- you have some freedom in interpreting what this means...

In [ ]: