Model Fitting

Before reading this section, it is recommended to go through the Getting Started page and check out the Orbital Decay of WASP-12 b and Long-Term Radial Velocity Trends examples.

To begin, a StarPlanet object is created that represents the star-planet system. The following code snippet, using WASP-12 b as an example, demonstrates the creation of a StarPlanet object:

from orbdot.star_planet import StarPlanet

# initialize the StarPlanet class
wasp12 = StarPlanet('examples/settings_files/WASP-12_settings.json')

See The Settings File for a description of the contents and structure of the input file.

Running Model Fits

To run the model fitting routines, one of the following methods is called on the StarPlanet object: run_joint_fit(), run_ttv_fit(), run_rv_fit(), or run_tdv_fit().

These methods have three key arguments in common:

Argument	Description
`free_params`	The free parameters for the model fit, formatted as a list of strings in any order.
`model`	The evolutionary model, must be `"constant"`, `"decay"`, or `"precession"`.
`file_suffix`	An optional string that is appended to the output file names.

The free parameters (free_params) can be given in any order, but an error will be raised if they are not part of the physical model. See Model Parameters for more information on the model parameters and their symbols.

The model argument specifies the evolutionary model to be fit, and must be either "constant" for the constant-period model, "decay" for orbital decay, or "precession" for apsidal precession. The default argument is "constant". See the Theory section for more information about the models.

The optional file_suffix argument enables separate fits of the same model to be differentiated. This functionality is nicely illustrated in Long-Term Radial Velocity Trends, which uses, for example, the arguments file_suffix="_circular" and file_suffix="_eccentric".

TTV Model Fits

The transit and eclipse timing fits are run by calling the run_ttv_fit() method on a StarPlanet object. The following code snippet shows an example call for each of the three evolutionary models:

wasp12.run_ttv_fit(['t0', 'P0'], model='constant')
wasp12.run_ttv_fit(['t0', 'P0', 'PdE'], model='decay')
wasp12.run_ttv_fit(['t0', 'P0', 'e0', 'w0', 'wdE'], model='precession')

TTV “clipping”

When fitting transit mid-times, there is an option to run a sigma-clipping routine to remove outliers in the transit mid-times, which may be useful for data with high variance [Hagey et al., 2022].

Passing sigma_clip=True to the run_ttv_fit() method runs the clip() method before the specified TTV model fit. Any subsequent model fits will use the cleaned data, so sigma_clip=True should only be specified once. For example,

wasp12.run_ttv_fit(['t0', 'P0'], model='constant', sigma_clip=True)
wasp12.run_ttv_fit(['t0', 'P0', 'PdE'], model='decay')
wasp12.run_ttv_fit(['t0', 'P0', 'e0', 'w0', 'wdE'], model='precession')

The clip() method operates by determining the best-fit transit timing model, subtracting it from the data, and then removing any data point with a nominal value that falls outside of a 3-\(\sigma\) range from the mean of the residuals. This process is repeated until no points fall outside the residuals, or until a maximum number of iterations has been reached.

RV Model Fits

The radial velocity model fits are run by calling the run_rv_fit() method on a StarPlanet object. The following code snippet shows an example call for each of the three evolutionary models:

wasp12.run_rv_fit(['t0', 'P0', 'K', 'v0', 'jit'], model='constant')
wasp12.run_rv_fit(['t0', 'P0', 'PdE', 'K', 'v0', 'jit'], model='decay')
wasp12.run_rv_fit(['t0', 'P0', 'e0', 'w0', 'wdE', 'K', 'v0', 'jit'], model='precession')

TDV Model Fits

Attention

The transit duration features of OrbDot have not been thoroughly tested and validated at this time. The methods are available to use, but the results should be treated with caution until this notice is removed.

The transit duration model fits are run by calling the run_tdv_fit() method on a StarPlanet object. The following code snippet shows an example call for each of the three evolutionary models:

wasp12.run_tdv_fit(['P0', 'ecosw', 'esinw', 'i0'], model='constant')
wasp12.run_tdv_fit(['P0', 'i0', 'PdE'], model='decay')
wasp12.run_tdv_fit(['P0', 'e0', 'w0', 'i0', 'wdE'], model='precession')

Joint Fits

Running a joint model fit is similar, but in this case the data types must be specified with at least two of the following arguments:

TTV=True: includes the transit and/or eclipse mid-times.

RV=True: includes the radial velocities.

TDV=True: includes the transit durations.

The following code snippet runs a joint fit of the mid-times and radial velocities for each of the three evolutionary models:

wasp12.run_joint_fit(['t0', 'P0', 'K', 'v0', 'jit'], model='constant', RV=True, TTV=True)
wasp12.run_joint_fit(['t0', 'P0', 'PdE', 'K', 'v0', 'jit'], model='decay', RV=True, TTV=True)
wasp12.run_joint_fit(['t0', 'P0', 'e0', 'w0', 'wdE', 'K', 'v0', 'jit'], model='precession', RV=True, TTV=True)

Fixed Parameter Values

The “fixed” values are assigned to any parameter that is not allowed to vary in a model fit. They are taken from the star-planet system info file, but may be updated at any time by calling the update_default() method. For example,

wasp12.update_default('P0', 3.14)

This is particularly useful for updating the fixed values between model fits. For example, the following code snippet runs a constant-period timing model fit, updates the fixed parameter values with the best-fit results, and then runs a radial velocity model fit:

# run the constant-period transit/eclipse timing model fit
ttv_fit = wasp12.run_ttv_fit(['t0', 'P0'], model='constant')

# update the default values for 'P0' and 't0'
wasp12.update_default('P0', ttv_fit['params']['P0'][0])
wasp12.update_default('t0', ttv_fit['params']['t0'][0])

# run the radial velocity model fit with 'P0' and 't0' fixed
wasp12.run_rv_fit(['K', 'v0', 'jit'], model='constant')

Priors

OrbDot currently supports three different prior distributions, the bounds of which are defined in the "prior" dictionary of the settings file.

The keys of "prior" are identical to the parameter symbols that are defined in the Model Parameters section. Every value is a list of three elements, the first being the type of prior ("uniform", "gaussian", or "log"), and the subsequent elements defining the distribution. This structure is illustrated in the following table:

Prior Type	Required Format	Example
Gaussian	`["gaussian", mean, std]`	`["gaussian", 2456305.5, 0.1]`
Uniform	`["uniform", min, max]`	`["uniform", -100, 100]`
Log-Uniform	`["log", min, max]`	`["log", -2, 1]`

There are default priors defined in the defaults/default_fit_settings.json file, but the user should, in general, specify them explicitly in the settings file. For example,

...
     "prior": {
        "t0": ["gaussian", 2456305.4555, 0.01],
        "P0": ["gaussian", 1.09142, 0.0001],
        "PdE": ["uniform", -1e-7, 0],
      }
}

Like the fixed values, the priors may be updated at any time by calling the update_prior() method. This is particularly useful for updating the priors in-between model fits. For example, the following code snippet runs a constant-period timing model fit, updates the priors with the best-fit results, and then runs a radial velocity model fit:

# run the constant-period transit/eclipse timing model fit
ttv_fit = wasp12.run_ttv_fit(['t0', 'P0'], model='constant')

# extract the best-fit results, structured as [value, upper_unc, lower_unc]
t0_best = ttv_fit['params']['t0']
P0_best = ttv_fit['params']['P0']

# update the priors for 'P0' and 't0'
wasp12.update_prior('P0', ['gaussian', P0_best[0], P0_best[1]])
wasp12.update_prior('t0', ['gaussian', t0_best[0], t0_best[1]])

# run the radial velocity model fit with 'P0' and 't0' as free parameters
wasp12.run_rv_fit(['t0', 'P0', 'K', 'v0', 'jit'], model='constant')

Output Files

At the end of every model fit, the following files are saved:

"*_summary.txt": a quick visual summary of the results

"*_results.json": the entire model fitting results dictionary

"*_corner.png": a corner plot

"*_weighted_samples.txt": the weighted posterior samples

"*_random_samples.json": a random set of 300 posterior samples

The `"*_summary.txt"` File

This file provides a concise overview of the results of the model fit in an easy-to-read text format. For example, the following output is from a fit of the orbital decay model to the transit and eclipse mid-times of WASP-12 b (see the WASP-12 b example for more):

Stats
-----
Sampler: nestle
Free parameters: ['t0' 'P0' 'PdE']
log(Z) = -104.47 ± 0.14
Run time (s): 7.04
Num live points: 1000
Evidence tolerance: 0.01
Eff. samples per second: 663

Results
-------
t0 = 2456305.4558077552 + 3.379490226507187e-05 - 3.208918496966362e-05
P0 = 1.0914201076440608 + 4.156631039364811e-08 - 4.3833844109997244e-08
PdE = -1.00348670058712e-09 + 6.98096735732343e-11 - 6.878773061871802e-11
dPdt (ms/yr) = -29.015070989305705 + 2.0184947476459363 - 1.9889460278124174

Fixed Parameters
----------------
e0 = 0.0
w0 = 0.0

The `"*_results.json"` File

This file stores a comprehensive summary of the model fit results and settings in .json format. It ensures that critical information about the model fit is not lost, but it is not designed for easy reading. Rather, the "*_summary.txt" file serves to quickly convey the results and should typically be examined first.

The following table lists the keys of the *_results.json file:

Key	Data Type	Description
`"stats"`	`dict`	A dictionary containing various model fit statistics and settings.
`"params"`	`dict`	A dictionary containing the best-fit parameters and their 68% credible intervals.
`"prior"`	`dict`	The dictionary of prior distributions from the settings file.
`"model"`	`str`	The model that was fit (e.g. `"ttv_constant"`, `"joint_precession"`, etc.).
`"suffix"`	`str`	The file suffix that was given to the model fit.
`"results_filename"`	`str`	The path to this file (recorded for the plotting functions).
`"samples_filename"`	`str`	The path to the `"*_random_samples.txt"` file (recorded for the plotting functions).

The "stats" dictionary records various model fit statistics and settings with the following keys:

`"logZ"`	`float`	The Bayesian evidence.
`"logZ_err"`	`float`	The Bayesian evidence uncertainty.
`"run_time"`	`float`	The run time of the model fit in seconds.
`"evidence_tolerance"`	`float`	The evidence tolerance given to the model fit.
`"n_live_points"`	`float`	The number of live points given to the model fit.
`"n_dims"`	`float`	The number of free parameters.
`"n_samples"`	`float`	The number of weighted posterior samples.
`"eff_samples_per_s"`	`float`	The effective samples per second.

The "params" dictionary contains key-value pairs that store the best-fit parameter values and their 68% credible intervals. The keys match the parameter symbols (see Model Parameters), and each value is a list of three elements: [best-fit value, upper uncertainty, lower uncertainty].

The following code snippet demonstrates how to access the best-fit parameters after a model fit:

# run the constant-period timing model fit
ttv_fit = wasp12.run_ttv_fit(['t0', 'P0'], model='constant')

# extract the best-fit parameter values and their uncertainties
t0_best, t0_upper_err, t0_lower_err = ttv_fit['params']['t0']
P0_best, P0_upper_err, P0_lower_err = ttv_fit['params']['P0']

If the free parameters include "ecosw" and "esinw" or "sq_ecosw" and "sq_esinw", the derived eccentricity "e0" and argument of pericenter "w0" can be accessed the same way.

The entire set of OrbDot parameters are included in the "params" dictionary for completeness, even if they are not part of the physical model, to ensure that no information is lost or overlooked. If a parameter was not allowed to vary in the model fit, its fixed value is given.

Interpreting the Results

The Analyzer class combines model fit results, star-planet system characteristics, and the data to compute and summarize analyses of various physical models, such as equilibrium tides, apsidal precession, systemic proper motion, and companion objects.

The initialization of an Analyzer class requires a StarPlanet object and the results of a model fit. The latter may be passed directly after a model fit, for example:

# run the orbital decay TTV model fit
decay_fit = wasp12.run_ttv_fit(['t0', 'P0', 'PdE'], model='decay')

# initialize the Analyzer class
analyzer = Analyzer(wasp12, decay_fit)

or loaded from a preexisting file:

import json

# load the orbital decay fit results
with open('results/WASP-12/ttv_fits/ttv_decay_results.json') as jf:
    decay_fit = json.load(jf)

# initialize the Analyzer class
analyzer = Analyzer(wasp12, decay_fit)

As soon as an Analyzer object is created, a text file is generated for recording the output of any methods that are called. For example, the code snippet above generates the file: results/WASP-12/analysis/ttv_decay_analysis.txt.

`Analyzer` Methods

The following sections summarize the main Analyzer methods, the output of which are appended to the *_analysis.txt file described above.

1. Model Comparison

The model_comparison() method compares the Bayesian evidence for the model fit with that of a different model. To compare the two models, the Bayes factor is calculated as:

\[\log{B_{12}} = \log{\mathrm{Z}}_{1} - \log{\mathrm{Z}}_{2}\]

where \(\log{\mathrm{Z}}\) is the Bayesian evidence, which is defined such that a lower magnitude signifies a superior fit to the observed data. The Bayes factor is then compared to the thresholds established by Kass and Raftery [1995], tabulated below:

Condition	Evidence for Model 1 (Model 1)
\(B_{12} \leq 1\)	Model 1 is not supported over Model 2
\(1 < B_{12} \leq 3\)	Evidence for Model 1 barely worth mentioning
\(3 < B_{12} \leq 20\)	Positive evidence for Model 1
\(20 < B_{12} \leq 150\)	Strong evidence for Model 1
\(150 < B_{12}\)	Very strong evidence for Model 1

The following code snippet calls this method after running two different TTV model fits:

# run the apsidal precession TTV model fit
decay_fit = wasp12.run_ttv_fit(['t0', 'P0', 'PdE'], model='decay')
precession_fit = wasp12.run_ttv_fit(['t0', 'P0', 'e0', 'w0', 'wdE'], model='precession')

# initialize the Analyzer class
analyzer = Analyzer(wasp12, decay_fit)

# compare the orbital decay and apsidal precession models
analyzer.model_comparison(precession_fit)

2. Orbital Decay Model Fit

The orbital_decay_fit() method produces a summary of various values derived from interpreting the results of an orbital decay model fit in the context of equilibrium tidal theory.

# run an analysis of the orbital decay model fit results
analyzer.orbital_decay_fit()