Note

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

# Fitting LHS 3844b¶

In this tutorial, we will reproduce the fits to the transiting planet in the LHS 3844 system discovered by Vanderspek et al. (2018). This is the same as the Pi Mensae on the exoplanet documentation with a few small changes.

import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt

tpf_url = "https://archive.stsci.edu/missions/tess/tid/s0001/0000/0004/1015/3553/tess2018206045859-s0001-0000000410153553-0120-s_tp.fits"
with fits.open(tpf_url) as hdus:
tpf = hdus[1].data

texp = tpf_hdr["FRAMETIM"] * tpf_hdr["NUM_FRM"]
texp /= 60.0 * 60.0 * 24.0
time = tpf["TIME"]
flux = tpf["FLUX"][:, :, 1:]
m = np.any(np.isfinite(flux), axis=(1, 2)) & (tpf["QUALITY"] == 0)
time = np.ascontiguousarray(time[m] - np.min(time[m]), dtype=np.float64)
flux = np.ascontiguousarray(flux[m], dtype=np.float64)

mean_img = np.median(flux, axis=0)
plt.imshow(mean_img.T, cmap="gray_r")
plt.title("TESS image of LHS 3844")
plt.xticks([])
plt.yticks([]);


## Aperture selection¶

Next, we’ll select an aperture using a hacky method that tries to minimizes the windowed scatter in the lightcurve (something like the CDPP).

from scipy.signal import savgol_filter

# Sort the pixels by median brightness
order = np.argsort(mean_img.flatten())[::-1]

# A function to estimate the windowed scatter in a lightcurve
smooth = savgol_filter(f, 1001, polyorder=5)
return 1e6 * np.sqrt(np.median((f / smooth - 1)**2))

# Loop over pixels ordered by brightness and add them one-by-one
# to the aperture
for i in range(10, 100):
msk = np.zeros_like(mean_img, dtype=bool)
msk[np.unravel_index(order[:i], mean_img.shape)] = True
scatters.append(scatter)

# Choose the aperture that minimizes the scatter

# Plot the selected aperture
plt.imshow(mean_img.T, cmap="gray_r")
plt.title("selected aperture")
plt.xticks([])
plt.yticks([]);


This aperture produces the following light curve:

plt.figure(figsize=(10, 5))
sap_flux = (sap_flux / np.median(sap_flux) - 1) * 1e3
plt.plot(time, sap_flux, "k")
plt.xlabel("time [days]")
plt.ylabel("relative flux [ppt]")
plt.title("raw light curve")
plt.xlim(time.min(), time.max());


## The transit model in PyMC3¶

The transit model, initialization, and sampling are all the same as in the TESS tutorial on the exoplanet docs.

import exoplanet as xo
import pymc3 as pm
import theano.tensor as tt

with pm.Model() as model:

# Parameters for the stellar properties
mean = pm.Normal("mean", mu=0.0, sd=10.0)

# Stellar parameters from Vanderspek et al (2018)
M_star_vanderspek = 0.151, 0.014
R_star_vanderspek = 0.189, 0.006
m_star = pm.Normal("m_star", mu=M_star_vanderspek[0], sd=M_star_vanderspek[1])
r_star = pm.Normal("r_star", mu=R_star_vanderspek[0], sd=R_star_vanderspek[1])

# Prior to require physical parameters
pm.Potential("m_star_prior", tt.switch(m_star > 0, 0, -np.inf))
pm.Potential("r_star_prior", tt.switch(r_star > 0, 0, -np.inf))

# Orbital parameters for the planets
logP = pm.Normal("logP", mu=np.log(bls_period), sd=1)
t0 = pm.Normal("t0", mu=bls_t0, sd=1)
testval_r=np.sqrt(1e-3)*np.sqrt(bls_depth),
testval_b=0.5)
ecc = pm.Beta("ecc", alpha=0.867, beta=3.03, testval=0.1)
omega = xo.distributions.Angle("omega")

# Log-uniform prior on ror
pm.Potential("ror_prior", -tt.log(ror))

# Transit jitter & GP parameters
logw0 = pm.Normal("logw0", mu=np.log(2*np.pi/10), sd=10)

# Tracking planet parameters
period = pm.Deterministic("period", tt.exp(logP))
r_pl = pm.Deterministic("r_pl", r_star * ror)

# Orbit model
orbit = xo.orbits.KeplerianOrbit(
r_star=r_star, m_star=m_star,
period=period, t0=t0, b=b,
ecc=ecc, omega=omega)

pm.Deterministic("a", orbit.a)
pm.Deterministic("rho_star", orbit.rho_star)

# Compute the model light curve using starry
light_curves = xo.StarryLightCurve(u_star, r_star=r_star).get_light_curve(
light_curve = pm.math.sum(light_curves, axis=-1) + mean
pm.Deterministic("light_curves", light_curves)

# GP model for the light curve
kernel = xo.gp.terms.SHOTerm(log_S0=logS0, log_w0=logw0, Q=1/np.sqrt(2))
pm.Deterministic("gp_pred", gp.predict())

# Fit for the maximum a posteriori parameters, I've found that I can get
# a better solution by trying different combinations of parameters in turn
if start is None:
start = model.test_point
map_soln = pm.find_MAP(start=start, vars=[logs2, logS0, logw0])
map_soln = pm.find_MAP(start=map_soln, vars=[model.rb])
map_soln = pm.find_MAP(start=map_soln)

return model, map_soln

model0, map_soln0 = build_model()

logp = -46,523, ||grad|| = 21.343: 100%|██████████| 25/25 [00:00<00:00, 74.26it/s]
logp = -46,492, ||grad|| = 0.032486: 100%|██████████| 37/37 [00:00<00:00, 81.48it/s]
logp = -46,459, ||grad|| = 1.1253e+05: 100%|██████████| 45/45 [00:00<00:00, 62.91it/s]


Here’s how we plot the initial light curve model:

def plot_light_curve(soln, mask=None):

fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)

ax = axes[0]
gp_mod = soln["gp_pred"] + soln["mean"]
ax.legend(fontsize=10)
ax.set_ylabel("relative flux [ppt]")

ax = axes[1]
for i, l in enumerate("b"):
mod = soln["light_curves"][:, i]
ax.legend(fontsize=10, loc=3)
ax.set_ylabel("de-trended flux [ppt]")

ax = axes[2]
mod = gp_mod + np.sum(soln["light_curves"], axis=-1)
ax.axhline(0, color="#aaaaaa", lw=1)
ax.set_ylabel("residuals [ppt]")
ax.set_xlabel("time [days]")

return fig

plot_light_curve(map_soln0);


We can do some sigma clipping to remove significant outliers.

mod = map_soln0["gp_pred"] + map_soln0["mean"] + np.sum(map_soln0["light_curves"], axis=-1)
resid = y - mod
rms = np.sqrt(np.median(resid**2))
mask = np.abs(resid) < 5 * rms

plt.figure(figsize=(10, 5))
plt.plot(x, resid, "k", label="data")
plt.axhline(0, color="#aaaaaa", lw=1)
plt.ylabel("residuals [ppt]")
plt.xlabel("time [days]")
plt.legend(fontsize=12, loc=3)
plt.xlim(x.min(), x.max());


And then we re-build the model using the data without outliers.

model, map_soln = build_model(mask, map_soln0)

logp = -46,137, ||grad|| = 232.01: 100%|██████████| 11/11 [00:00<00:00, 75.27it/s]
logp = -46,134, ||grad|| = 17.169: 100%|██████████| 6/6 [00:00<00:00, 89.54it/s]
logp = -46,134, ||grad|| = 1,490.3: 100%|██████████| 8/8 [00:00<00:00, 74.52it/s]


Now that we have the model, we can sample it using a exoplanet.PyMC3Sampler:

np.random.seed(42)
sampler = xo.PyMC3Sampler(window=100, start=200, finish=200)
with model:
burnin = sampler.tune(tune=3000, start=map_soln, step_kwargs=dict(target_accept=0.9))

Sampling 4 chains: 100%|██████████| 808/808 [08:36<00:00,  2.26s/draws]
Sampling 4 chains: 100%|██████████| 408/408 [04:12<00:00,  2.01s/draws]
Sampling 4 chains: 100%|██████████| 808/808 [09:02<00:00,  5.70s/draws]
The chain contains only diverging samples. The model is probably misspecified.
The chain contains only diverging samples. The model is probably misspecified.
Sampling 4 chains: 100%|██████████| 1608/1608 [23:19<00:00,  1.54s/draws]
Sampling 4 chains: 100%|██████████| 8408/8408 [1:24:57<00:00,  1.12s/draws]

with model:
trace = sampler.sample(draws=2000)

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [logw0, logS0, logs2, omega, ecc, rb, t0, logP, r_star, m_star, u_star, mean]
Sampling 4 chains: 100%|██████████| 8800/8800 [42:39<00:00,  1.12s/draws]
There were 137 divergences after tuning. Increase target_accept or reparameterize.
There were 90 divergences after tuning. Increase target_accept or reparameterize.
There were 102 divergences after tuning. Increase target_accept or reparameterize.
There were 65 divergences after tuning. Increase target_accept or reparameterize.

pm.summary(trace, varnames=["logw0", "logS0", "logs2", "omega", "ecc", "rb", "t0", "logP", "r_star", "m_star", "u_star", "mean"])

mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
logw0 0.789575 0.186396 2.727383e-03 0.417512 1.147283 4950.234236 1.000526
logS0 -0.464405 0.442049 6.192979e-03 -1.277101 0.407132 4461.844862 1.000312
logs2 2.258889 0.010632 1.458397e-04 2.238001 2.279675 5598.519190 0.999966
omega -0.207643 1.897765 3.291704e-02 -3.141438 2.939382 3208.214587 1.000389
ecc 0.106674 0.109967 2.098153e-03 0.000015 0.336140 2605.775464 1.000295
rb__0_0 0.062219 0.001485 2.264017e-05 0.059169 0.064892 3999.789649 1.000333
rb__1_0 0.255344 0.146670 2.821448e-03 0.000422 0.504333 2464.240599 1.000604
t0 0.426613 0.000292 4.046495e-06 0.426040 0.427182 4708.325766 0.999975
logP -0.770190 0.000017 2.241095e-07 -0.770222 -0.770157 4780.473804 1.000041
r_star 0.189080 0.005672 7.409012e-05 0.178541 0.200236 4668.210096 1.000097
m_star 0.150736 0.013453 1.960888e-04 0.124422 0.176856 4459.340468 1.000334
u_star__0 0.284330 0.189327 2.876975e-03 0.000067 0.629356 4670.050619 1.000213
u_star__1 0.144143 0.270250 3.936996e-03 -0.306457 0.699236 4287.876925 0.999807
mean -0.095182 0.217488 3.095587e-03 -0.545166 0.320977 4838.564147 0.999803

## Results¶

After sampling, we can make the usual plots. First, let’s look at the folded light curve plot:

# Compute the GP prediction
gp_mod = np.median(trace["gp_pred"] + trace["mean"][:, None], axis=0)

# Get the posterior median orbital parameters
p = np.median(trace["period"])
t0 = np.median(trace["t0"])

# Plot the folded data
x_fold = (x[mask] - t0 + 0.5*p) % p - 0.5*p
plt.plot(x_fold * 24., y[mask] - gp_mod, ".k", label="data",
alpha=0.5, zorder=-1000, mec="none")

# Overplot the phase binned light curve
bins = np.linspace(-0.11, 0.11, 50)
denom, _ = np.histogram(x_fold, bins)
num, _ = np.histogram(x_fold, bins, weights=y[mask])
denom[num == 0] = 1.0
plt.plot(0.5*(bins[1:] + bins[:-1]) *  24., num / denom, "o", color="C2",
label="binned")

# Plot the folded model
inds = np.argsort(x_fold)
inds = inds[np.abs(x_fold)[inds] < 0.3]
pred = trace["light_curves"][:, inds, 0]
pred = np.percentile(pred, [16, 50, 84], axis=0)
plt.plot(x_fold[inds] * 24., pred[1], color="C1", label="model")
art = plt.fill_between(x_fold[inds] * 24., pred[0], pred[2], color="C1", alpha=0.5,
zorder=1000)
art.set_edgecolor("none")

# Annotate the plot with the planet's period
txt = "period = {0:.5f} +/- {1:.5f} hr".format(
np.mean(trace["period"] * 24.), np.std(trace["period"] * 24.))
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.xlabel("time since transit [hours]")
plt.ylabel("de-trended flux")
plt.xlim(-2, 2);


And a corner plot of some of the key parameters:

import corner
import astropy.units as u
varnames = ["period", "b", "ecc", "r_pl"]
samples = pm.trace_to_dataframe(trace, varnames=varnames)