Planets of every kind¶
This example demonstrates kima's ability to constrain the orbital parameters of different "types" of planets: RV-detected planets, known objects, and transiting planets. We will create a (contrived) simulated dataset to showcase the combination of all these components.
The "normal", known object, and transiting planets have their specific use cases. The mishmash shown in this example is not necessarily something we would do in a real analysis.
Running the example¶
This example can be used directly as
from kima.examples import every_kind_of_planet
model, res = every_kind_of_planet(run=True, load=True)
but here we'll go through it in more detail.
Simulating a dataset¶
First, we get some standard imports out of the way
import numpy as np
import matplotlib.pyplot as plt
import kima
We'll create a function to generate an RV dataset containing three planet signals:
we assume the first one is transiting, and define the time of periastron and then calculate the time of conjunction
for the second planet, we have information about the orbital period
finally, we also include a third short-period planet for which we will assume no extra information
To make things interesting, we remove a few data points so that the orbits are slightly harder to constrain.
def create_data():
np.random.seed(24) # to be reproducible
# random times and uncertainties
t = np.sort(np.random.uniform(0, 100, 77))
err = np.random.uniform(0.1, 0.3, t.size)
# the transiting planet
P1, K1, e1, w1 = 30, 2, 0.6, 0.1
Tp1 = 15 # define time of periastron
M0_epoch1 = t[0]
M01 = 2 * np.pi * (M0_epoch1 - Tp1) / P1
v = np.array(kima.keplerian(t, P1, K1, e1, w1, M01, M0_epoch1))
# calculate time of transit
f = np.pi/2 - w1
E = 2 * np.arctan(np.tan(f/2) * np.sqrt((1-e1)/(1+e1)))
Tc = Tp1 + P1/(2*np.pi) * (E - e1*np.sin(E))
# the known object
P2, K2, e2, w2, M02 = 45, 1.3, 0.1, 0.1, 1.0
M0_epoch2 = t[0]
v += np.array(kima.keplerian(t, P2, K2, e2, w2, M02, M0_epoch2))
# yet another planet
P3, K3, e3, w3, M03 = 2.1, 1.0, 0.0, 0.3, 0.0
M0_epoch3 = t[0]
v += np.array(kima.keplerian(t, P3, K3, e3, w3, M03, M0_epoch3))
# add some small Gaussian noise
v += np.random.normal(loc=0.0, scale=0.05, size=t.size)
# remove some points to make things interesting
t = t[v < 2]
err = err[v < 2]
v = v[v < 2]
np.savetxt('data.txt', np.c_[t, v, err], fmt='%10.5f')
return (
(t, v, err),
(P1, K1, e1, w1, M01, M0_epoch1, Tc),
(P2, K2, e2, w2, M02, M0_epoch2, None),
(P3, K3, e3, w3, M03, M0_epoch3, None),
)
(t, v, err), par1, par2, par3 = create_data()
Let's plot the data together with the true Keplerian model
tt = np.linspace(t[0], t[-1], 1000)
vv = np.c_[
kima.keplerian(tt, *par1[:-1]),
kima.keplerian(tt, *par2[:-1]),
kima.keplerian(tt, *par3[:-1]),
]
fig, axs = plt.subplots(2, 1)
axs[0].plot(tt, vv, label=[f'planet {i+1}' for i in range(3)])
axs[0].legend()
axs[0].set(xlabel='Time [days]', ylabel='RV [m/s]');
axs[1].errorbar(t, v, err, fmt='o', label='data')
axs[1].plot(tt, vv.sum(axis=1), color='k', label='planet signals')
axs[1].legend()
axs[1].set(xlabel='Time [days]', ylabel='RV [m/s]');
Fitting all planets simultaneously¶
Let's now fit these data with a model that accounts for the different planet types.
from kima import RVData, RVmodel, distributions
First, we load the data
data = RVData(t, v, err, instrument='data')
In the baseline model, we will allow for up to 2 planets, setting the number of "normal" Keplerians free. For these Keplerians, we'll use all the default priors.
model = RVmodel(fix=False, npmax=2, data=data)
we include the transiting planet, setting an informative prior for the time of conjunction
Tc = par1[-1]
model.set_transiting_planet(1)
model.TR_Pprior = [distributions.LogUniform(1.0, data.get_timespan())]
model.TR_Kprior = [distributions.Uniform(0.0, 5*data.get_RV_span())]
model.TR_eprior = [distributions.Uniform(0, 1)]
model.TR_Tcprior = [distributions.Gaussian(Tc, 1e-3)]
model.TR_wprior = [distributions.Uniform(0, 2*np.pi)]
and a known object, with an informative prior for the orbital period as well as slightly informative priors for the semi-amplitude and eccentricity
P = par2[0]
model.set_known_object(1)
model.KO_Pprior = [distributions.Gaussian(P, 0.3)]
model.KO_Kprior = [distributions.Uniform(0, 2)]
model.KO_eprior = [distributions.Uniform(0, 0.2)]
model.KO_phiprior = [distributions.UniformAngle()]
model.KO_wprior = [distributions.UniformAngle()]
Finally, we fix the jitter parameter to zero (since we didn't include any in the simulated data) and we enforce AMD stability of the planet system
model.Jprior = distributions.Fixed(0.0)
model.enforce_stability = True
Enforcing AMD stability is ...
Run the model for a few thousand steps
kima.run(model, steps=1000, num_threads=8)
and load the results
res = kima.load_results(model)
log(Z) = -585.3736357283343 Information = 15.853541227368396 nats. Effective sample size = 1.021353814925894
100%|██████████| 1/1 [00:00<00:00, 215.27it/s]
Results¶
Plotting a few posterior samples together with the data already shows that we didn't recover the true Keplerian signal very well
fig = res.plot_random_samples(
isolate_known_object=False,
isolate_transiting_planet=False,
)
res.plot2();
and indeed the posteriors for the orbital parameters show that, while the period was recovered, the semi-amplitude and eccentricity are a bit off from the true values
res1.corner_planet_parameters(true_values=[P, K, e, M0, w],
wrap_M0=True, wrap_w=True);
Assuming $T_c$ is known¶
In our convoluted example, we do have information about the time of conjunction $T_c$. If we assume we don't know the orbital period, then we are basically in a similar situation to a mono-transit or a duo-transit.
Let's fit the data again with a model that uses one transiting planet, allowing us to assign a prior to $T_c$ directly.
#same as before
data = RVData(t, v, err)
# no 'normal' Keplerians, as we'll use a 'transiting planet'
model = RVmodel(fix=True, npmax=0, data=data)
# same as before
model.Jprior = distributions.Fixed(0.0)
# assume there is one planet
model.set_transiting_planet(1)
# set standard priors (same as before)
model.TR_Pprior = [distributions.LogUniform(1.0, data.get_timespan())]
model.TR_Kprior = [distributions.Uniform(0.0, 5*data.get_RV_span())]
model.TR_eprior = [distributions.Uniform(0, 1)]
model.TR_wprior = [distributions.UniformAngle()]
# but for Tc, use the known value
model.TR_Tcprior = [distributions.Gaussian(Tc, 1e-3)]
Let's run the model again
kima.run(model, steps=30_000, num_threads=4)
res2 = kima.load_results(model)
log(Z) = 15.826318304580681 Information = 15.872243383699903 nats. Effective sample size = 4353.03798580868
100%|██████████| 4353/4353 [00:00<00:00, 32585.88it/s]
and now the true Keplerian signal is much better recovered
fig = res2.plot_random_samples(legend=False, isolate_transiting_planet=False)
fig.axes[0].plot(tt, vv, color='g');
The posteriors for the orbital parameters (of the transiting planet) show that just by setting a prior on $T_c$, we are now able to recover the true values of the semi-amplitude and eccentricity
res2.corner_planet_parameters(true_values=[P, K, e, Tc, w],
include_transiting_planet=True,
wrap_w=True);