This web page describes the latest procedures used for
fitting light curve (LC) submissions to the AXA.
Amatuer LC data is almost always "noisier" than professional data. Often
the level of noise does not permit a solution to be sought for some of the
subtle aspects of LC shape. Starting April 14, 2008 I'll use a procedure
that allows for these limitations without sacrificing the ability to extract
mid-transit time from noisy data.
Another difference between amatuer and professional LC data is the size
of systematic errors. There are two main systematics: 1) temporal trends
(produced by polar axis mis-alignment that cause image rotation, which in
turn produces trends as the star field moves across an imperfect flat field),
and 2) "air mass curvature" (produced differences in star color between the
target star and reference stars). The first effect is usually greater for
amateurs because our polar axis alignment is rarely as good as to professional
observastories, nor are our flat fields as good, usually. The second effect
is usually greater for amateurs because our observing sites are at lower
altitudes, where extinction is greater, and where the effect of star color
differences is greater. Note that each star's extinction coefficient is slightly
different depending on its color.
Because the two main systematic effects on LCs are greater for amateurs
there's a greater need to address them and remove their effect when processing
amateur data. In order for the AXA viewer to see how much systematics may
be affecting the final LC I have decided to include two LC versions with
each observer's data submission: 1) with systematics removed, and 2) without
I use the term "model" reluctantly because my model is primitive compared
with the models used by professionals. The following diagram illustrates
what I'm calling a "model."
Figure 1. LC shape model used in fitting AXA submitted transit
data. Trend and air mass curvature are present. First contact occurs at
t1, etc. Descriptions are in the text.
This "model" includes 7 parameters (plust an offset), listed here:
1) ingress time, t1
2) egress time, t4
3) depth at mid-transit, D [mmag]
4) fraction of transit that's partial, Fp = [ (t2-t1) + (t4-t3) ] / (t4-t1)
5) relative depth at contact 2, F2 (depth at contact 2 divided by depth at mid-transit),
6) temporal trend [mmag/hour], and
7) air mass curvature [mmag/airmass].
Mid-transit is "flattened" for a time interval given by (t3-t2)/4. Including
an "air mass curvature" term requires that air mass versus time be
represented by an equation including terms for powers of (t-to). I use a
4th power expansion with to set to the mid-point of the observations.
The interval t1 to t2 is "partial transit." Same for t3 to t4. The parameter
"Fp" is the fraction of the total transit time that's partial. When a transit
is "central" (impact parameter b = 0) Fp can be used to estimate the size
of the transiting exoplanet (or brown dwarf star) compared with the size
of the star being obscured.
The parameter "F2" is used to specify the amount of this limb darkening
The interval t2 to t3 is not flat due to stellar "limb darkening."
Depth at mid-transit, D [mmag], can be used to estimate the size of the
transiting exoplanet (or brown dwarf).
Least-Squares Fitting Procedure
Many steps are needed to perform a LS fitting of a model to data. Sorry
for the detail in what follows, but some readers may want to know.
Deriving SE versus time
Every measurement has an associated uncertainty. If a region of data share
the same SE (standard error) uncertainty, it is possible to deduce the value
for SE by comparing a measured value with the average of its neighbors.
This assumes that the thing being measured changes slowly compared with
the time spanned by the neighbor values being averaged. Let's clarify this:
the average of neighbors will be a good representation of the true value
being measured provided the second derivative of the thing being measured
is small (I won't belabor this since if you know about second derivatives
you'll know what I'm getting at). I like using the 4 preceding neighbors
and the 4 following neighbors, and I also prefer using their median instead
of average (to be less affected by outliers). The median of 8 neighbors will
itself have an uncertainty, but it will be much smaller than that for the
measurement under consideration. Using the median of 8 neighbors leads to
a SE on the 8 neighbors of ~ SEi ×1.15 × sqrt(1/7), where
SEi is the SE per individual measurement (the 1.15 accounts for use of median
instead of average). When a measurement is compared with this slightly uncertain
version of "truth" the difference is slightly more uncertain than the measurement's
actual SE. The histogram of many such differences will have a spread that
is larger than SEi by 9%. My procedure for calculating SEi versus time during
the observing session is to calculate the standard deviation of the "neighbor
differences" data at each time step, then divide the entire sequence by 1.09.
This description left out one data quality checking step. The "neighbor
differences" data can be used to identify outliers. When the absolute value
of a measurement's neighbor difference exceeds a user-specified rejection
criterion the measurement is rejected! Although the rejection criterion
should depend on the total number of measurements, and should be set to
a multiple of the SEi determined to exist at the timeof the measurement
under question, I prefer to set the criterion to about 2.5 times the typical
SEi for the observing session. This usually leads to a rejection of about
3% of the data (it would be 1.2% if only Gaussian noise were present). There's
one other data quality checking step that is applied to some submitted data.
When check stars (or even just one check star) are included in the submission
it's possible to assess the amount of loss due to cirrus clouds, or dew on
the corrector plate, or PSF spreading outsdie the aperture circle (due to
seeing changes or focus drift), etc. I usually reject all data that has losses
exceeding ~0.05 magnitude.
Solving for Model Parameter Values
Minimizing chi-square is a subset of least squares. Both assume that a
measurement is actually a probability function having width given by SEi.
The goal is to adjust model parameters to maximize the combined probability
that the model "belongs" to the data. I suppose the word "belongs" should
For a single measurement with SEi specified there is just one parameter
that can be adjusted for assessing "belong." The key concept is that every
parameter value under consideration has a probability of "belonging" to
the measurment, and this probability is simply the value of the Gaussian
function at the parameter value's location. With two measured values the
combined probability is the product of each probability. Since the Gaussian
probability function is proportional to the inverse square of the difference,
divided by SEi, the maximum combined probability will be achieved by a model
that also minimizes the sum of squares of "model minus measurement divided
by SEi." Chi-squared analysis involves finding a set of parameter values
that minimize this sum, and it is equivalent to maximizing probability products.
(Probability prodcuts is a Bayesian technique, and it lends itself to the
case when probability functions are not Gaussian.)
There's one more feature available when minimizing chi-squared (or maximizing
probability products). Outside information about likely values for model
parameters can be explicitly made use of by treating them in a way similar
to the way measurements are treated. If a model parameter has a known a
priori probability function it can be included in the analysis. For example,
if a Bayesian approach (product of probabilities) is used then whenever a
model parameter is being evaluated the a priori probability that a value
for it that is being considered can be included in the set of probability
products (most of which come from comparing measurements with model predicted
values). If it is assumed that all probability functions are Gaussian, then
the chi-squared procedure can be used and the model parameter's difference
from it's most likely value can be included in the sum of chi-squares. For
each model parameter the additional term to include in the chi-square sum
is "model parameter value minus most likely value divided by SE for the parameter."
The effect of including this term is to reward model parameter values that
agree with expectations, and punish model parameter values that stray too
far from what is expected. For want of a better term I'll say that my auto-fit
procedure minimizes "Super Chi-Squared" where:
where measured magnitudes Oi and model computed magnitudes Ci are differenced
and divided by measured magnitude SEi, then squared and summed for all measurements;
followed by a set of model terms consisting of model parameter value Pi differenced
with an a priori (most likely) value for the parameter Pm (based on
a history of fitted light curves for the exoplanet in question) and divided
by an a priori uncertainty of this parameter SEp, which are also squared.
Of course there is a potential danger in assigning a priori SE for
model parameter values and punishing solutions that involve parameters that
stray too far from what's expected. When an exoplanet transit has never been
observed there is little information besides some basic physics for bounding
parameter values (e.g., depth can't exceed 100% would be one example). There
are occasions when using a priori SE for model parametsr is appropriate.
Consider the following most common one for exoplanet transits. XO-1 has been
observed on a few occasions with high precision, so the shape and depth and
length of the LC is well established. For example, Fp = 0.26 ±
0.02, F2 = 0.83 ± 0.05, D = 23.8 ± 0.5 mmag (for R-band) and
L = 2.91 ± 0.03 hours. Whenever a new high-precision complete
transit observation is being fitted it is not necessary to makes use of this
a priori information since the new LC data can be sued to solve for
all parameters. However, when a low precision LC data set is to be fitted,
or when only an ingress (or egress) was observed, it is not possible to solve
for Fp and F2, and sometimes it is also not possible to solve for L or D
with an expectation of reasonable results. Should this noisy or partial transit
data set be rejected? No, if we constrain the parameters for whcih we have
confidence and solve for the parameter(s) that may be changing, such as mid-transit
time, we will be making legitimate use of prior knowledge as we try to learn
about something for which we do not have prior knowledge. In the case of
XO-1, for example, this permits us to make use of noisy LC data, or partial
transit data, to solve for mid-transit time - a parameter that can be useful
in a search for transit timing variations (TTV) produced by another exoplanet
in a resonant orbit.
Professional LC data is usually very precise, for several reasons: large
apertures have higher SNR and lower scintillation effects, professional
observatories are at high altitude where extinction is less and image rotation
is non-existent because polar axes are well aligned. Professional LCs therefore
should have small temporal trends, small air mass curvatures and should
never need to have model parameters constrained by a priori information.
We poor amateurs, on the other hand, must work harder to overcome all the
handicaps of not benefitting from the above list. This web page section is
long because I need to explain the many procedure for overcoming amateur
I have decided to present two versions of every amateur LC so that the
severity of sytematics and my success at removing them can be assessed.
One LC version has systematics removed (which is what the professionals
do) and immediately below it is another LC without the systematics removed.
This dual presentation will be done for all submissions after April 15, 2008.
I will use the first of these to illustrate what information is contained
in the two LC versions.
Amateur LC with a chi-square fit that includes a priori values
for model parameters Fp, F2, D and L. Other details are described in the text.
(fp is the same as Fp, f2 is the same as F2.)
Magnitudes from individual images are shown by small orange dots. Groups
of 5 non-overlapping individual measurements ar shown with large orange
symbols. A running average of the individual measurements is shown by a
wavy light blue trace. The thick gray trace is a model fit. The gray dotted
trace during the transit protion shows how the model would plot if the transit
depth were zero.
GJ 436 is a very red star, and no reference stars are available that are equally red (and are bright enough for use as reference). Therefore it is common for LCs of this object to exhibit "air mass curvature" (since the target star and reference stars will have slightly different extinction coefficients). In addition to the air mass curvature this LC plot exhibits a temporal trend (usually due to imperfect polar axis alignment that produces image rotation). As noted near the lower-right corner the solved-for temporal trend coefficient (labeled "Slope") is -0.85 mmag/hour. The air mass curvature is also presented as -68 mmag/airmass (labeled "ExtCurv"). The same information group gives the single-image SE to be 3.30 mmag for a neighbor noise rejection criterion that leaves 98% of the data as accepted. The exposure time is 120 seconds.
At the bottom central is a list of a priori model parameter assumptions.
D = 8.0 ± 1.5 mmag, L = 0.94 ± 0.10, etc. My choice of SE
for these parameters is twice their established SE from previous high-precision
The lower-left corner shows solved-for parameters whose solution was subject
to the a priori constraints. The SE value for each parameter is obtained
by adjusting the parameter value until chi-square increases by 2. Normally
a chi-sqaure increase of 1 is required, but this assumes that either all
parameters are "orthogonal" (do not mimic each others effects) or when a
parameter is adjusted all other parameters are free to find their minimum
chi-square solution. Since neither condition is met, I have adopted the custom,
provisionally, of using a chi-square increase criterion of 2 instead of 1.
The ingress and egree times call for a mid-transit time that is 10.8 ±
3.4 minutes later than the ephemeris used for predicting the transit. (It
is known that the published ephemeris has a period that is too short, and
LCs like this one are helping to establish a new period.)
The lower pael shows air mass. If reference stars had been included in the submission another curve would show "extra losses" due to cirrus clouds, etc.
Note that this LC was made using an 8-inch telescope. This is the smallest
aperture used for observers who have submitted data to the AXA, and I salute
Manuel Mendez for his successful observing.
The next figure shows the same data with the two systematics removed.
Systematics removed version of an amateur LC. Descritpion is given
below. (fp is the same as Fp, f2 is the same as F2.)
The dashed blue trace is an offset LC shape showing when the transit was
expected. The same grouping of 5 non-overlapping data are shown using the
same orange symbols as in the previous plot. These data have SE bars that
are based on the SEi versus UT (not shown) and divided by square-root of
4. The model fit is horizontal outsdie of transit, by definition. The small
red crosses at the bottom are residuals of the individual image measurements
with respect to the model. The information groups are a repeat from the
On the web pages devoted to each exoplanet I will present both graphs.
Professionals probably don't publish the graph with systemtics present because
their systematics are small and it is therefore not necessary to assess
their importance. Amateur observations, however, should have the systematics
illustrated so that anyone considering using the analysis results can be
warned about errors that may remain after attempting to remove the systematics.
Plans for Future Automation
If you think this fitting procedure is labor intensive, you're right! I
do it using a spreadsheet (LCx_template.xls, which you can download and use
in the same way I use it). If Caltech's IPAC/NStED decides to include and
AXA-like archive with my help, then I will work with Caltech to implement
this procedure using a set of programs that chain to each other (a "pipeline").
Then the entire process will be painless, provided users have the discipline
to submit their data in the required format (something only one of you has
done for me). Wish me luck in getting Caltech to take over the AXA. It needs
to be the responsibility of an institution, because I'm an old man!
If Caltech's NStED decides to not take over the AXA then I don't know
what will happen to this site. In a parallel life I do work as a consultant
for unrelated matters, so if anyone has an idea for funding me to automate
the process I'll consider it.
When a well-observed exoplanet is under study for transit timing variations
(TTV) there is no need to solve for Fp and F2, as this would just add to the
uncertainty of the mid-transit time. It is also not necessary to allow length
and depth be free parameters. Length and depth are usually well established
but not with perfect accuracy. Therefore, for TTV studies I assign L and D
their best estiamte SE, and assing Fp and F2 zero error, and proceed with
the Super Chi-Square minimization procedure inspired by Bayesian Estimation
Theory to solve for mid-transit time. SE on mid-transit time is evaluated
by manually adjusting it so that Super Chi-Square increases by ~1.5, then
finding a best set of parameter values for L and D (and the systematic parameters
trend and AMC), and usually Super Chi-Square comes down to ~1.0. If it doesn't,
then I adjust mid-transit time again, and solve for the other parameters.
This process actually goes quite fast, and the resulting SE for mid-transit
time is an objective evaluation that is valid if stochastic uncertainties
dominate systematic errors.
One systematic that can't be dealt with using this procedure is a clock
error that is not corrected before the observations, or at least measured
afterwards and reported in the data submission file. So, please, observers,
be diligent in monitoring your computer clock accuracy.
There will be an illustration of the TTV search in a couple months on
the AXA for an exoplanet that appears to exhibit it.
Return to AXA home page.
This site opened: April 15, 2008. Last Update: July 03, 2008