Multiple instruments¶

It is quite typical to have RV observations from more than one instrument. These data should be analysed together for the detection of planets.

Note

• Considering RVs from multiple instruments means adding new offset parameters to the model. By default, all these parameters share the same prior, but we'll see below how to change this.
• The same applies to the additional white noise (jitter). One jitter parameter is added per instrument, all sharing the same prior.

New

Newer versions of kima now allow for a stellar jitter parameter which is shared between all the instruments. By default, this parameter is fixed to zero to recover the old behaviour.

First of all, let's import the package

In [1]:

import kima

import kima

Detecting the planet around HD106252¶

First announced by Fischer et al. (2002), the giant planet orbiting HD106252 has been detected with a number of different instruments (Perrier et al. 2003, Butler et al. 2006, Wittenmyer et al.2009). In this example, we look at data from ELODIE, HET, HJS, and Lick.

Each of the datasets can be imported like this

In [2]:

from kima.examples.multi_instruments import HD106252_ELODIE

from kima.examples.multi_instruments import HD106252_ELODIE

and a simple visualisation of the RVs is provided by the plot method

In [3]:

HD106252_ELODIE.plot();

HD106252_ELODIE.plot();

The combined data is also readily available

In [4]:

from kima.examples.multi_instruments import HD106252_combined

HD106252_combined.plot();

from kima.examples.multi_instruments import HD106252_combined HD106252_combined.plot();

Note the clear offset between ELODIE data (in blue above) and the other datasets. The average RV has been subtracted from the HET, HJS, and Lick data, but not from ELODIE. This slightly convoluted situation is ideal to demonstrate how kima deals with multiple instruments, especially with the default priors.

We can now use one of the available models to analyse this dataset. We can either create the model

In [5]:

from kima import RVmodel

model = RVmodel(fix=False, npmax=1, data=HD106252_combined)

from kima import RVmodel model = RVmodel(fix=False, npmax=1, data=HD106252_combined)

or import the multi_instruments function from the examples, which does the same thing (but can also run the model directly)

In [6]:

from kima.examples import multi_instruments

model = multi_instruments(run=False)

from kima.examples import multi_instruments model = multi_instruments(run=False)

Notice how we set the number of planets to be free, implicitly assigning to it a uniform prior between 0 and 1 (by setting npmax=1).

This is all that is necessary to setup the data and model, and we are now ready to run the analysis. Notice that we didn't explicitly assign any priors (except for the number of planets). This means kima will use all the default priors.

Let's run the model for a few thousand steps and load the results

In [7]:

model, res = multi_instruments(run=True, load=True, steps=30_000)

model, res = multi_instruments(run=True, load=True, steps=30_000)

log(Z) = -498.67164265769134
Information = 67.6310091586351 nats.
Effective sample size = 2984.7236142899683

100%|██████████| 2984/2984 [00:00<00:00, 8048.22it/s]

Did we find HD106252 b?

In [8]:

res.plot_posterior_np();

res.plot_posterior_np();

Np probability ratios:  []

Yes! In fact, even though the $N_p$ parameter was free, all the samples have $N_p=1$ because the posterior probability is so much higher than for $N_p=0$. So we're pretty sure about this detection. How does the fit look?

In [9]:

res.plot_random_samples();

res.plot_random_samples();

The plot_random_samples function just shows the Keplerian curves from a few random posterior samples, together with the data. Note that the RV offsets for each instrument are not subtracted from the data. This is on purpose: since we're showing several posterior samples (50 by default), there are actually several values for the offsets. Zooming in on one particular dataset helps (note the several blue curves):

In [10]:

fig = res.plot_random_samples();
fig.axes[0].set_ylim(9600, 10200);

fig = res.plot_random_samples(); fig.axes[0].set_ylim(9600, 10200);

Looking at the posterior for the orbital period, semi-amplitude, and eccentricity of the planet clearly shows that the default priors are a bit too wide

In [11]:

res.plot_posterior_PKE(show_prior=True);

res.plot_posterior_PKE(show_prior=True);

The situation is similar for the systemic velocity and for the RV offsets, as evidenced in the plots below, where we might not even see the posterior!

In [12]:

res.hist_vsys(show_prior=True)

res.hist_vsys(show_prior=True)

Out[12]:

[<Figure size 640x480 with 1 Axes>, <Figure size 1100x500 with 3 Axes>]