Transit Fitting¶

Welcome to the transit fitting tutorial of chromatic_fitting! In this tutorial we will go through how to create a simulated transit using the handy chromatic package and then fit it using the transit model (TransitModel) in chromatic_fitting.

In [1]:

# import chromatic_fitting of course!
from chromatic_fitting import *

# import any prior distributions we want to use for our parameters - I've chosen Normal and Uniform from pymc3
# and QuadLimbDark and ImpactParameter from exoplanet
from pymc3 import Normal, Uniform
from exoplanet import QuadLimbDark, ImpactParameter

# import astropy table just for this tutorial (pandas dataframes don't show up nicely on the docs pages!)
from astropy.table import Table

plt.matplotlib.style.use('default')

# import chromatic_fitting of course! from chromatic_fitting import * # import any prior distributions we want to use for our parameters - I've chosen Normal and Uniform from pymc3 # and QuadLimbDark and ImpactParameter from exoplanet from pymc3 import Normal, Uniform from exoplanet import QuadLimbDark, ImpactParameter # import astropy table just for this tutorial (pandas dataframes don't show up nicely on the docs pages!) from astropy.table import Table plt.matplotlib.style.use('default')

Running chromatic_fitting v0.0.2!

This program is running on:
Python v3.9.12 (main, Jun  1 2022, 06:34:44) 
[Clang 12.0.0 ]
numpy v1.22.1
chromatic v0.3.14
pymc3 v3.11.5
pymc3_ext v0.1.1
exoplanet v0.5.2

Create Synthetic Rainbow + Transit¶

To create our simulated data set we will use SimulatedRainbow() from within chromatic. This creates only basic data with time and wavelength axes. Then we want to inject some noise (to make it realistic) and then inject a transit (the whole point of this tutorial). By default chromatic will inject a planet with Rp/R* of 0.1, but we want to make things a little more realistic so we will make the radius ratio vary linearly with wavelength.

In [2]:

# create transit rainbow:
r = SimulatedRainbow(dt=1 * u.minute, R=50).inject_noise(signal_to_noise=100)

# add transit (with depth varying with wavelength):
r = r.inject_transit(
        planet_radius=np.linspace(0.2, 0.15, r.nwave))

# bin into 5 wavelength bins:
rbin5 = r.bin(nwavelengths=int(r.nwave/5), dt=5 * u.minute)

# show the simulated dataset
rbin5.imshow_quantities();

# create transit rainbow: r = SimulatedRainbow(dt=1 * u.minute, R=50).inject_noise(signal_to_noise=100) # add transit (with depth varying with wavelength): r = r.inject_transit( planet_radius=np.linspace(0.2, 0.15, r.nwave)) # bin into 5 wavelength bins: rbin5 = r.bin(nwavelengths=int(r.nwave/5), dt=5 * u.minute) # show the simulated dataset rbin5.imshow_quantities();

  0%|          | 0/116 [00:00<?, ?it/s]

  0%|          | 0/61 [00:00<?, ?it/s]

Define a PyMC3 Transit Model¶

Now onto chromatic_fitting. First we need to create a transit model.

In [3]:

# create transit model:
t = TransitModel()

# create transit model: t = TransitModel()

Then we will want to decide how our parameters will vary. We can see the parameters that we will need to provide to the model. Warning: if we dont set any of these parameters then they will take on default values.

In [4]:

t.required_parameters

t.required_parameters

Out[4]:

['stellar_radius',
 'stellar_mass',
 'radius_ratio',
 'period',
 'epoch',
 'baseline',
 'impact_parameter',
 'limb_darkening']

In [5]:

t.defaults

t.defaults

Out[5]:

{'stellar_radius': 1.0,
 'stellar_mass': 1.0,
 'radius_ratio': 1.0,
 'period': 1.0,
 'epoch': 0.0,
 'baseline': 1.0,
 'impact_parameter': 0.5,
 'eccentricity': 0.0,
 'omega': 1.5707963267948966,
 'limb_darkening': [0.2, 0.2]}

Usually when we're fitting a transit we'll have some idea about the transit parameters (from previous literature or just looking at the lightcurves by eye), so it's a good idea to give good initial estimates to help our sampling converge nicely. Here we're defining the prior distributions for each parameter. There are four options for parameters in chromatic_fitting: Fixed, WavelikeFixed, Fitted, and WavelikeFitted. Fixed is one value fixed across all wavelengths, WavelikeFixed are fixed values that can vary between wavelengths. Fitted determines a prior distribution (e.g. Uniform, Normal, TruncatedNormal) that we will use to fit one value for the parameter across all wavelengths. Similarly, WavelikeFitted is a prior distribution that we will use to fit for a different value for every wavelength.

In [6]:

# add our parameters:
t.setup_parameters(
                  period=1, # a fixed value!
                   epoch=Fitted(Uniform,lower=-0.05,upper=0.05), # one fitted value across all wavelengths
                   stellar_radius = Fitted(Uniform, lower=0.8, upper=1.2,testval=1),
                   stellar_mass =Fitted(Uniform, lower=0.8, upper=1.2,testval=1),
                   radius_ratio=WavelikeFitted(Normal, mu=0.5, sigma=0.05), # a different value fitted for every wavelength!
                   impact_parameter=Fitted(ImpactParameter,ror=0.15,testval=0.44),
                   limb_darkening=WavelikeFitted(QuadLimbDark,testval=[0.05,0.35]),
                    baseline = WavelikeFitted(Normal, mu=1.0, sigma=0.05), 
                )

# print a summary of all params:
t.summarize_parameters()

# add our parameters: t.setup_parameters( period=1, # a fixed value! epoch=Fitted(Uniform,lower=-0.05,upper=0.05), # one fitted value across all wavelengths stellar_radius = Fitted(Uniform, lower=0.8, upper=1.2,testval=1), stellar_mass =Fitted(Uniform, lower=0.8, upper=1.2,testval=1), radius_ratio=WavelikeFitted(Normal, mu=0.5, sigma=0.05), # a different value fitted for every wavelength! impact_parameter=Fitted(ImpactParameter,ror=0.15,testval=0.44), limb_darkening=WavelikeFitted(QuadLimbDark,testval=[0.05,0.35]), baseline = WavelikeFitted(Normal, mu=1.0, sigma=0.05), ) # print a summary of all params: t.summarize_parameters()

transit_stellar_radius =
  <🧮 Fitted Uniform(lower=0.8, upper=1.2, testval=1, name='transit_stellar_radius') 🧮>

transit_stellar_mass =
  <🧮 Fitted Uniform(lower=0.8, upper=1.2, testval=1, name='transit_stellar_mass') 🧮>

transit_radius_ratio =
  <🧮 WavelikeFitted Normal(mu=0.5, sigma=0.05, name='transit_radius_ratio') for each wavelength 🧮>

transit_period =
  <🧮 Fixed | 1 🧮>

transit_epoch =
  <🧮 Fitted Uniform(lower=-0.05, upper=0.05, name='transit_epoch') 🧮>

transit_baseline =
  <🧮 WavelikeFitted Normal(mu=1.0, sigma=0.05, name='transit_baseline') for each wavelength 🧮>

transit_impact_parameter =
  <🧮 Fitted ImpactParameter(ror=0.15, testval=0.44, name='transit_impact_parameter') 🧮>

transit_eccentricity =
  <🧮 Fixed | 0.0 🧮>

transit_omega =
  <🧮 Fixed | 1.5707963267948966 🧮>

transit_limb_darkening =
  <🧮 WavelikeFitted QuadLimbDark(testval=[0.05, 0.35], name='transit_limb_darkening') for each wavelength 🧮>

Attach the Rainbow Object and Set-up the Model¶

The next step is to attach the actual data to the model and setup the lightcurves!

In [7]:

# attach the Rainbow object to the model:
t.attach_data(rbin5)

# setup the lightcurves for the transit model:
t.setup_lightcurves()

# relate the "actual" data to the model (using a Normal likelihood function)
t.setup_likelihood()

# attach the Rainbow object to the model: t.attach_data(rbin5) # setup the lightcurves for the transit model: t.setup_lightcurves() # relate the "actual" data to the model (using a Normal likelihood function) t.setup_likelihood()

  0%|          | 0/5 [00:00<?, ?it/s]

If we look at our PyMC3 model we can see that it has a lot of parameters to optimize! If we've chosen the separate wavelength fitting method (.choose_optimization_method("separate")) then ._pymc3_model will return a list of PyMC3 models (one for each wavelength).

In [8]:

print(t._pymc3_model)

print(t._pymc3_model)

                transit_epoch_interval__ ~ TransformedDistribution
       transit_impact_parameter_impact__ ~ TransformedDistribution
       transit_stellar_radius_interval__ ~ TransformedDistribution
         transit_stellar_mass_interval__ ~ TransformedDistribution
transit_limb_darkening_w0_quadlimbdark__ ~ TransformedDistribution
                 transit_radius_ratio_w0 ~ Normal
                     transit_baseline_w0 ~ Normal
transit_limb_darkening_w1_quadlimbdark__ ~ TransformedDistribution
                 transit_radius_ratio_w1 ~ Normal
                     transit_baseline_w1 ~ Normal
transit_limb_darkening_w2_quadlimbdark__ ~ TransformedDistribution
                 transit_radius_ratio_w2 ~ Normal
                     transit_baseline_w2 ~ Normal
transit_limb_darkening_w3_quadlimbdark__ ~ TransformedDistribution
                 transit_radius_ratio_w3 ~ Normal
                     transit_baseline_w3 ~ Normal
transit_limb_darkening_w4_quadlimbdark__ ~ TransformedDistribution
                 transit_radius_ratio_w4 ~ Normal
                     transit_baseline_w4 ~ Normal
                           transit_epoch ~ Uniform
                transit_impact_parameter ~ ImpactParameter
                  transit_stellar_radius ~ Uniform
                    transit_stellar_mass ~ Uniform
                            transit_a_R* ~ Deterministic
               transit_limb_darkening_w0 ~ QuadLimbDark
               transit_limb_darkening_w1 ~ QuadLimbDark
               transit_limb_darkening_w2 ~ QuadLimbDark
               transit_limb_darkening_w3 ~ QuadLimbDark
               transit_limb_darkening_w4 ~ QuadLimbDark
                       wavelength_0_data ~ Normal
                       wavelength_1_data ~ Normal
                       wavelength_2_data ~ Normal
                       wavelength_3_data ~ Normal
                       wavelength_4_data ~ Normal

We've got our transit parameters that are wavelength-independant in this case (epoch, impact_parameter, stellar_radius, stellar_mass, a_R*, limb_darkening) and the wavelength-dependant parameters we've defined to be radius_ratio_w{N} and baseline_w{N} only. The wavelength_{N}_data parameter just represents the fit of the data to the model at each wavelength (which we've defined to be a Normal distribution). If we've set store_models = True at the .setup_lightcurves() stage then we will also see a bunch of models!

Now we can plot a couple of priors (samples from our prior distribution) - do they look OK?

In [9]:

t.plot_priors()

t.plot_priors()

/Users/camu5866/opt/anaconda3/envs/chromaticfitting/lib/python3.9/site-packages/exoplanet/distributions/transforms.py:35: RuntimeWarning: invalid value encountered in log
  return np.log(q) - np.log(1 - q)
/Users/camu5866/opt/anaconda3/envs/chromaticfitting/lib/python3.9/site-packages/exoplanet/distributions/transforms.py:35: RuntimeWarning: invalid value encountered in log
  return np.log(q) - np.log(1 - q)

Plotting the priors can be reassuring for two reasons: (1) we're not giving priors that are vastly off from the true values, and (2) we're not over-constraining our model by giving it the exact solution and priors that are too tight.

Another check is does this planetary system actually transit given our parameters?

In [10]:

t.plot_orbit()

t.plot_orbit()

Looks good! And a final check of what the actual lightcurves look like:

In [11]:

t.plot_lightcurves()

t.plot_lightcurves()

The summarize step has not been run yet. To include the 'best-fit' model please run {self}.summarize() before calling this step!

PyMC3 sampling¶

Now we can run the NUTS sampling for our light curves (first by optimizing our initial values)

In [12]:

# optimize for initial values!
opt = t.optimize(plot=False)

# put those initial values into the sampling and define the number of tuning and draw steps, 
# as well as the number of chains. NOTE: if you do separate wavelength fitting then the number of steps 
# is per wavelengths, not divided between the wavelengths!
t.sample(start=opt, tune=4000, draws=4000, chains=4, cores=4)

# optimize for initial values! opt = t.optimize(plot=False) # put those initial values into the sampling and define the number of tuning and draw steps, # as well as the number of chains. NOTE: if you do separate wavelength fitting then the number of steps # is per wavelengths, not divided between the wavelengths! t.sample(start=opt, tune=4000, draws=4000, chains=4, cores=4)

optimizing logp for variables: [transit_baseline_w4, transit_radius_ratio_w4, transit_limb_darkening_w4, transit_baseline_w3, transit_radius_ratio_w3, transit_limb_darkening_w3, transit_baseline_w2, transit_radius_ratio_w2, transit_limb_darkening_w2, transit_baseline_w1, transit_radius_ratio_w1, transit_limb_darkening_w1, transit_baseline_w0, transit_radius_ratio_w0, transit_limb_darkening_w0, transit_stellar_mass, transit_stellar_radius, transit_impact_parameter, transit_epoch]

100.00% [27/27 00:00<00:00 logp = nan]

message: Desired error not necessarily achieved due to precision loss.
logp: -1976751.023912184 -> -197432.9221621303

/Users/camu5866/opt/anaconda3/envs/chromaticfitting/lib/python3.9/site-packages/deprecat/classic.py:215: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.
  return wrapped_(*args_, **kwargs_)
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [transit_baseline_w4, transit_radius_ratio_w4, transit_limb_darkening_w4, transit_baseline_w3, transit_radius_ratio_w3, transit_limb_darkening_w3, transit_baseline_w2, transit_radius_ratio_w2, transit_limb_darkening_w2, transit_baseline_w1, transit_radius_ratio_w1, transit_limb_darkening_w1, transit_baseline_w0, transit_radius_ratio_w0, transit_limb_darkening_w0, transit_stellar_mass, transit_stellar_radius, transit_impact_parameter, transit_epoch]

100.00% [32000/32000 04:20<00:00 Sampling 4 chains, 238 divergences]

/Users/camu5866/opt/anaconda3/envs/chromaticfitting/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf
  return _boost._beta_ppf(q, a, b)
/Users/camu5866/opt/anaconda3/envs/chromaticfitting/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf
  return _boost._beta_ppf(q, a, b)
/Users/camu5866/opt/anaconda3/envs/chromaticfitting/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf
  return _boost._beta_ppf(q, a, b)
/Users/camu5866/opt/anaconda3/envs/chromaticfitting/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf
  return _boost._beta_ppf(q, a, b)
Sampling 4 chains for 4_000 tune and 4_000 draw iterations (16_000 + 16_000 draws total) took 261 seconds.
There were 21 divergences after tuning. Increase `target_accept` or reparameterize.
There were 78 divergences after tuning. Increase `target_accept` or reparameterize.
There were 138 divergences after tuning. Increase `target_accept` or reparameterize.
There was 1 divergence after tuning. Increase `target_accept` or reparameterize.
The number of effective samples is smaller than 25% for some parameters.

At this stage the sampler may print out some warnings that we don't have enough tuning steps! Also be aware this step can be quite slow depending on how many parameters we're trying to fit at once (which may increase with the number of wavelengths). We can then see the results of our sampling by running .summarize():

In [13]:

t.summarize(round_to=7, hdi_prob=0.68, fmt='wide')

t.summarize(round_to=7, hdi_prob=0.68, fmt='wide')

                                  mean        sd   hdi_16%   hdi_84%  \
transit_radius_ratio_w0       0.196850  0.001054  0.195735  0.197820   
transit_baseline_w0           1.000135  0.000177  0.999955  1.000308   
transit_radius_ratio_w1       0.185004  0.001068  0.183938  0.186037   
transit_baseline_w1           0.999951  0.000176  0.999782  1.000131   
transit_radius_ratio_w2       0.175489  0.001032  0.174520  0.176575   
transit_baseline_w2           0.999971  0.000176  0.999798  1.000150   
transit_radius_ratio_w3       0.167636  0.001032  0.166691  0.168739   
transit_baseline_w3           1.000046  0.000169  0.999875  1.000210   
transit_radius_ratio_w4       0.157092  0.001048  0.156066  0.158183   
transit_baseline_w4           1.000005  0.000170  0.999833  1.000176   
transit_epoch                 0.000010  0.000073 -0.000063  0.000081   
transit_impact_parameter      0.072298  0.048729  0.000032  0.092909   
transit_stellar_radius        1.179374  0.012408  1.170083  1.197520   
transit_stellar_mass          0.844964  0.026998  0.805358  0.865231   
transit_a_R*                  3.373038  0.016508  3.360682  3.391459   
transit_limb_darkening_w0[0]  0.282839  0.087861  0.200358  0.377082   
transit_limb_darkening_w0[1]  0.144619  0.158613 -0.017603  0.301598   
transit_limb_darkening_w1[0]  0.302008  0.095797  0.206296  0.396989   
transit_limb_darkening_w1[1]  0.206545  0.171916  0.026470  0.369993   
transit_limb_darkening_w2[0]  0.419366  0.085170  0.350372  0.517058   
transit_limb_darkening_w2[1] -0.028803  0.146806 -0.226166  0.062104   
transit_limb_darkening_w3[0]  0.331411  0.087270  0.262827  0.434720   
transit_limb_darkening_w3[1] -0.012111  0.144503 -0.205404  0.065244   
transit_limb_darkening_w4[0]  0.399725  0.091249  0.331058  0.506786   
transit_limb_darkening_w4[1] -0.052149  0.144996 -0.243752  0.020698   

                                 mcse_mean       mcse_sd      ess_bulk  \
transit_radius_ratio_w0       1.130000e-05  8.000000e-06   8964.806137   
transit_baseline_w0           1.400000e-06  1.000000e-06  15602.497523   
transit_radius_ratio_w1       1.030000e-05  7.300000e-06  10817.294878   
transit_baseline_w1           1.500000e-06  1.000000e-06  14688.887443   
transit_radius_ratio_w2       1.030000e-05  7.300000e-06   9979.889853   
transit_baseline_w2           1.300000e-06  1.000000e-06  17114.559639   
transit_radius_ratio_w3       9.400000e-06  6.600000e-06  11976.130095   
transit_baseline_w3           1.300000e-06  9.000000e-07  16035.331257   
transit_radius_ratio_w4       9.900000e-06  7.000000e-06  11204.001252   
transit_baseline_w4           1.400000e-06  1.000000e-06  14605.592146   
transit_epoch                 5.000000e-07  6.000000e-07  19301.645809   
transit_impact_parameter      6.890000e-04  5.826000e-04   4621.771437   
transit_stellar_radius        1.807000e-04  1.278000e-04   4030.901344   
transit_stellar_mass          3.962000e-04  2.802000e-04   3502.906550   
transit_a_R*                  2.309000e-04  1.633000e-04   6187.964172   
transit_limb_darkening_w0[0]  9.613000e-04  7.369000e-04   8031.357267   
transit_limb_darkening_w0[1]  1.784500e-03  1.261900e-03   7314.561667   
transit_limb_darkening_w1[0]  1.014000e-03  7.171000e-04   8742.894845   
transit_limb_darkening_w1[1]  1.842500e-03  1.460100e-03   8375.544704   
transit_limb_darkening_w2[0]  9.557000e-04  7.188000e-04   7172.844669   
transit_limb_darkening_w2[1]  1.821600e-03  1.288100e-03   5241.889446   
transit_limb_darkening_w3[0]  9.236000e-04  6.855000e-04   8326.935897   
transit_limb_darkening_w3[1]  1.662000e-03  1.175300e-03   5802.919728   
transit_limb_darkening_w4[0]  9.346000e-04  7.048000e-04   8545.725944   
transit_limb_darkening_w4[1]  1.607700e-03  1.209700e-03   5851.202906   

                                  ess_tail     r_hat  
transit_radius_ratio_w0        5963.959924  1.000038  
transit_baseline_w0           10103.415532  1.000604  
transit_radius_ratio_w1        6507.848928  1.000541  
transit_baseline_w1            8830.209029  1.000364  
transit_radius_ratio_w2        7769.737289  1.000677  
transit_baseline_w2           11270.736585  1.001561  
transit_radius_ratio_w3        8664.438824  1.000265  
transit_baseline_w3            9692.552933  1.000007  
transit_radius_ratio_w4        8520.385442  1.000551  
transit_baseline_w4           10234.023422  1.000196  
transit_epoch                 10860.905012  1.000366  
transit_impact_parameter       3419.281968  1.000693  
transit_stellar_radius         3579.652342  1.001794  
transit_stellar_mass           1762.094829  1.001168  
transit_a_R*                   3772.370721  1.000964  
transit_limb_darkening_w0[0]   4694.432452  1.000342  
transit_limb_darkening_w0[1]   3786.230563  1.000531  
transit_limb_darkening_w1[0]   5336.131956  1.000500  
transit_limb_darkening_w1[1]   4727.750183  1.000340  
transit_limb_darkening_w2[0]   7165.335530  1.000304  
transit_limb_darkening_w2[1]   3300.049279  1.000720  
transit_limb_darkening_w3[0]   7883.285443  1.000342  
transit_limb_darkening_w3[1]   3557.528409  1.000540  
transit_limb_darkening_w4[0]   7320.797781  1.000565  
transit_limb_darkening_w4[1]   4054.146420  1.000761  

An important parameter to look out for here to check whether your samplings have converged is r_hat. This rank normalized R-hat diagnostic tests for lack of convergence by comparing the variance between multiple chains to the variance within each chain. If convergence has been achieved, the between-chain and within-chain variances should be identicalm therefore, the closer it is to 1 the better the chance that your sampling successfully converged. If you're interested in the sampling see the PyMC3 docs for much more detail!

We might also want to see a couple of posterior samples as a "quick-look" check! But beware if you've chosen the "separate" optimization method then it will plot the posteriors for every wavelength!):

In [14]:

# t.plot_posteriors()

# t.plot_posteriors()

But what are the results?? We can easily see the results using the handy .get_results() function:

In [15]:

results = t.get_results(uncertainty=['hdi_16%','hdi_84%'])

# results is a pandas dataframe, however, it doesn't show up properly on Git docs so I'll convert it to an
# astropy table
from astropy.table import Table
Table.from_pandas(results)

results = t.get_results(uncertainty=['hdi_16%','hdi_84%']) # results is a pandas dataframe, however, it doesn't show up properly on Git docs so I'll convert it to an # astropy table from astropy.table import Table Table.from_pandas(results)

Out[15]:

Table length=5

transit_baseline	transit_baseline_hdi_16%	transit_baseline_hdi_84%	transit_eccentricity	transit_eccentricity_hdi_16%	transit_eccentricity_hdi_84%	transit_epoch	transit_epoch_hdi_16%	transit_epoch_hdi_84%	transit_impact_parameter	transit_impact_parameter_hdi_16%	transit_impact_parameter_hdi_84%	transit_limb_darkening	transit_limb_darkening_hdi_16%	transit_limb_darkening_hdi_84%	transit_omega	transit_omega_hdi_16%	transit_omega_hdi_84%	transit_period	transit_period_hdi_16%	transit_period_hdi_84%	transit_radius_ratio	transit_radius_ratio_hdi_16%	transit_radius_ratio_hdi_84%	transit_stellar_mass	transit_stellar_mass_hdi_16%	transit_stellar_mass_hdi_84%	transit_stellar_radius	transit_stellar_radius_hdi_16%	transit_stellar_radius_hdi_84%	wavelength
object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object	object
1.0001349	0.9999549	1.0003083	0.0	0.0	0.0	1.02e-05	-6.31e-05	8.07e-05	0.072298	3.18e-05	0.0929086	[0.2828392, 0.1446192]	[0.2003576, -0.0176025]	[0.3770824, 0.3015977]	1.5707963267948966	1.5707963267948966	1.5707963267948966	1	1	1	0.1968495	0.1957345	0.1978196	0.8449638	0.8053577	0.8652308	1.1793739	1.1700833	1.1975196	0.639572482934883 micron
0.9999506	0.9997816	1.0001306	0.0	0.0	0.0	1.02e-05	-6.31e-05	8.07e-05	0.072298	3.18e-05	0.0929086	[0.3020083, 0.206545]	[0.2062965, 0.0264702]	[0.3969892, 0.369993]	1.5707963267948966	1.5707963267948966	1.5707963267948966	1	1	1	0.1850039	0.183938	0.1860373	0.8449638	0.8053577	0.8652308	1.1793739	1.1700833	1.1975196	1.013209338074884 micron
0.999971	0.9997975	1.0001496	0.0	0.0	0.0	1.02e-05	-6.31e-05	8.07e-05	0.072298	3.18e-05	0.0929086	[0.4193659, -0.0288029]	[0.3503719, -0.2261662]	[0.5170579, 0.0621045]	1.5707963267948966	1.5707963267948966	1.5707963267948966	1	1	1	0.1754891	0.1745199	0.1765749	0.8449638	0.8053577	0.8652308	1.1793739	1.1700833	1.1975196	1.604998553797903 micron
1.000046	0.9998747	1.0002097	0.0	0.0	0.0	1.02e-05	-6.31e-05	8.07e-05	0.072298	3.18e-05	0.0929086	[0.3314109, -0.0121108]	[0.2628265, -0.2054041]	[0.4347198, 0.0652441]	1.5707963267948966	1.5707963267948966	1.5707963267948966	1	1	1	0.1676364	0.1666909	0.1687389	0.8449638	0.8053577	0.8652308	1.1793739	1.1700833	1.1975196	2.542436455025025 micron
1.0000054	0.9998326	1.0001761	0.0	0.0	0.0	1.02e-05	-6.31e-05	8.07e-05	0.072298	3.18e-05	0.0929086	[0.3997247, -0.0521489]	[0.331058, -0.2437521]	[0.5067859, 0.0206982]	1.5707963267948966	1.5707963267948966	1.5707963267948966	1	1	1	0.1570924	0.1560659	0.1581833	0.8449638	0.8053577	0.8652308	1.1793739	1.1700833	1.1975196	4.027407446906737 micron

We can also make a transmission spectrum table using .make_transmission_spectrum_table()

In [16]:

# this will show errors for the models that aren't transit and return a list!
transmission_spectrum = t.make_transmission_spectrum_table(uncertainty=['hdi_16%','hdi_84%'])

Table.from_pandas(transmission_spectrum)

# this will show errors for the models that aren't transit and return a list! transmission_spectrum = t.make_transmission_spectrum_table(uncertainty=['hdi_16%','hdi_84%']) Table.from_pandas(transmission_spectrum)

/Users/camu5866/opt/anaconda3/envs/chromaticfitting/lib/python3.9/site-packages/chromatic_fitting/models/transit.py:366: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  trans_table[f"{self.name}_radius_ratio_neg_error"] = (

Out[16]:

Table length=5

wavelength	transit_radius_ratio	transit_radius_ratio_neg_error	transit_radius_ratio_pos_error
object	object	object	object
0.639572482934883 micron	0.1968495	0.0011150000000000049	0.0009701000000000015
1.013209338074884 micron	0.1850039	0.0010659000000000085	0.0010333999999999899
1.604998553797903 micron	0.1754891	0.0009692000000000034	0.0010857999999999979
2.542436455025025 micron	0.1676364	0.000945499999999988	0.0011025000000000063
4.027407446906737 micron	0.1570924	0.0010264999999999858	0.0010909000000000058

We can also create (or call if you set store_models=True) the final best-fit models using the .get_models() function. This can take a minute if generating the models from parameters for lots of wavelengths.

In [17]:

models = t.get_model()
models.keys()

models = t.get_model() models.keys()

Out[17]:

dict_keys(['w0', 'w1', 'w2', 'w3', 'w4'])

As the .get_model() process can take time to generate the models, we store the models for use later...

In [18]:

t._fit_models

t._fit_models

Out[18]:

{'w0': array([1.0001349 , 1.0001349 , 1.0001349 , 1.0001349 , 1.0001349 ,
        1.0001349 , 1.0001349 , 1.0001349 , 1.0001349 , 1.0001349 ,
        1.0001349 , 1.0001349 , 1.0001349 , 1.0001349 , 0.99865364,
        0.99336154, 0.9860587 , 0.97773319, 0.96944652, 0.96314242,
        0.9613083 , 0.96004426, 0.95908471, 0.95833581, 0.9577462 ,
        0.95728393, 0.95692778, 0.9566633 , 0.95648061, 0.95637332,
        0.95633781, 0.9563729 , 0.95647976, 0.95666199, 0.95692597,
        0.95728154, 0.95774314, 0.95833192, 0.95907976, 0.96003786,
        0.96129967, 0.96312328, 0.96940077, 0.97768318, 0.98601186,
        0.99332344, 0.99863094, 1.0001349 , 1.0001349 , 1.0001349 ,
        1.0001349 , 1.0001349 , 1.0001349 , 1.0001349 , 1.0001349 ,
        1.0001349 , 1.0001349 , 1.0001349 , 1.0001349 , 1.0001349 ,
        1.0001349 ]),
 'w1': array([0.9999506 , 0.9999506 , 0.9999506 , 0.9999506 , 0.9999506 ,
        0.9999506 , 0.9999506 , 0.9999506 , 0.9999506 , 0.9999506 ,
        0.9999506 , 0.9999506 , 0.9999506 , 0.9999506 , 0.99916651,
        0.99470675, 0.98810614, 0.98040434, 0.97271675, 0.96736782,
        0.96551083, 0.96421232, 0.96323561, 0.96248111, 0.96189299,
        0.96143608, 0.96108689, 0.96082929, 0.96065228, 0.9605487 ,
        0.96051448, 0.9605483 , 0.96065147, 0.96082802, 0.96108512,
        0.96143373, 0.96188995, 0.96247721, 0.96323059, 0.96420578,
        0.96550194, 0.96735409, 0.97267467, 0.98035784, 0.98806316,
        0.99467308, 0.99914936, 0.9999506 , 0.9999506 , 0.9999506 ,
        0.9999506 , 0.9999506 , 0.9999506 , 0.9999506 , 0.9999506 ,
        0.9999506 , 0.9999506 , 0.9999506 , 0.9999506 , 0.9999506 ,
        0.9999506 ]),
 'w2': array([0.999971  , 0.999971  , 0.999971  , 0.999971  , 0.999971  ,
        0.999971  , 0.999971  , 0.999971  , 0.999971  , 0.999971  ,
        0.999971  , 0.999971  , 0.999971  , 0.999971  , 0.99950008,
        0.99529593, 0.9889944 , 0.98174186, 0.97472388, 0.97065505,
        0.96930987, 0.96825175, 0.96738414, 0.96666445, 0.9660687 ,
        0.9655818 , 0.96519376, 0.96489777, 0.96468921, 0.9645651 ,
        0.96452375, 0.96456462, 0.96468824, 0.96489629, 0.96519176,
        0.96557924, 0.96606553, 0.96666061, 0.96737951, 0.96824615,
        0.96930296, 0.97064587, 0.97468673, 0.98169855, 0.98895357,
        0.99526369, 0.99948513, 0.999971  , 0.999971  , 0.999971  ,
        0.999971  , 0.999971  , 0.999971  , 0.999971  , 0.999971  ,
        0.999971  , 0.999971  , 0.999971  , 0.999971  , 0.999971  ,
        0.999971  ]),
 'w3': array([1.000046  , 1.000046  , 1.000046  , 1.000046  , 1.000046  ,
        1.000046  , 1.000046  , 1.000046  , 1.000046  , 1.000046  ,
        1.000046  , 1.000046  , 1.000046  , 1.000046  , 0.99981994,
        0.99587146, 0.98973456, 0.9827221 , 0.97613391, 0.97296487,
        0.97200729, 0.97125484, 0.97063981, 0.97013118, 0.96971125,
        0.9693688 , 0.96909637, 0.96888886, 0.96874279, 0.96865593,
        0.968627  , 0.96865559, 0.96874211, 0.96888782, 0.96909497,
        0.969367  , 0.96970902, 0.97012847, 0.97063653, 0.97125087,
        0.97200237, 0.97295836, 0.97610029, 0.98268062, 0.98969482,
        0.99584018, 0.99980796, 1.000046  , 1.000046  , 1.000046  ,
        1.000046  , 1.000046  , 1.000046  , 1.000046  , 1.000046  ,
        1.000046  , 1.000046  , 1.000046  , 1.000046  , 1.000046  ,
        1.000046  ]),
 'w4': array([1.0000054 , 1.0000054 , 1.0000054 , 1.0000054 , 1.0000054 ,
        1.0000054 , 1.0000054 , 1.0000054 , 1.0000054 , 1.0000054 ,
        1.0000054 , 1.0000054 , 1.0000054 , 1.0000054 , 0.9999991 ,
        0.99678932, 0.9912162 , 0.98475427, 0.97878585, 0.97635879,
        0.9754208 , 0.97466387, 0.97403457, 0.9735075 , 0.97306801,
        0.9727068 , 0.97241768, 0.9721964 , 0.97204013, 0.97194698,
        0.97191592, 0.97194662, 0.97203939, 0.9721953 , 0.97241618,
        0.9727049 , 0.97306567, 0.97350468, 0.97403119, 0.97465983,
        0.9754159 , 0.97635254, 0.97875657, 0.98471609, 0.99117968,
        0.99676151, 0.99999538, 1.0000054 , 1.0000054 , 1.0000054 ,
        1.0000054 , 1.0000054 , 1.0000054 , 1.0000054 , 1.0000054 ,
        1.0000054 , 1.0000054 , 1.0000054 , 1.0000054 , 1.0000054 ,
        1.0000054 ])}

Visualizing Results¶

We have several different methods (mostly wrappers to chromatic functions) to plot the modelled results. I'll demonstrate several of them below:

In [19]:

t.plot_lightcurves()

t.plot_lightcurves()

Check Residuals¶

In [20]:

t.plot_with_model_and_residuals()

t.plot_with_model_and_residuals()

No model attached to data. Running `add_model_to_rainbow` now. You can access this data later using [self].data_with_model

In [21]:

t.imshow_with_models()

t.imshow_with_models()

🌈🤖 'systematics_model' doesn't exist and will be skipped.

Plot the Transmission Spectrum¶

We can also plot the transmission spectrum (and we know what the true values are):

In [57]:

t.plot_transmission_spectrum(uncertainty=['hdi_16%','hdi_84%'])
plt.plot(r.wavelength, rbin5.metadata['transit_parameters']['rp_unbinned'], label="True Rp/R*")

plt.legend();

t.plot_transmission_spectrum(uncertainty=['hdi_16%','hdi_84%']) plt.plot(r.wavelength, rbin5.metadata['transit_parameters']['rp_unbinned'], label="True Rp/R*") plt.legend();

/Users/camu5866/opt/anaconda3/envs/chromaticfitting/lib/python3.9/site-packages/chromatic_fitting/models/transit.py:366: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  trans_table[f"{self.name}_radius_ratio_neg_error"] = (