Getting started

After installing kima, we'll go through one of the examples to start playing with the package.

Let's find 51 Peg b!
In the kima/examples/51Peg folder you will find a dataset of radial-velocity observations and the code to analyse it with kima.

51 Peg b was the first exoplanet discovered around a solar-type star. The original paper, by Mayor & Queloz (1995), used data from the ELODIE spectrograph. Here we will use another dataset, obtained with the Hamilton echelle spectrograph, situated in the Lick Observatory, California (see Butler et al. 2006).

The data is in the 51Peg.rv file. It has 256 observations separated in three columns: time (in modified Julian Days), measured radial velocity (in m/s) and its uncertainty (also in m/s).
This type of file is the typical input to kima (see this page for more help).

# 51Peg.rv
50002.665695	-52.9	 4.1
50002.684340	-45.8	 4.8
50002.800220	-60.8	 4.6
50002.815961	-53.3	 5.0
...

Our goal is to fit a Keplerian model to these radial-velocity (RV) observations.
To do this with kima, we prepare a kima_setup.cpp file (there's one already in the folder, we'll just go through it).

This C++ file starts with a general include statement

#include "kima.h"

Afterwards, we define general settings for the model

const bool GP = false;
const bool MA = false;
const bool hyperpriors = false;
const bool trend = false;
const int degree = 0;
const bool multi_instrument = false;
const bool known_object = false;
const int n_known_object = 0;
const bool studentt = false;

These settings tell kima that

• we will use a Keplerian model without considering correlated noise (instead of a Gaussian process model or a moving average model),
• we will not consider hyper-priors for the orbital parameters
• no RV linear trend
• since we only have RV data from one instrument (Lick) we also set the flag multi_instrument = false
• the known_object mode is set to false
• and the likelihood is a Gaussian (instead of a Student-t distribution)

Next comes the RVmodel constructor, where we set a few more options for the model, and the priors.

RVmodel::RVmodel():fix(false),npmax(1)
{
    // systemic velocity; m/s
    Cprior = make_prior<Uniform>(-10., 10.);
    // additional white noise; m/s
    Jprior = make_prior<ModifiedLogUniform>(1.0, 1000.);

    auto conditional = planets.get_conditional_prior();
    // orbital period; days
    conditional->Pprior = make_prior<LogUniform>(0.2, 2000.); // days
    // semi-amplitude; m/s
    conditional->Kprior = make_prior<ModifiedLogUniform>(1., 1000.);
    // eccentricity
    conditional->eprior = make_prior<Uniform>(0., 1.);
    // orbital angles, w and phi; radians
    conditional->phiprior = make_prior<Uniform>(0., 2*PI);
    conditional->wprior = make_prior<Uniform>(0., 2*PI);
}

The first line sets the maximum number of planets in the model to 1, and tells kima not to fix it during the analysis (by setting fix to false). That is, the number of planets is a free parameter, having a uniform prior between 0 and npmax (1 in this case).

We then define priors for all the other parameters. In this example we chose very wide priors for demonstration purposes. More details about how to setup the priors can be found in this page. If we don't redefine a parameter's prior in the RVModel() function, it will be assigned a default one. Those defaults are described here.

We're almost ready. All that's left is to read the data file and run the model, which is done inside a main function

int main(int argc, char** argv)
{
    // define the filename with the RV data
    datafile = "51Peg.rv";

    // and load it
    // the second argument is the units of the RVs
    // the third (optional) argument, tells kima 
    // not to skip any line in the header of the file
    load(datafile, "ms", 0);
    
    // set the sampler and run it!
    Sampler<RVmodel> sampler = setup<RVmodel>(argc, argv);
    sampler.run(50); // print to terminal every 50 saves
    return 0;
}

And that's it, we can now compile and run the model

kima-run

kima will start (using 4 threads, by default) and print some information to the terminal. The full example, with the default options, will take around 5 min to finish.

In the terminal, kima will say that it read the RVs from the file, is creating levels at given likelihood values and saving particles. The default option is to stop after saving 5000 particles.

But we don't have to wait for it to finish!

If we open another terminal in the same folder (kima/examples/51Peg), we can use the kima-showresults script included in the pykima package to look at the results.

You can run

kima-showresults -h

to see all the available options.

Calling the script with all

kima-showresults all

will display a number of plots with results from the run.
Note: there's no problem if you run kima-showresults while kima is still running!

A number of figures will show up:

• the posterior distribution for the extra white noise,
  which should peak around 3 m/s
• the posterior distribution for the systemic velocity,
  which should peak around -1.8 m/s
• a plot of posterior samples on top of the RV data,
  you will need to zoom in to see the Keplerian orbit
• the phased RV curve
• the posteriors for semi-amplitude, eccentricity, and orbital period
• the posterior for the orbital period,
  which should show a very strong peak around 4.22 days
• and the posterior for the number of planets

This last plot will likely show only samples with Np=1, even though the number of planets is a free parameter. In this case, the evidence for one planet is so high that the ratio of probabilities
p(Np=1) / p(Np=0) is (much) higher than the number of posterior samples we got.

So, kima tells us that a planet at a period of 4.22 days and an amplitude of 55 m/s is the best model given the RV data. If it was 1995, we'd be going about changing the history of astronomy!