Fitting transit times#

import exoplanet

exoplanet.utils.docs_setup()
print(f"exoplanet.__version__ = '{exoplanet.__version__}'")
exoplanet.__version__ = '0.5.4.dev5+g6f498a8'

Fitting for or marginalizing over the transit times or transit timing variations (TTVs) can be useful for several reasons, and it is a compelling use case for exoplanet becuase the number of parameters in the model increases significantly because there will be a new parameter for each transit. The performance of the NUTS sampler used by exoplanet scales well with the number of parameters, so a TTV model should be substantially faster to run to convergence with exoplanet than with other tools. There are a few definitions and subtleties that should be considered before jumping in.

In this tutorial, we will be using a “descriptive” model orbits.TTVOrbit to fit the light curve where the underlying motion is still Keplerian, but the time coordinate is warped to make t0 a function of time. All of the other orbital elements besides t0 are shared across all orbits, but the t0 for each transit will be a parameter. This means that other variations (like transit duration variations) are not currently supported, but it would be possible to include more general effects. exoplanet also supports photodynamics modeling using the orbits.ReboundOrbit for more detailed analysis, but that is a topic for a future tutorial.

It is also important to note that “transit time” within exoplanet (and most other transit fitting software) is defined as the time of conjunction (called t0 in the code): the time when the true anomaly is \(\pi/2 - \omega\). Section 18 of the EXOFASTv2 paper includes an excellent discussion of some of the commonly used definitions of “transit time” in the literature.

Finally, there is a subtlety in the definition of the “period” of an orbit with TTVs. Two possible definitions are: (1) the average time between transits, or (2) the slope of a least squares fit to the transit times as a function of transit number. In exoplanet, we use the latter definition and call this parameter the ttv_period to distinguish it from the period of the underlying Keplerian motion which sets the shape and duration of the transit. By default, these two periods are constrained to be equal, but it can be useful to fit for both parameters since the shape of the transit might not be perfectly described by the same period. That being said, if you fit for both periods, make sure that you constrain ttv_period and period to be similar or things can get a bit ugly.

To get started, let’s generate some simulated transit times. We’ll use the orbits.ttv.compute_expected_transit_times function to get the expected transit times for a linear ephemeris within some observation baseline:

import numpy as np
import matplotlib.pyplot as plt

import exoplanet as xo

np.random.seed(3948)
true_periods = np.random.uniform(8, 12, 2)
true_t0s = true_periods * np.random.rand(2)
t = np.arange(0, 80, 0.01)
texp = 0.01
yerr = 5e-4

# Compute the transit times for a linear ephemeris
true_transit_times = xo.orbits.ttv.compute_expected_transit_times(
    t.min(), t.max(), true_periods, true_t0s
)

# Simulate transit timing variations using a simple model
true_ttvs = [
    (0.05 - (i % 2) * 0.1) * np.sin(2 * np.pi * tt / 23.7)
    for i, tt in enumerate(true_transit_times)
]
true_transit_times = [tt + v for tt, v in zip(true_transit_times, true_ttvs)]

# Plot the true TTV model
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
ax1.plot(true_transit_times[0], true_ttvs[0], ".k")
ax1.axhline(0, color="k", lw=0.5)
ax1.set_ylim(np.max(np.abs(ax1.get_ylim())) * np.array([-1, 1]))
ax1.set_ylabel("$O-C$ [days]")

ax2.plot(true_transit_times[1], true_ttvs[1], ".k")
ax2.axhline(0, color="k", lw=0.5)
ax2.set_ylim(np.max(np.abs(ax2.get_ylim())) * np.array([-1, 1]))
ax2.set_ylabel("$O-C$ [days]")

ax2.set_xlabel("transit time [days]")
_ = ax1.set_title("true TTVs")
../../_images/c66c0d372651ce048428c720cc47160e99e821bb658e795d0e6b140fd85f3126.png

Now, like in the Transit fitting tutorial, we’ll set up the the model using PyMC3 and exoplanet, and then simulate a data set from that model.

import pymc3 as pm
import pymc3_ext as pmx
import aesara_theano_fallback.tensor as tt

np.random.seed(9485023)

with pm.Model() as model:

    # This part of the model is similar to the model in the `transit` tutorial
    mean = pm.Normal("mean", mu=0.0, sd=1.0)
    u = xo.QuadLimbDark("u", testval=np.array([0.3, 0.2]))
    logr = pm.Uniform(
        "logr",
        lower=np.log(0.01),
        upper=np.log(0.1),
        shape=2,
        testval=np.log([0.04, 0.06]),
    )
    r = pm.Deterministic("r", tt.exp(logr))
    b = xo.ImpactParameter(
        "b", ror=r, shape=2, testval=0.5 * np.random.rand(2)
    )

    # Now we have a parameter for each transit time for each planet:
    transit_times = []
    for i in range(2):
        transit_times.append(
            pm.Normal(
                "tts_{0}".format(i),
                mu=true_transit_times[i],
                sd=1.0,
                shape=len(true_transit_times[i]),
            )
        )

    # Set up an orbit for the planets
    orbit = xo.orbits.TTVOrbit(b=b, transit_times=transit_times)

    # It will be useful later to track some parameters of the orbit
    pm.Deterministic("t0", orbit.t0)
    pm.Deterministic("period", orbit.period)
    for i in range(2):
        pm.Deterministic("ttvs_{0}".format(i), orbit.ttvs[i])

    # The rest of this block follows the transit fitting tutorial
    light_curves = xo.LimbDarkLightCurve(u).get_light_curve(
        orbit=orbit, r=r, t=t, texp=texp
    )
    light_curve = pm.math.sum(light_curves, axis=-1) + mean
    pm.Deterministic("light_curves", light_curves)

    # ******************************************************************* #
    # On the folowing lines, we simulate the dataset that we will fit     #
    #                                                                     #
    # NOTE: if you are fitting real data, you shouldn't include this line #
    #       because you already have data!                                #
    # ******************************************************************* #
    y = pmx.eval_in_model(light_curve)
    y += yerr * np.random.randn(len(y))
    # ******************************************************************* #
    # End of fake data creation; you want to include the following lines  #
    # ******************************************************************* #

    pm.Normal("obs", mu=light_curve, sd=yerr, observed=y)

    map_soln = model.test_point
    map_soln = pmx.optimize(start=map_soln, vars=transit_times)
    map_soln = pmx.optimize(start=map_soln, vars=[r, b])
    map_soln = pmx.optimize(start=map_soln, vars=transit_times)
    map_soln = pmx.optimize(start=map_soln)
optimizing logp for variables: [tts_1, tts_0]
100.00% [93/93 00:00<00:00 logp = 4.946e+04]

message: Desired error not necessarily achieved due to precision loss.
logp: 49454.87884962603 -> 49461.281726403526
optimizing logp for variables: [b, logr]
100.00% [17/17 00:00<00:00 logp = 4.946e+04]

message: Optimization terminated successfully.
logp: 49461.281726403526 -> 49463.56218609475
optimizing logp for variables: [tts_1, tts_0]
100.00% [104/104 00:00<00:00 logp = 4.946e+04]

message: Desired error not necessarily achieved due to precision loss.
logp: 49463.56218609475 -> 49463.62677318899
optimizing logp for variables: [tts_1, tts_0, b, logr, u, mean]
100.00% [84/84 00:00<00:00 logp = 4.946e+04]

message: Optimization terminated successfully.
logp: 49463.62677318899 -> 49464.0047187812

Here’s our simulated light curve and the initial model:

plt.plot(t, y, ".k", ms=4, label="data")
for i, l in enumerate("bc"):
    plt.plot(
        t, map_soln["light_curves"][:, i], lw=1, label="planet {0}".format(l)
    )
plt.xlim(t.min(), t.max())
plt.ylabel("relative flux")
plt.xlabel("time [days]")
plt.legend(fontsize=10)
_ = plt.title("map model")
../../_images/99d8e7f7c7e7fc88e0b1ef859bec9891a7affad39616ec31ec8c50a3d1111302.png

This looks similar to the light curve from the Transit fitting tutorial, but if we try plotting the folded transits, we can see that something isn’t right: these transits look pretty smeared out!

for n, letter in enumerate("bc"):
    plt.figure()

    # Get the posterior median orbital parameters
    p = map_soln["period"][n]
    t0 = map_soln["t0"][n]

    # Compute the median of posterior estimate of the contribution from
    # the other planet. Then we can remove this from the data to plot
    # just the planet we care about.
    other = map_soln["light_curves"][:, (n + 1) % 2]

    # Plot the folded data
    x_fold = (t - t0 + 0.5 * p) % p - 0.5 * p
    plt.errorbar(
        x_fold, y - other, yerr=yerr, fmt=".k", label="data", zorder=-1000
    )

    plt.legend(fontsize=10, loc=4)
    plt.xlim(-0.5 * p, 0.5 * p)
    plt.xlabel("time since transit [days]")
    plt.ylabel("relative flux")
    plt.title("planet {0}".format(letter))
    plt.xlim(-0.3, 0.3)
../../_images/4123354fb4a05446e1802a45b34b4fc99af983f57a64f81522ade3237c6d8270.png ../../_images/eea2340293120137caa164e0d4a52ea98e91551d7b47565ca36528597ac5e87d.png

Instead, we can correct for the transit times by subtracting these estimates and plot that instead:

with model:
    t_warp = pmx.eval_in_model(orbit._warp_times(t), map_soln)

for n, letter in enumerate("bc"):
    plt.figure()

    p = map_soln["period"][n]
    other = map_soln["light_curves"][:, (n + 1) % 2]

    # NOTE: 't0' has already been subtracted!
    x_fold = (t_warp[:, n] + 0.5 * p) % p - 0.5 * p
    plt.errorbar(
        x_fold, y - other, yerr=yerr, fmt=".k", label="data", zorder=-1000
    )

    plt.legend(fontsize=10, loc=4)
    plt.xlim(-0.5 * p, 0.5 * p)
    plt.xlabel("time since transit [days]")
    plt.ylabel("relative flux")
    plt.title("planet {0}".format(letter))
    plt.xlim(-0.3, 0.3)
../../_images/11c89d79bd85a92dc5c4d889de180a43e6387df435c59f66d650c53ce37a5c7d.png ../../_images/e5b948f80e5683aad137305496ba531df60f67f1c6f84b89930e2571e73c5b1c.png

That looks better!

Sampling#

Now let’s run some MCMC as usual:

np.random.seed(230948)
with model:
    trace = pmx.sample(
        tune=1000,
        draws=1000,
        start=map_soln,
        cores=2,
        chains=2,
        return_inferencedata=True,
    )
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [tts_1, tts_0, b, logr, u, mean]
100.00% [4000/4000 01:29<00:00 Sampling 2 chains, 0 divergences]
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 91 seconds.

Then check the convergence diagnostics:

import arviz as az

az.summary(trace, var_names=["mean", "u", "logr", "b", "tts_0", "tts_1"])
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mean -0.000 0.000 -0.000 0.000 0.000 0.000 2908.0 1633.0 1.00
u[0] 0.341 0.175 0.010 0.612 0.004 0.003 1423.0 1041.0 1.00
u[1] 0.186 0.303 -0.302 0.724 0.009 0.007 1032.0 764.0 1.00
logr[0] -3.224 0.019 -3.260 -3.189 0.000 0.000 1745.0 1366.0 1.00
logr[1] -2.814 0.013 -2.840 -2.791 0.000 0.000 1251.0 768.0 1.00
b[0] 0.391 0.044 0.311 0.470 0.001 0.001 1615.0 1051.0 1.00
b[1] 0.353 0.032 0.288 0.406 0.001 0.001 1265.0 604.0 1.00
tts_0[0] 6.963 0.005 6.955 6.971 0.000 0.000 1884.0 1673.0 1.00
tts_0[1] 15.105 0.007 15.090 15.115 0.000 0.000 1710.0 1776.0 1.00
tts_0[2] 23.380 0.002 23.376 23.385 0.000 0.000 2217.0 1562.0 1.00
tts_0[3] 31.667 0.004 31.659 31.673 0.000 0.000 2826.0 1642.0 1.00
tts_0[4] 39.813 0.003 39.808 39.819 0.000 0.000 3633.0 1532.0 1.00
tts_0[5] 48.105 0.003 48.099 48.111 0.000 0.000 2899.0 1730.0 1.00
tts_0[6] 56.371 0.003 56.366 56.378 0.000 0.000 3055.0 1683.0 1.01
tts_0[7] 64.525 0.002 64.522 64.529 0.000 0.000 3346.0 1397.0 1.00
tts_0[8] 72.835 0.003 72.830 72.840 0.000 0.000 3096.0 1519.0 1.00
tts_1[0] 1.971 0.002 1.967 1.973 0.000 0.000 3102.0 1270.0 1.00
tts_1[1] 13.395 0.001 13.393 13.398 0.000 0.000 2606.0 1451.0 1.00
tts_1[2] 24.744 0.001 24.741 24.747 0.000 0.000 3289.0 1748.0 1.00
tts_1[3] 36.144 0.001 36.141 36.146 0.000 0.000 2958.0 1359.0 1.00
tts_1[4] 47.512 0.001 47.509 47.515 0.000 0.000 3431.0 1582.0 1.00
tts_1[5] 58.889 0.001 58.886 58.891 0.000 0.000 2927.0 1532.0 1.00
tts_1[6] 70.284 0.001 70.282 70.287 0.000 0.000 2623.0 1393.0 1.00

And plot the corner plot of the physical parameters:

import corner

names = ["period", "r", "b"]

with model:
    truths = dict(zip(names, pmx.eval_in_model([orbit.period, r, b])))

_ = corner.corner(
    trace,
    truths=truths,
    var_names=names,
)
../../_images/359461607ec785ac267e602dcfb7351b4840f9e558f2acea0e7be6bf0bbe7a71.png

We could also plot corner plots of the transit times, but they’re not terribly enlightening in this case so let’s skip it.

Finally, let’s plot the posterior estimates of the the transit times in an O-C diagram:

flat_samps = trace.posterior.stack(sample=("chain", "draw"))

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5), sharex=True)

q = np.percentile(flat_samps["ttvs_0"], [16, 50, 84], axis=-1)
ax1.fill_between(
    np.mean(flat_samps["tts_0"], axis=-1),
    q[0],
    q[2],
    color="C0",
    alpha=0.4,
    edgecolor="none",
)
ref = np.polyval(
    np.polyfit(true_transit_times[0], true_ttvs[0], 1), true_transit_times[0]
)
ax1.plot(true_transit_times[0], true_ttvs[0] - ref, ".k")
ax1.axhline(0, color="k", lw=0.5)
ax1.set_ylim(np.max(np.abs(ax1.get_ylim())) * np.array([-1, 1]))

ax1.set_ylabel("$O-C$ [days]")

q = np.percentile(flat_samps["ttvs_1"], [16, 50, 84], axis=-1)
ax2.fill_between(
    np.mean(flat_samps["tts_1"], axis=-1),
    q[0],
    q[2],
    color="C1",
    alpha=0.4,
    edgecolor="none",
)
ref = np.polyval(
    np.polyfit(true_transit_times[1], true_ttvs[1], 1), true_transit_times[1]
)
ax2.plot(true_transit_times[1], true_ttvs[1] - ref, ".k", label="truth")
ax2.axhline(0, color="k", lw=0.5)
ax2.set_ylim(np.max(np.abs(ax2.get_ylim())) * np.array([-1, 1]))

ax2.legend(fontsize=10)
ax2.set_ylabel("$O-C$ [days]")
ax2.set_xlabel("transit time [days]")
_ = ax1.set_title("posterior inference")
../../_images/a8caa1ef2c4055df91a61e4009a18bf82e03babac68373f33613ed6dfa6bf1a7.png

Citations#

As described in the citation tutorial, we can use citations.get_citations_for_model to construct an acknowledgement and BibTeX listing that includes the relevant citations for this model.

with model:
    txt, bib = xo.citations.get_citations_for_model()
print(txt)
This research made use of \textsf{exoplanet} \citep{exoplanet:joss,
exoplanet:zenodo} and its dependencies \citep{exoplanet:agol20,
exoplanet:arviz, exoplanet:astropy13, exoplanet:astropy18, exoplanet:kipping13,
exoplanet:luger18, exoplanet:pymc3, exoplanet:theano}.
print(bib.split("\n\n")[0] + "\n\n...")
@article{exoplanet:joss,
       author = {{Foreman-Mackey}, Daniel and {Luger}, Rodrigo and {Agol}, Eric
                and {Barclay}, Thomas and {Bouma}, Luke G. and {Brandt},
                Timothy D. and {Czekala}, Ian and {David}, Trevor J. and
                {Dong}, Jiayin and {Gilbert}, Emily A. and {Gordon}, Tyler A.
                and {Hedges}, Christina and {Hey}, Daniel R. and {Morris},
                Brett M. and {Price-Whelan}, Adrian M. and {Savel}, Arjun B.},
        title = "{exoplanet: Gradient-based probabilistic inference for
                  exoplanet data \& other astronomical time series}",
      journal = {arXiv e-prints},
         year = 2021,
        month = may,
          eid = {arXiv:2105.01994},
        pages = {arXiv:2105.01994},
archivePrefix = {arXiv},
       eprint = {2105.01994},
 primaryClass = {astro-ph.IM},
       adsurl = {https://ui.adsabs.harvard.edu/abs/2021arXiv210501994F},
      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}

...