Monitoring the flux of a star might reveal the presence of a exoplanets when the planet "transits" across the front of the start as seen from our perspective. It causes a reduction of the apparent brightness of the star, and conducting large monitoring programs for the transit method has become one of the most successful methods to detect exoplanets.
We will use a simple model to describe the transit signal: a symmetric trapezoid. While it is geometrically simple, it is a non-linear model, that means you canot write down a design matrix for it So we need a more general minimizer to fit any given data. By the end of this notebook you will have all you need to fit transit data.
Write a function read_data that takes the object numer (e.g. 7016.01) as an argument and returns two arrays: time and flux, read from the data file ('data/7016.01.txt').
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
Write a function plot_data that takes the object number as an argument and plots time vs.
flux of the data returned by read_data.
In [1]:
Write a function trapezoid implementing the trapezoid model. It should take four parameters (delta, T, tau, t0) describing depth, duration, ingress duration, and center time) and return a trapezoid as a function of time.
In [ ]:
Make four different plots that show how the trapezoid shape changes when you vary each parameter independently (maybe 10 examples per plot).
In [ ]:
Eyeball the plot of the actual 7016.01 transit signal, and try to guess what parameters might fit the data best. Overplot the model on top of the data, and make a plot of the residuals (data - model) in a subplot. Your first guess doesn't have to be spot-on, the residuals should tell you where you're off the most.
In [ ]:
Use your best guess for the parameters to write a function plot_fit that takes the object number and a parameter vector, and then makes the two plots: one with data and the model, and another one with the residuals (i.e. data minus model).
In [ ]:
Define the function fit_trapezoid that uses scipy.optimize.minimize to find the best-fit parameters for the 7016.01 data set, and display these results using plot_fit.
The function you want to minimize to find the "best-fit" parameters is the one that maximizes the likelihood. Since we're dealing with (assumed) Gaussian errors, this is the usual quadratic error form
$$ \chi^2(T,t_0,\tau,d;\ \lbrace (t_i,f_i) \rbrace) = \sum_i^N \frac{\left(f_i - \mathtt{trapezoid}(t_i;\ T,t_0,\tau,d)\right)^2}{\sigma_i^2}, $$there $t_i$ and $f_i$ are the individual measurements of time and flux.
Note: Minimizing $\chi^2$ means altering the parameters at fixed data.
In [ ]: