# Incorporating linear models¶

In previous `juliet`

tutorials for transits (Lightcurve fitting with juliet) and radial-velocities (Fitting radial-velocities), we have so far
assumed that the only deterministic signals under consideration in the models \(\mathcal{M}_i(t)\) for instrument
\(i\) are composed of underlying physical processes. For transits, we assume the function is a transit model
distorted both by a normalization constant and a dilution factor, whereas for the radial-velocities we assume
this is an addition between a Keplerian signal, a systemic radial-velocity and a long-term trend. Typically, however,
these are not the only components that make up a model. For transits, systematics in the data (e.g., airmass trends,
meridian flips, etc.) can distort the signals further — for radial-velocities some linear models might help out
constrain activity signals.

Within `juliet`

one can model, in addition to the deterministic signal for transits and radial-velocities,
\(\mathcal{M}_i(t)\), a linear model such that the full data-generating process can be written as

\(\mathcal{M}_i(t) + \textrm{LM}_i(t) + \epsilon_i(t)\),

where the terms \(\mathcal{M}_i(t)\) is the transit or radial-velocity model, \(\epsilon_i(t)\) is the noise model (for details on those, see previous tutorials on transits and radial-velocities), and where \(\textrm{LM}_i(t)\) is a linear model given by:

\(\textrm{LM}_i(t) = \sum_{n=0}^{p_i}x_{n,i}(t) \theta_{n,i}^{\textrm{LM}}\).

Here, the \(x_{n,i}(t)\) are the \(p_i+1\) linear regressors at time \(t\) for instrument \(i\), and the \(\theta_{n,i}^{\textrm{LM}}\) are the coefficients of those regressors (e.g., \(x_{n,i}(t) = t^n\) would model a polynomial trend for instrument \(i\)).

## Linear models in transit fits¶

Adding linear terms to a model within `juliet`

is very simple, and can be done in two ways. One way is to simply pack the lightcurve
and regressors in a text file of the form:

```
2458459.7999999998 1.0126748331 0.0030000000 CHAT 1.2107127967
2458459.8013377925 1.0127453892 0.0030000000 CHAT 1.2107915485
2458459.8026755853 1.0158682599 0.0030000000 CHAT 1.2108919775
2458459.8040133780 1.0117892069 0.0030000000 CHAT 1.2110140837
2458459.8053511707 1.0125201749 0.0030000000 CHAT 1.2111578671
2458459.8066889634 1.0133562197 0.0030000000 CHAT 1.2113233277
.
.
.
```

where, the first column saves the times, second the relative fluxes, third errors on these relative fluxes,
fourth the instrument names and the \(p_i+1\) subsequent columns store the \(p_i+1\) linear regressors
to be fitted to the data (in the above example, 1). Once this file is created, the filename can be simply given to the
`juliet.load`

call with the `lcfilename`

parameter — this will store the times, lightcurves and linear regressors
in a given dataset. The second way is to simply pass all the linear regressors using the `linear_regressors_lc`

variable
of the `juliet.load`

call — the input should be a dictionary, where each key is a different instrument and contains
an array of dimensions \((N_i, p_i+1)\), where \(N_i\) is the number of datapoints for instrument \(i\). In this
tutorial, we will use the former way of fitting linear models.

In this tutorial we will use the dataset uploaded [here] —
this dataset has one linear regressor. For each linear regressor, we must define the prior for the coefficient \(\theta_{n,i}\);
these are expected to be of the form `thetaN_i`

, where `N`

is the numbering of the linear regressor (as given in the
file or dictionary) and `i`

is the instrument name. In our case, we have data from the *CHAT* telescope — let’s fit it assuming a
linear model:

```
import juliet
import numpy as np
priors = {}
# Name of the parameters to be fit:
params = ['P_p1','t0_p1','r1_p1','r2_p1','q1_CHAT','q2_CHAT','ecc_p1','omega_p1',\
'rho', 'mdilution_CHAT', 'mflux_CHAT', 'sigma_w_CHAT', 'theta0_CHAT']
# Distributions:
dists = ['fixed','normal','uniform','uniform','uniform','uniform','fixed','fixed',\
'loguniform', 'fixed', 'normal', 'loguniform', 'uniform']
# Hyperparameters
hyperps = [3.1, [2458460,0.1], [0.,1], [0.,1.], [0., 1.], [0., 1.], 0.0, 90.,\
[100., 10000.], 1.0, [0.,0.1], [0.1, 1000.],[-100,100]]
# Populate the priors dictionary:
for param, dist, hyperp in zip(params, dists, hyperps):
priors[param] = {}
priors[param]['distribution'], priors[param]['hyperparameters'] = dist, hyperp
# Load dataset:
dataset = juliet.load(priors=priors, lcfilename = 'lc_lm.dat', out_folder = 'lm_transit_fit')
results = dataset.fit(n_live_points = 300)
```

Now let’s plot it:

```
t0 = np.median(results.posteriors['posterior_samples']['t0_p1'])
# Plot. First extract model:
transit_model, transit_up68, transit_low68, components = results.lc.evaluate('CHAT', return_err=True, \
return_components = True, \
all_samples = True)
import matplotlib.pyplot as plt
plt.errorbar(dataset.times_lc['CHAT']-t0, dataset.data_lc['CHAT'], \
yerr = dataset.errors_lc['CHAT'], fmt = '.' , alpha = 0.1)
# Out-of-transit flux:
oot_flux = np.median(1./(1. + results.posteriors['posterior_samples']['mflux_CHAT']))
# Plot non-transit model::
plt.plot(dataset.times_lc['CHAT']-t0, oot_flux + components['lm'], color='grey', lw = 3, label = 'Linear model + oot flux')
plt.plot(dataset.times_lc['CHAT']-t0, transit_model, color='black', label = 'Full model')
plt.fill_between(dataset.times_lc['CHAT']-t0,transit_up68,transit_low68,\
color='cornflowerblue',alpha=0.5,zorder=5)
plt.xlabel('Time from mid-transit (days)')
plt.ylabel('Relative flux')
plt.legend()
```