runtransiterModel.py
Using runtransiterModel.py to fit data using a Levenburg-Marquardt least squares algorithm.
Tutorial Contents
So the first thing we want to do is generate fake data using the transiterFit.fake_data function. First, we'll import the necessary packages. Currently you must be in the OSCAAR-extras/examples directory for all this to work and the osccarForTransiter is hard-coded with some of the parameters. We're working on changing this around for ease of access.
from matplotlib import pyplot as plt
import oscaar
from scipy import optimize
from numpy.random import shuffle
import transiterFit
from oscaar.extras.knownSystemParameters import returnSystemParams
Next thing we want to do is to set the planetary parameters. For more information of each one see the
tutorial for modelLightCurve.py
fk_RpRs=0.122
fk_aRs=12.7
fk_per=1.58
fk_inc=89.5
fk_t0=2454344.30867
fk_gam1=0.23
fk_gam2=0.30
fk_ecc=0.0
fk_argper=0.0
stddev=0.0005 #Standard Deviation of data.
The last parameter, stddev, is used to set the standard deviation of the data. The tool uses a gaussian distribution to randomize what would be perfect data. Typical amateur observatories have values of stddev less than 0.001 whereas the Kepler Space Telescope is around 0.0001.
The next step is to create the data using the transiterFit.fake_data function which uses the data we just set doing something like,
timeObs,NormFlux=transiterFit.fake_data(stddev,fk_RpRs,fk_aRs,fk_per,fk_inc,fk_t0,fk_gam1,fk_gam2,fk_ecc,fk_argper)
flux_error=stddev*np.ones(np.size(timeObs))
we create times based on the mid-transit time, and a set exposure time of 45 seconds (typical for exoplanet observations). Additionally, we use the stddev as the uncertainty in each of the data points.
Now we wanna fit our data to a model and determine the planetary parameters. The first step is to set initial guesses for the exoplanet. When not comparing to fake data you will want to set this to the archival results for the observed exoplanet, which you will have the option to pull from the exoplanet website in the near future. It is assumed that orbital parameters such as eccentricity and the argument of Pericenter (periapsis?) are known through other observational methods such as radial velocity. We'll set initial guesses slightly different than the one's set for the fake planet.
# Initial Guesses for Planetary Parameters
RpOverRs = 0.117
aOverRs = 12.9
inclination = 89.5
epoch = 2454344.30767
gamma1 = 0.23
gamma2 = 0.30
period = 1.58
eccentricity = 0.0
longPericenter = 0.0
Now we're all ready to perform the fit! We have data, and initial guesses so all we have to do is pass them through the transiterFit.run_LMfit function like so,
#Run the intial fit with input parameters stated above.
fit,success = transiterFit.run_LMfit(timeObs,NormFlux,flux_error,RpOverRs,aOverRs,inclination,epoch,period,
gamma1,gamma2,eccentricity,longPericenter,
fitLimbDark='quadratic',
plotting=True)
The function fits to the data using the scipy fitting tool, optimize.curve_fit which uses a Levenburg-Marquardt least squares algorithm. Additionally, there is the keyword argument, 'fitLimbDark', which is set as either, False, 'linear', or 'quadratic' and sets the type of limb-darkening law we are using in the fit. If it's set to False we assume uniform stellar brightness. For data w/ low signal to noise it may be best to exclude the limb-darkening parameters from the fit since generally are poorly constrained.
Fitting is not always straight-forward and this is just one way to approach the fit. The LM method typically finds a local minimum so care should be taken to the input parameters for the fit. Since we usually know the planetary/orbital parameters from other observations this shouldn't be a problem. To double check we visually check the output to make sure it converges. You should get something like this,
Looks good! You'll also see in the command terminal the fit results which might be something like this,
Results from the inital fit w/ uncertainties based on the sq. root of the covariance matrix
Rp/Rs 0.122473959378 +/- 0.00158263245467
a/Rs 12.8087441716 +/- 0.335239157429
inc 89.9825837636 +/- 0.9284553196
Mid-Tran Time 2454344.30866 +/- 0.00013559865912
Gamma 1 0.249322920785 +/- 0.0694516408044
Gamma 2 0.441650826105 +/- 0.176312922825
which estimates the uncertainty by calculating the square root of the covariance matrix created in the fitting tool. It's at this point you should check to make sure everything looks reasonable.
To best estimate the accuracy of the model parameters we use the transiterFit.run_MCfit function. This function estimates the uncertainty in each model parameter by using a bootstrap Monte Carlo method where we create multiple data sets based on the scrambled residuals of the initial fit and adding back the model fit. Alternatively, we could use the Prayer-Bead method but the benefit with this method is we generate an arbitrary amount of random datasets. We run the MC error estimation by,
n_iter = 1000 #Number of random datasets to fit to.
transiterFit.run_MCfit(n_iter,timeObs,NormFlux,flux_error,fit,success,plotting=True)
It may take a bit as it's fitting 1000 different datasets but when it's done you should get a plot that looks something like,
Awesome! What we see is the fake data with associated uncertainty, the initial fit which is plotted as the thick black line, and just for fun we can see all the different fit results for each randomly scrambled dataset. In addition we'll get the results shown in the command terminal which show up as something like,
Results from random MC fit . . . . .
Planetary to Stellar Radius: 0.122488498263 +/- 0.000397513327873
Semi-major Axis to Stellar Radius: 12.7930148874 +/- 0.0702191639686
Inclination of Orbit: 89.8848653841 +/- 0.249828874008
Mid-Transit Time [JD]: 2454344.30866 +/- 4.15759531477e-05
Gamma 1: 0.249369710137 +/- 0.0257376041574
Gamma 2: 0.438208135565 +/- 0.0578404395176
The whole script can be executed from the command line after initial parameters are set by running
$python runtransiterModel.py