Models
Types of models
kima currently implements dedicated models for different analyses of a given dataset. The models share a common organization, but each has its own parameters (and thus priors) and settings.

RVmodel
models the RVs with a sumofKeplerians 
GPmodel
models the RVs with a sumofKeplerians plus a correlated noise component given by a Gaussian process 
RVFWHMmodel
models the RVs together with the FWHM as an activity indicator, including a Gaussian process for the activity signal 
BINARIESmodel
this includes tidal and relativistic effects as well as apsidal precession of the binary’s orbit (applicable for circumbinary planets), and can also fit doublelined binary data.
The type of model to be used is set by defining a constructor in the
kima_setup
file
or
Info
A constructor in C++ is a special function that instantiates an object, in
this case of type RVmodel
or RVFWHMmodel
. If this sounds weird, think of
it as the __init__
function in a Python class. If this also sounds weird,
just think of it as a simple function where we implement the model.
Since each model has its own parameters and settings, they are described separately below.
This is the most basic model, but already allows for an indepth analysis of an RV dataset. In terms of general settings, we can define
RVmodel::RVmodel() : fix(true), npmax(1) // (4)
{
trend = false / true;
degree = 0 / 1 / 2 / 3; // (1)
known_object = false / true; // (2)
n_known_object = 0 / ... ;
studentt = false / true; // (3)
}

Sets up a longterm trend in the RVs, of a given degree. This is a linear trend for
degree=1
, a quadratic fordegree=2
, and a cubic fordegree=3
. Higher degree trends are not implemented. 
In known object (KO!) mode, kima considers some Keplerians as part of a "background" model (see below). This can be useful when modelling transiting planets or simply planets with betterknown orbital parameters.

Whether to use a Studentt distribution as the likelihood, instead of a Gaussian. If this is set to
true
, the degrees of freedom of the Studentt distribution can be estimated from the data. This is sometimes useful when the RV data is suspected to contain (a few) strong outliers. See the example for more information. 
See below for more information
model = RVmodel(fix=True, npmax=1) # (4)
model.trend = False / True
model.degree = 0 / 1 / 2 / 3 # (1)
model.known_object = False / True # (2)
model.n_known_object = 0 / ...
model.studentt = False / True # (3)

sets up a longterm trend in the RVs, of a given degree. This is a linear trend for
degree=1
, a quadratic fordegree=2
, and a cubic fordegree=3
. Higher degree trends are not implemented. 
In known object (KO!) mode, kima considers some Keplerians as part of a "background" model (see below). This can be useful when modelling transiting planets or simply planets with betterknown orbital parameters.

Whether to use a Studentt distribution as the likelihood, instead of a Gaussian. If this is set to
True
, the degrees of freedom of the Studentt distribution can be estimated from the data. This is sometimes useful when the RV data is suspected to contain (a few) strong outliers. See the example for more information. 
See below for more information
This model considers a Gaussian process as a surrogate model for the stellar activity signal present in the RVs. We can specify different kernels each with its own hyperparameters and infer their posterior as well. By default, the standard quasiperiodic kernel is used, with four hyperparameters:
Most of the same settings as for the RVmodel are available, except for the
studentt
which doesn't apply in this case.
GPmodel::GPmodel() : fix(true), npmax(1) // (4)
{
trend = false / true;
degree = 0 / 1 / 2 / 3; // (1)
known_object = false / true; // (2)
n_known_object = 0 / ... ;
}

Sets up a longterm trend in the RVs, of a given degree. This is a linear trend for
degree=1
, a quadratic fordegree=2
, and a cubic fordegree=3
. Higher degree trends are not implemented. 
In known object (KO!) mode, kima considers some Keplerians as part of a "background" model (see below). This can be useful when modelling transiting planets or simply planets with betterknown orbital parameters.

See below for more information
model = GPmodel(fix=True, npmax=1) # (4)
model.trend = False / True
model.degree = 0 / 1 / 2 / 3 # (1)
model.known_object = False / True # (2)
model.n_known_object = 0 / ...

sets up a longterm trend in the RVs, of a given degree. This is a linear trend for
degree=1
, a quadratic fordegree=2
, and a cubic fordegree=3
. Higher degree trends are not implemented. 
In known object (KO!) mode, kima considers some Keplerians as part of a "background" model (see below). This can be useful when modelling transiting planets or simply planets with betterknown orbital parameters.

See below for more information
The Np planets
kima can model Keplerian functions, the RV signals induced by planets. But the distinguishing feature is that the code can actually model a number \(N_p\) of Keplerians, and this number does not need to be fixed a priori. In other words, kima estimates the joint posterior distribution for \(N_p\) and the orbital parameters of the planets \(\theta\), given some data \(\mathcal{D}\):
So \(N_p\) is actually a free parameter in the model (if we want). It is somewhat of a special parameter, in that it only takes discrete values, but it's otherwise similar to other parameters in the model. In particular, it needs a prior distribution as well.
This part of the model in kima can be set with a special pair of parameters for the model constructors
which define whether of not \(N_p\) is fixed and its value or prior distribution
By default, each of the \(N_p\) planets has 5 orbital parameters
corresponding to the orbital period, semiamplitude, eccentricity, mean anomaly at the epoch, and argument of periastron.
Some models include additional perplanet parameters: the BINARIESmodel considers a linear precession parameter \(\dot\omega\) and the BDmodel infers a mixture probability \(\lambda\) for each planet.
units
\(P\) is in [days], \(K\) is in [m/s], \(e\) is unitless, and \(M_0\) and \(ω\) are in radians.
By default, the epoch is set to the first observed time, but it can be changed
inside the constructor by setting the variable M0_epoch
.
The "background" model
Besides the \(N_p\) planets, kima models the RV observations as follows
 by default
 if using a Studentt likelihood
 if using a Gaussian process
where \(j\) is a jitter parameter that represents additional white noise variations not accounted for by the individual uncertainties.
When considering a longterm trend (by setting trend=true
), instead of a
constant systemic velocity, we have
where \(t_m\) is the mean of the times and slope, quadr, and cubic are free parameters. See the example for more information.
When using data from multiple instruments, kima adds RV offsets between pairs of instruments as well as individual jitter parameters per instrument.
Model parameters
How many parameters?
The number of free parameters can change substantially with each model, with the different settings, and with the specific parameter priors.
Prior distributions
As in any Bayesian analysis, kima needs a set of priors for the model
parameters. If we don't explicitly set the priors in the kima_setup
file, a
default prior will be used and the model will run. But sometimes (most times?)
we will want to set custom priors for some parameters.
To change specific priors, we just need to redefine some variables inside the model constructor.
Warning
Note below that, on the C++ side, the priors for the planet orbital
parameters require the conditional
object, while the priors for the other
parameters do not. In Python, this is automatically taken care of.
For example, let's redefine the priors for the jitter and for the orbital periods of the \(N_p\) planets
RVmodel::RVmodel()
{
// define new prior for the extra white noise parameter (jitter), j
Jprior = make_prior<Uniform>(1.0, 20.0);
// define new prior for the orbital period(s)
// first, get the conditional prior object
auto conditional = planets.get_conditional_prior();
conditional>Pprior = make_prior<Gaussian>(15.0, 0.1); // in days
}
from scipy.stats import uniform, norm
model = RVmodel()
# define new prior for the extra white noise parameter (jitter), j
model.priors.J = uniform(1, 20) # (1)
# define new prior for the orbital period(s)
model.priors.P = norm(15, 0.1)
 Note that this is a uniform distribution from 1 m/s to 20 m/s,
unlike the
scipy
distribution which would have support [1, 21] !
While it might seem counterintuitive at first, in kima the priors for the orbital parameters are the same for all \(N_p\) planets. This is intentional! But if this is not what you want, consider using the known object feature.
In the example above, we specify the prior distribution and its parameters. kima (and DNest4) currently implement the following distributions:
in C++  in Python  arguments  link 

Uniform 
scipy.stats.uniform 
lower, upper  wiki, scipy^{!} 
LogUniform 
scipy.stats.loguniform 
lower, upper  scipy 
ModifiedLogUniform 
knee, upper  
Rayleigh 
scipy.stats.rayleigh 
scale  
Triangular 
lower, center, upper  wiki, scipy  
Laplace 
center, width  wiki  
Kumaraswamy 
pykima.priors.kumaraswamy 
alpha, beta  wiki 
Gaussian 
scipy.stats.norm 
mean, stdev  wiki, scipy 
Exponential 
mean  wiki  
Cauchy 
scipy.stats.cauchy 
location, scale  wiki, scipy 
Fixed 
value  
TruncatedGaussian 
mean, stdev, lower, upper  
TruncatedCauchy 
location, scale, lower, upper 
If you need a distribution that is not yet implemented, consider opening an issue.
Checking if the priors are correct
For any kima model, it's very easy to check if the priors are set up
correctly. In the OPTIONS
file, set the maximum number of levels to 1, and run
the model as usual. With this setting, DNest4 will sample directly from the
priors and ignore the likelihood.
The kimacheckpriors
script is provided with the pykima
package to make
histograms of the samples of each parameter. Simply call it with the column or
the name of the parameter you want to check (look inside sample.txt
to see
which one is which)
For example
might print
Histogram of column 1: jitter
number of samples: 10000
max value: 1999.744937
min value: 0.000128
and display the histogram of the prior samples:
Default priors
Below is the list of default priors which are used if not explicitly redefined in the model constructor.
expand

orbital periods and semiamplitudes
units: days and m/s, respectively 
orbital eccentricities
another good option could be
which approximates the Beta distribution proposed by Kipping (2013). 
orbital phase and argument of periastron
units: radians```c++ phiprior = make_shared
(0, 2*PI); wprior = make_shared
(0, 2*PI); ``` 
systemic velocity
units: m/s 
additional white noise, or jitter
units: m/s 
ccoefficients of the trend (if trend=true)
units: m/s/day, m/s/day², m/s/day³ 
betweeninstrument offsets units: m/s
{% katexmm %}

GP hyperparameters
units: \(\eta_1\) in m/s, \(\eta_2\) and \(\eta_3\) in days
{% endkatexmm %}

degrees of freedom for Studentt likelihood

moving average parameters

hyperparameters for orbital period and semiamplitude hierachical priors
only ifhyperpriors=true

correlation coefficients with activity indicators