Note

This tutorial was generated from an IPython notebook that can be downloaded here.

# Transit fitting¶

exoplanet includes methods for computing the light curves transiting planets. In its simplest form this can be used to evaluate a light curve like you would do with batman, for example:

import numpy as np
import matplotlib.pyplot as plt

import exoplanet as xo

# The light curve calculation requires an orbit
orbit = xo.orbits.KeplerianOrbit(period=3.456)

# Compute a limb-darkened light curve using starry
t = np.linspace(-0.1, 0.1, 1000)
u = [0.3, 0.2]
light_curve = (
xo.LimbDarkLightCurve(u)
.get_light_curve(orbit=orbit, r=0.1, t=t, texp=0.02)
.eval()
)
# Note: the eval is needed because this is using Theano in
# the background

plt.plot(t, light_curve, color="C0", lw=2)
plt.ylabel("relative flux")
plt.xlabel("time [days]")
_ = plt.xlim(t.min(), t.max()) But the real power comes from the fact that this is defined as a Theano operation so it can be combined with PyMC3 to do transit inference using Hamiltonian Monte Carlo.

## The transit model in PyMC3¶

In this section, we will construct a simple transit fit model using PyMC3 and then we will fit a two planet model to simulated data. To start, let’s randomly sample some periods and phases and then define the time sampling:

np.random.seed(123)
periods = np.random.uniform(5, 20, 2)
t0s = periods * np.random.rand(2)
t = np.arange(0, 80, 0.02)
yerr = 5e-4


Then, define the parameters. In this simple model, we’ll just fit for the limb darkening parameters of the star, and the period, phase, impact parameter, and radius ratio of the planets (note: this is already 10 parameters and running MCMC to convergence using emcee would probably take at least an hour). For the limb darkening, we’ll use a quadratic law as parameterized by Kipping (2013). This reparameterizations is implemented in exoplanet as custom PyMC3 distribution exoplanet.distributions.QuadLimbDark.

import pymc3 as pm

with pm.Model() as model:

# The baseline flux
mean = pm.Normal("mean", mu=0.0, sd=1.0)

# The time of a reference transit for each planet
t0 = pm.Normal("t0", mu=t0s, sd=1.0, shape=2)

# The log period; also tracking the period itself
logP = pm.Normal("logP", mu=np.log(periods), sd=0.1, shape=2)
period = pm.Deterministic("period", pm.math.exp(logP))

# The Kipping (2013) parameterization for quadratic limb darkening paramters

r = pm.Uniform(
"r", lower=0.01, upper=0.1, shape=2, testval=np.array([0.04, 0.06])
)
b = xo.distributions.ImpactParameter(
"b", ror=r, shape=2, testval=np.random.rand(2)
)

# Set up a Keplerian orbit for the planets
orbit = xo.orbits.KeplerianOrbit(period=period, t0=t0, b=b)

# Compute the model light curve using starry
light_curves = xo.LimbDarkLightCurve(u).get_light_curve(
orbit=orbit, r=r, t=t
)
light_curve = pm.math.sum(light_curves, axis=-1) + mean

# Here we track the value of the model light curve for plotting
# purposes
pm.Deterministic("light_curves", light_curves)

# In this line, we simulate the dataset that we will fit
y = xo.eval_in_model(light_curve)
y += yerr * np.random.randn(len(y))

# The likelihood function assuming known Gaussian uncertainty
pm.Normal("obs", mu=light_curve, sd=yerr, observed=y)

# Fit for the maximum a posteriori parameters given the simuated
# dataset
map_soln = xo.optimize(start=model.test_point)

optimizing logp for variables: [b, r, u, logP, t0, mean]
100it [00:05, 16.93it/s, logp=2.479354e+04]
message: Desired error not necessarily achieved due to precision loss.
logp: 24787.977771807487 -> 24793.53942561119


Now we can plot the simulated data and the maximum a posteriori model to make sure that our initialization looks ok.

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") ## Sampling¶

Now, let’s sample from the posterior defined by this model. As usual, there are strong covariances between some of the parameters so we’ll use exoplanet.get_dense_nuts_step().

np.random.seed(42)
with model:
trace = pm.sample(
tune=3000,
draws=3000,
start=map_soln,
cores=2,
chains=2,
step=xo.get_dense_nuts_step(target_accept=0.9),
)

Multiprocess sampling (2 chains in 2 jobs)
NUTS: [b, r, u, logP, t0, mean]
Sampling 2 chains, 0 divergences: 100%|██████████| 12000/12000 [01:30<00:00, 132.79draws/s]


After sampling, it’s important that we assess convergence. We can do that using the pymc3.summary function:

pm.summary(trace, varnames=["period", "t0", "r", "b", "u", "mean"])

mean sd hpd_3% hpd_97% mcse_mean mcse_sd ess_mean ess_sd ess_bulk ess_tail r_hat
period 15.448 0.002 15.443 15.452 0.000 0.000 1973.0 1973.0 2219.0 3072.0 1.0
period 9.292 0.000 9.292 9.293 0.000 0.000 5339.0 5339.0 5522.0 4097.0 1.0
t0 3.503 0.006 3.494 3.515 0.000 0.000 1927.0 1926.0 2080.0 3193.0 1.0
t0 5.121 0.001 5.119 5.124 0.000 0.000 6331.0 6330.0 6492.0 3902.0 1.0
r 0.040 0.002 0.037 0.042 0.000 0.000 2948.0 2948.0 3004.0 3249.0 1.0
r 0.058 0.001 0.057 0.060 0.000 0.000 2580.0 2580.0 2609.0 3801.0 1.0
b 0.671 0.046 0.581 0.739 0.001 0.001 1709.0 1709.0 1957.0 2716.0 1.0
b 0.402 0.038 0.334 0.472 0.001 0.001 2606.0 2606.0 2745.0 2906.0 1.0
u 0.375 0.209 0.002 0.717 0.003 0.002 4195.0 4155.0 4108.0 3492.0 1.0
u 0.273 0.343 -0.327 0.842 0.007 0.005 2663.0 2348.0 2584.0 2841.0 1.0
mean 0.000 0.000 -0.000 0.000 0.000 0.000 7178.0 3124.0 7182.0 3887.0 1.0

That looks pretty good! Fitting this without exoplanet would have taken a lot more patience.

Now we can also look at the corner plot of some of that parameters of interest:

import corner

samples = pm.trace_to_dataframe(trace, varnames=["period", "r"])
truth = np.concatenate(
xo.eval_in_model([period, r], model.test_point, model=model)
)
_ = corner.corner(
samples,
truths=truth,
) ## Phase plots¶

Like in the radial velocity tutorial (Radial velocity fitting), we can make plots of the model predictions for each planet.

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

# Get the posterior median orbital parameters
p = np.median(trace["period"][:, n])
t0 = np.median(trace["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 = np.median(trace["light_curves"][:, :, (n + 1) % 2], axis=0)

# 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
)

# Plot the folded model
inds = np.argsort(x_fold)
inds = inds[np.abs(x_fold)[inds] < 0.3]
pred = trace["light_curves"][:, inds, n] + trace["mean"][:, None]
pred = np.median(pred, axis=0)
plt.plot(x_fold[inds], pred, color="C1", label="model")

# Annotate the plot with the planet's period
txt = "period = {0:.4f} +/- {1:.4f} d".format(
np.mean(trace["period"][:, n]), np.std(trace["period"][:, n])
)
plt.annotate(
txt,
(0, 0),
xycoords="axes fraction",
xytext=(5, 5),
textcoords="offset points",
ha="left",
va="bottom",
fontsize=12,
)

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)  