Note
This tutorial was generated from an IPython notebook that can be downloaded here.
Note
You will need exoplanet version 0.3.1 or later to run this tutorial.
In this case study, we’ll go through the steps required to fit the light curve and radial velocity measurements for the detached eclipsing binary system HD 23642. This is a bright system that has been fit by many authors (1, 2, 3, 4, and 5 to name a few) so this is a good benchmark for our demonstration.
The light curve that we’ll use is from K2 and we’ll use the same radial velocity measurements as David+ (2016) compiled from here and here. We’ll use a somewhat simplified model for the eclipses that treats the stars as spherical and ignores the phase curve (we’ll model it using a Gaussian process instead of a more physically motivated model). But, as you’ll see, these simplifying assumptions are sufficient for this case of a detached and well behaved system. Unlike some previous studies, we will fit an eccentric orbit instead of fixing the eccentricity to zero. This probably isn’t really necessary here, but it’s useful to demonstrate how you would fit a more eccentric system. Finally, we model the phase curve and other triends in both the light curve and radial velocities using Gaussian processes. This will account for unmodeled stellar variability and residual systematics, drifts, and other effects left over from the data reduction procedure.
First, let’s define some values from the literature that will be useful below. Here we’re taking the period and eclipse time from David+ (2016) as initial guesses for these parameters in our fit. We’ll also include the same prior on the flux ratio of the two stars that was computed for the Kepler bandpass by David+ (2016).
lit_period = 2.46113408
lit_t0 = 119.522070 + 2457000 - 2454833
# Prior on the flux ratio for Kepler
lit_flux_ratio = (0.354, 0.035)
Then we’ll download the Kepler data. In this case, the pipeline aperture photometry isn’t very good (because this star is so bright!) so we’ll just download the target pixel file and co-add all the pixels.
import numpy as np
import matplotlib.pyplot as plt
import lightkurve as lk
tpf = lk.search_targetpixelfile("EPIC 211082420").download()
lc = tpf.to_lightcurve(aperture_mask="all")
lc = lc.remove_nans().normalize()
hdr = tpf.hdu[1].header
texp = hdr["FRAMETIM"] * hdr["NUM_FRM"]
texp /= 60.0 * 60.0 * 24.0
x = np.ascontiguousarray(lc.time, dtype=np.float64)
y = np.ascontiguousarray(lc.flux, dtype=np.float64)
mu = np.median(y)
y = (y / mu - 1) * 1e3
plt.plot((x - lit_t0 + 0.5 * lit_period) % lit_period - 0.5 * lit_period, y, ".k")
plt.xlim(-0.5 * lit_period, 0.5 * lit_period)
plt.xlabel("time since primary eclipse [days]")
_ = plt.ylabel("relative flux [ppt]")
Then we’ll enter the radial velocity data. I couldn’t find these data online anywhere so I manually transcribed the data from the referenced papers (typos are my own!).
ref1 = 2453000
ref2 = 2400000
rvs = np.array(
[
# https://arxiv.org/abs/astro-ph/0403444
(39.41273 + ref1, -85.0, 134.5),
(39.45356 + ref1, -88.0, 139.0),
(39.50548 + ref1, -91.0, 143.0),
(43.25049 + ref1, 105.5, -136.0),
(46.25318 + ref1, 29.5, -24.5),
# https://ui.adsabs.harvard.edu/abs/2007A%26A...463..579G/abstract
(52629.6190 + ref2, 88.8, -127.0),
(52630.6098 + ref2, -48.0, 68.0),
(52631.6089 + ref2, -9.5, 13.1),
(52632.6024 + ref2, 63.6, -90.9),
(52633.6162 + ref2, -94.5, 135.0),
(52636.6055 + ref2, 10.3, -13.9),
(52983.6570 + ref2, 18.1, -25.1),
(52987.6453 + ref2, -80.6, 114.5),
(52993.6322 + ref2, 49.0, -70.7),
(53224.9338 + ref2, 39.0, -55.7),
(53229.9384 + ref2, 57.2, -82.0),
]
)
rvs[:, 0] -= 2454833
rvs = rvs[np.argsort(rvs[:, 0])]
x_rv = np.ascontiguousarray(rvs[:, 0], dtype=np.float64)
y1_rv = np.ascontiguousarray(rvs[:, 1], dtype=np.float64)
y2_rv = np.ascontiguousarray(rvs[:, 2], dtype=np.float64)
fold = (rvs[:, 0] - lit_t0 + 0.5 * lit_period) % lit_period - 0.5 * lit_period
plt.plot(fold, rvs[:, 1], ".", label="primary")
plt.plot(fold, rvs[:, 2], ".", label="secondary")
plt.legend(fontsize=10)
plt.xlim(-0.5 * lit_period, 0.5 * lit_period)
plt.ylabel("radial velocity [km / s]")
_ = plt.xlabel("time since primary eclipse [days]")
Then we define the probabilistic model using PyMC3 and exoplanet. This
is similar to the other tutorials and case studies, but here we’re using
a exoplanet.SecondaryEclipseLightCurve
to generate the model
light curve and we’re modeling the radial velocity trends using a
Gaussian process instead of a polynomial. Otherwise, things should look
pretty familiar!
After defining the model, we iteratively clip outliers in the light curve using sigma clipping and then estimate the maximum a posteriori parameters.
import pymc3 as pm
import theano.tensor as tt
import exoplanet as xo
def build_model(mask):
with pm.Model() as model:
# Systemic parameters
mean_lc = pm.Normal("mean_lc", mu=0.0, sd=5.0)
mean_rv = pm.Normal("mean_rv", mu=0.0, sd=50.0)
u1 = xo.QuadLimbDark("u1")
u2 = xo.QuadLimbDark("u2")
# Parameters describing the primary
M1 = pm.Lognormal("M1", mu=0.0, sigma=10.0)
R1 = pm.Lognormal("R1", mu=0.0, sigma=10.0)
# Secondary ratios
k = pm.Lognormal("k", mu=0.0, sigma=10.0) # radius ratio
q = pm.Lognormal("q", mu=0.0, sigma=10.0) # mass ratio
s = pm.Lognormal("s", mu=np.log(0.5), sigma=10.0) # surface brightness ratio
# Prior on flux ratio
pm.Normal(
"flux_prior", mu=lit_flux_ratio[0], sigma=lit_flux_ratio[1], observed=k * s
)
# Parameters describing the orbit
b = xo.ImpactParameter("b", ror=k, testval=1.5)
period = pm.Lognormal("period", mu=np.log(lit_period), sigma=1.0)
t0 = pm.Normal("t0", mu=lit_t0, sigma=1.0)
# Parameters describing the eccentricity: ecs = [e * cos(w), e * sin(w)]
ecs = xo.UnitDisk("ecs", testval=np.array([1e-5, 0.0]))
ecc = pm.Deterministic("ecc", tt.sqrt(tt.sum(ecs ** 2)))
omega = pm.Deterministic("omega", tt.arctan2(ecs[1], ecs[0]))
# Build the orbit
R2 = pm.Deterministic("R2", k * R1)
M2 = pm.Deterministic("M2", q * M1)
orbit = xo.orbits.KeplerianOrbit(
period=period,
t0=t0,
ecc=ecc,
omega=omega,
b=b,
r_star=R1,
m_star=M1,
m_planet=M2,
)
# Track some other orbital elements
pm.Deterministic("incl", orbit.incl)
pm.Deterministic("a", orbit.a)
# Noise model for the light curve
sigma_lc = pm.InverseGamma(
"sigma_lc", testval=1.0, **xo.estimate_inverse_gamma_parameters(0.1, 2.0)
)
S_tot_lc = pm.InverseGamma(
"S_tot_lc", testval=2.5, **xo.estimate_inverse_gamma_parameters(1.0, 5.0)
)
ell_lc = pm.InverseGamma(
"ell_lc", testval=2.0, **xo.estimate_inverse_gamma_parameters(1.0, 5.0)
)
kernel_lc = xo.gp.terms.SHOTerm(
S_tot=S_tot_lc, w0=2 * np.pi / ell_lc, Q=1.0 / 3
)
# Noise model for the radial velocities
sigma_rv1 = pm.InverseGamma(
"sigma_rv1", testval=1.0, **xo.estimate_inverse_gamma_parameters(0.5, 5.0)
)
sigma_rv2 = pm.InverseGamma(
"sigma_rv2", testval=1.0, **xo.estimate_inverse_gamma_parameters(0.5, 5.0)
)
S_tot_rv = pm.InverseGamma(
"S_tot_rv", testval=2.5, **xo.estimate_inverse_gamma_parameters(1.0, 5.0)
)
ell_rv = pm.InverseGamma(
"ell_rv", testval=2.0, **xo.estimate_inverse_gamma_parameters(1.0, 5.0)
)
kernel_rv = xo.gp.terms.SHOTerm(
S_tot=S_tot_rv, w0=2 * np.pi / ell_rv, Q=1.0 / 3
)
# Set up the light curve model
lc = xo.SecondaryEclipseLightCurve(u1, u2, s)
def model_lc(t):
return (
mean_lc
+ 1e3 * lc.get_light_curve(orbit=orbit, r=R2, t=t, texp=texp)[:, 0]
)
# Condition the light curve model on the data
gp_lc = xo.gp.GP(
kernel_lc, x[mask], tt.zeros(mask.sum()) ** 2 + sigma_lc ** 2, mean=model_lc
)
gp_lc.marginal("obs_lc", observed=y[mask])
# Set up the radial velocity model
def model_rv1(t):
return mean_rv + 1e-3 * orbit.get_radial_velocity(t)
def model_rv2(t):
return mean_rv - 1e-3 * orbit.get_radial_velocity(t) / q
# Condition the radial velocity model on the data
gp_rv1 = xo.gp.GP(
kernel_rv, x_rv, tt.zeros(len(x_rv)) ** 2 + sigma_rv1 ** 2, mean=model_rv1
)
gp_rv1.marginal("obs_rv1", observed=y1_rv)
gp_rv2 = xo.gp.GP(
kernel_rv, x_rv, tt.zeros(len(x_rv)) ** 2 + sigma_rv2 ** 2, mean=model_rv2
)
gp_rv2.marginal("obs_rv2", observed=y2_rv)
# Optimize the logp
map_soln = model.test_point
# First the RV parameters
map_soln = xo.optimize(map_soln, [mean_rv, q])
map_soln = xo.optimize(
map_soln, [mean_rv, sigma_rv1, sigma_rv2, S_tot_rv, ell_rv]
)
# Then the LC parameters
map_soln = xo.optimize(map_soln, [mean_lc, R1, k, s, b])
map_soln = xo.optimize(map_soln, [mean_lc, R1, k, s, b, u1, u2])
map_soln = xo.optimize(map_soln, [mean_lc, sigma_lc, S_tot_lc, ell_lc])
map_soln = xo.optimize(map_soln, [t0, period])
# Then all the parameters together
map_soln = xo.optimize(map_soln)
model.gp_lc = gp_lc
model.model_lc = model_lc
model.gp_rv1 = gp_rv1
model.model_rv1 = model_rv1
model.gp_rv2 = gp_rv2
model.model_rv2 = model_rv2
model.x = x[mask]
model.y = y[mask]
return model, map_soln
def sigma_clip():
mask = np.ones(len(x), dtype=bool)
num = len(mask)
for i in range(10):
model, map_soln = build_model(mask)
with model:
mdl = xo.eval_in_model(
model.model_lc(x[mask]) + model.gp_lc.predict(), map_soln
)
resid = y[mask] - mdl
sigma = np.sqrt(np.median((resid - np.median(resid)) ** 2))
mask[mask] = np.abs(resid - np.median(resid)) < 7 * sigma
print("Sigma clipped {0} light curve points".format(num - mask.sum()))
if num == mask.sum():
break
num = mask.sum()
return model, map_soln
model, map_soln = sigma_clip()
optimizing logp for variables: [q, mean_rv]
13it [00:03, 3.85it/s, logp=-1.307282e+04]
message: Optimization terminated successfully.
logp: -23380.5842277239 -> -13072.816791759677
optimizing logp for variables: [ell_rv, S_tot_rv, sigma_rv2, sigma_rv1, mean_rv]
40it [00:00, 87.85it/s, logp=-8.482576e+03]
message: Optimization terminated successfully.
logp: -13072.816791759677 -> -8482.576126912147
optimizing logp for variables: [b, k, s, R1, mean_lc]
39it [00:00, 90.36it/s, logp=-4.696693e+03]
message: Optimization terminated successfully.
logp: -8482.576126912147 -> -4696.693195993793
optimizing logp for variables: [u2, u1, b, k, s, R1, mean_lc]
40it [00:00, 46.45it/s, logp=-4.693636e+03]
message: Optimization terminated successfully.
logp: -4696.693195993793 -> -4693.636423265372
optimizing logp for variables: [ell_lc, S_tot_lc, sigma_lc, mean_lc]
84it [00:00, 112.79it/s, logp=-3.640328e+03]
message: Desired error not necessarily achieved due to precision loss.
logp: -4693.6364232653705 -> -3640.327844208639
optimizing logp for variables: [period, t0]
101it [00:00, 162.06it/s, logp=-3.636692e+03]
message: Desired error not necessarily achieved due to precision loss.
logp: -3640.3278442086394 -> -3636.6924367320316
optimizing logp for variables: [ell_rv, S_tot_rv, sigma_rv2, sigma_rv1, ell_lc, S_tot_lc, sigma_lc, ecs, t0, period, b, s, q, k, R1, M1, u2, u1, mean_rv, mean_lc]
171it [00:01, 147.70it/s, logp=-3.539690e+03]
message: Desired error not necessarily achieved due to precision loss.
logp: -3636.6924367320316 -> -3539.689989722382
Sigma clipped 34 light curve points
optimizing logp for variables: [q, mean_rv]
13it [00:00, 59.37it/s, logp=-1.271094e+04]
message: Optimization terminated successfully.
logp: -22965.237266602366 -> -12710.941974601861
optimizing logp for variables: [ell_rv, S_tot_rv, sigma_rv2, sigma_rv1, mean_rv]
40it [00:00, 127.35it/s, logp=-8.141907e+03]
message: Optimization terminated successfully.
logp: -12710.941974601868 -> -8141.907237100719
optimizing logp for variables: [b, k, s, R1, mean_lc]
42it [00:00, 125.83it/s, logp=-4.369324e+03]
message: Optimization terminated successfully.
logp: -8141.907237100712 -> -4369.324119375993
optimizing logp for variables: [u2, u1, b, k, s, R1, mean_lc]
42it [00:00, 119.98it/s, logp=-4.366186e+03]
message: Optimization terminated successfully.
logp: -4369.324119375993 -> -4366.186197582646
optimizing logp for variables: [ell_lc, S_tot_lc, sigma_lc, mean_lc]
23it [00:00, 39.73it/s, logp=-1.789877e+03]
message: Optimization terminated successfully.
logp: -4366.186197582645 -> -1789.876738524621
optimizing logp for variables: [period, t0]
134it [00:00, 190.84it/s, logp=-1.784739e+03]
message: Desired error not necessarily achieved due to precision loss.
logp: -1789.8767385246208 -> -1784.7391024505662
optimizing logp for variables: [ell_rv, S_tot_rv, sigma_rv2, sigma_rv1, ell_lc, S_tot_lc, sigma_lc, ecs, t0, period, b, s, q, k, R1, M1, u2, u1, mean_rv, mean_lc]
211it [00:01, 146.54it/s, logp=-1.624762e+03]
message: Desired error not necessarily achieved due to precision loss.
logp: -1784.7391024505662 -> -1624.7623102375767
Sigma clipped 1 light curve points
optimizing logp for variables: [q, mean_rv]
13it [00:00, 54.23it/s, logp=-1.270829e+04]
message: Optimization terminated successfully.
logp: -22960.077690259175 -> -12708.286451975617
optimizing logp for variables: [ell_rv, S_tot_rv, sigma_rv2, sigma_rv1, mean_rv]
40it [00:00, 118.44it/s, logp=-8.140053e+03]
message: Optimization terminated successfully.
logp: -12708.286451975624 -> -8140.052575477539
optimizing logp for variables: [b, k, s, R1, mean_lc]
43it [00:00, 123.85it/s, logp=-4.367856e+03]
message: Desired error not necessarily achieved due to precision loss.
logp: -8140.052575477532 -> -4367.855852535798
optimizing logp for variables: [u2, u1, b, k, s, R1, mean_lc]
127it [00:00, 177.62it/s, logp=-4.364686e+03]
message: Desired error not necessarily achieved due to precision loss.
logp: -4367.855852535798 -> -4364.686321204608
optimizing logp for variables: [ell_lc, S_tot_lc, sigma_lc, mean_lc]
78it [00:00, 153.05it/s, logp=-1.780554e+03]
message: Desired error not necessarily achieved due to precision loss.
logp: -4364.686321204609 -> -1780.5541237576306
optimizing logp for variables: [period, t0]
64it [00:00, 138.25it/s, logp=-1.775471e+03]
message: Desired error not necessarily achieved due to precision loss.
logp: -1780.5541237576322 -> -1775.4707961536694
optimizing logp for variables: [ell_rv, S_tot_rv, sigma_rv2, sigma_rv1, ell_lc, S_tot_lc, sigma_lc, ecs, t0, period, b, s, q, k, R1, M1, u2, u1, mean_rv, mean_lc]
423it [00:02, 177.80it/s, logp=-1.614994e+03]
message: Desired error not necessarily achieved due to precision loss.
logp: -1775.4707961536694 -> -1614.9941877137726
Sigma clipped 0 light curve points
At these best fit parameters, let’s make some plots of the model predictions compared to the observations to make sure that things look reasonable. First the phase-folded radial velocities:
period = map_soln["period"]
t0 = map_soln["t0"]
mean = map_soln["mean_rv"]
x_fold = (x_rv - t0 + 0.5 * period) % period - 0.5 * period
plt.plot(fold, y1_rv - mean, ".", label="primary")
plt.plot(fold, y2_rv - mean, ".", label="secondary")
x_phase = np.linspace(-0.5 * period, 0.5 * period, 500)
with model:
y1_mod, y2_mod = xo.eval_in_model(
[model.model_rv1(x_phase + t0), model.model_rv2(x_phase + t0)], map_soln
)
plt.plot(x_phase, y1_mod - mean, "C0")
plt.plot(x_phase, y2_mod - mean, "C1")
plt.legend(fontsize=10)
plt.xlim(-0.5 * period, 0.5 * period)
plt.ylabel("radial velocity [km / s]")
plt.xlabel("time since primary eclipse [days]")
_ = plt.title("HD 23642; map model", fontsize=14)
And then the light curve. In the top panel, we show the Gaussian process model for the phase curve. It’s clear that there’s a lot of information there that we could take advantage of, but that’s a topic for another day. In the bottom panel, we’re plotting the phase folded light curve and we can see the ridiculous signal to noise that we’re getting on the eclipses.
with model:
gp_pred = xo.eval_in_model(model.gp_lc.predict(), map_soln) + map_soln["mean_lc"]
lc = xo.eval_in_model(model.model_lc(model.x), map_soln) - map_soln["mean_lc"]
fig, (ax1, ax2) = plt.subplots(2, sharex=True, figsize=(12, 7))
ax1.plot(model.x, model.y, "k.", alpha=0.2)
ax1.plot(model.x, gp_pred, color="C1", lw=1)
ax2.plot(model.x, model.y - gp_pred, "k.", alpha=0.2)
ax2.plot(model.x, lc, color="C2", lw=1)
ax2.set_xlim(model.x.min(), model.x.max())
ax1.set_ylabel("raw flux [ppt]")
ax2.set_ylabel("de-trended flux [ppt]")
ax2.set_xlabel("time [KBJD]")
ax1.set_title("HD 23642; map model", fontsize=14)
fig.subplots_adjust(hspace=0.05)
fig, ax1 = plt.subplots(1, figsize=(12, 3.5))
x_fold = (model.x - map_soln["t0"]) % map_soln["period"] / map_soln["period"]
inds = np.argsort(x_fold)
ax1.plot(x_fold[inds], model.y[inds] - gp_pred[inds], "k.", alpha=0.2)
ax1.plot(x_fold[inds] - 1, model.y[inds] - gp_pred[inds], "k.", alpha=0.2)
ax2.plot(x_fold[inds], model.y[inds] - gp_pred[inds], "k.", alpha=0.2, label="data!")
ax2.plot(x_fold[inds] - 1, model.y[inds] - gp_pred, "k.", alpha=0.2)
yval = model.y[inds] - gp_pred
bins = np.linspace(0, 1, 75)
num, _ = np.histogram(x_fold[inds], bins, weights=yval)
denom, _ = np.histogram(x_fold[inds], bins)
ax2.plot(0.5 * (bins[:-1] + bins[1:]) - 1, num / denom, ".w")
args = dict(lw=1)
ax1.plot(x_fold[inds], lc[inds], "C2", **args)
ax1.plot(x_fold[inds] - 1, lc[inds], "C2", **args)
ax1.set_xlim(-1, 1)
ax1.set_ylabel("de-trended flux [ppt]")
ax1.set_xlabel("phase")
_ = ax1.set_title("HD 23642; map model", fontsize=14)
Finally we can run the MCMC:
np.random.seed(23642)
with model:
trace = pm.sample(
tune=3500,
draws=3000,
start=map_soln,
chains=4,
step=xo.get_dense_nuts_step(start=map_soln, target_accept=0.95),
)
Multiprocess sampling (4 chains in 4 jobs) NUTS: [ell_rv, S_tot_rv, sigma_rv2, sigma_rv1, ell_lc, S_tot_lc, sigma_lc, ecs, t0, period, b, s, q, k, R1, M1, u2, u1, mean_rv, mean_lc] Sampling 4 chains, 6 divergences: 100%|██████████| 26000/26000 [41:28<00:00, 10.45draws/s] There were 6 divergences after tuning. Increase target_accept or reparameterize. The number of effective samples is smaller than 25% for some parameters.
As usual, we can check the convergence diagnostics for some of the key parameters.
pm.summary(trace, var_names=["M1", "M2", "R1", "R2", "ecs", "incl", "s"])
mean | sd | hpd_3% | hpd_97% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
M1 | 2.237 | 0.030 | 2.179 | 2.292 | 0.000 | 0.000 | 8331.0 | 8331.0 | 8762.0 | 7520.0 | 1.0 |
M2 | 1.564 | 0.023 | 1.520 | 1.608 | 0.000 | 0.000 | 8737.0 | 8737.0 | 9346.0 | 6936.0 | 1.0 |
R1 | 1.828 | 0.052 | 1.730 | 1.925 | 0.001 | 0.001 | 3153.0 | 3125.0 | 3123.0 | 2677.0 | 1.0 |
R2 | 1.384 | 0.095 | 1.208 | 1.564 | 0.002 | 0.001 | 2782.0 | 2782.0 | 2900.0 | 2438.0 | 1.0 |
ecs[0] | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 15228.0 | 14571.0 | 15204.0 | 8731.0 | 1.0 |
ecs[1] | -0.003 | 0.005 | -0.012 | 0.006 | 0.000 | 0.000 | 7756.0 | 6066.0 | 7829.0 | 7870.0 | 1.0 |
incl | 1.371 | 0.005 | 1.363 | 1.380 | 0.000 | 0.000 | 2529.0 | 2522.0 | 2921.0 | 2371.0 | 1.0 |
s | 0.471 | 0.020 | 0.434 | 0.509 | 0.000 | 0.000 | 6753.0 | 6753.0 | 6733.0 | 7604.0 | 1.0 |
It can be useful to take a look at some diagnostic corner plots to see how the sampling went. First, let’s look at some observables:
import corner
samples = pm.trace_to_dataframe(trace, varnames=["k", "q", "ecs"])
_ = corner.corner(
samples,
labels=["$k = R_2 / R_1$", "$q = M_2 / M_1$", "$e\,\cos\omega$", "$e\,\sin\omega$"],
)
And then we can look at the physical properties of the stars in the system. In this figure, we’re comparing to the results from David+ (2016). We find that the radius of the primary is marginally inconsistent with the David+ (2016) measurement, but more consistent with other estimates from the literature (that are overplotted in different colors).
samples = pm.trace_to_dataframe(trace, varnames=["R1", "R2", "M1", "M2"])
fig = corner.corner(samples, truths=[1.727, 1.503, 2.203, 1.5488])
for ax in np.reshape(fig.axes, (4, 4))[:, 0]:
ax.axvline(1.727, color="C0", label="David+ (2016)")
ax.axvline(1.890, color="C1", label="Groenewegen+ (2007)")
ax.axvline(1.831, color="C2", label="Southworth+ (2005)")
ax.axvline(1.81, color="C3", label="Munari+ (2004)")
_ = fig.axes[0].legend(fontsize=10, loc="center left", bbox_to_anchor=(1.1, 0.5))
As described in the Citing exoplanet & its dependencies tutorial, we can use
exoplanet.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} and its dependencies citep{exoplanet:agol19, exoplanet:astropy13, exoplanet:astropy18, exoplanet:exoplanet, exoplanet:foremanmackey17, exoplanet:foremanmackey18, exoplanet:kipping13, exoplanet:luger18, exoplanet:pymc3, exoplanet:theano}.
print("\n".join(bib.splitlines()[:10]) + "\n...")
@misc{exoplanet:exoplanet,
author = {Daniel Foreman-Mackey and Ian Czekala and Rodrigo Luger and
Eric Agol and Geert Barentsen and Tom Barclay},
title = {exoplanet-dev/exoplanet v0.3.0},
month = apr,
year = 2020,
doi = {10.5281/zenodo.1998447},
url = {https://doi.org/10.5281/zenodo.1998447}
}
...