MCMC Basics

Here is a quick outline of how to use the MCMC tool in henrietta.

import henrietta as hsl

Before we can create a light curve model using an MCMC, we need to download some light curve data to work with.

lc = hsl.download_kepler_lc('Kepler-10', quarter=1)

Then, we need to define our custom BATMAN model, including which parameters we want to fix and which parameters we want to fit.

astropy_model = hsl.setup_transit_model()

Let’s check out what the default model parameters are:

astropy_model

Kepler-10b has a period of 0.8375 days and an impact parameter of 0.3. Let’s fix those values, but let the radius, t0, and a be free parameters, with ranges defined around the anticipated true value:

astropy_model = hsl.setup_transit_model(period = 0.8375,
                                        b = 0.3,
                                        t0 = [0.09,0.12],
                                        radius = [0.005,0.02],
                                        a = [1.0,4.0])

Now that we have our astropy_model object and our light curve data in the form of a lightkurve object, we are ready to model these parameters with a Markov-Chain Monte Carlo.

The mcmc_fit function takes the astropy model, the light curve data, and one additional argument that will determine whether or not output plots will be saved to the local directory. You might want to start with a small number of steps in the MCMC chain (say 100 or 1000) to see how long it takes, and then scale up to longer MCMC runs from there. This mcmc_fit function returns two useful variables:

max_likelihood is a dictionary with keys equal to the names of the free parameters. This dictionary contains 3 values for each free parameter - the maximum likelihood value (detered as the median sampled parameter) and the upper and lower 1-sigma uncertainty parameter values.

samples is an object that contains many different tools for examining the MCMC results. For an in-depth look at the capabilities of this object, the user should consult the emcee handbook: http://dfm.io/emcee/current/

max_likelihood, samples = hsl.mcmc_fit(astropy_model,lc,saveplots=True, nsteps=1000)

If we open up our current directory, we should see some .pdf files that contain summaries of diagnostic plots made by this MCMC. We can also look at the max_likelihood variable now, to see the central 1-sigma ranges for each of our fitted parameters.

max_likelihood

There’s lots to learn about MCMC fitting of transits, but this tool might be a handy way to get started!