Running The Joker with custom prior samples

By default, The Joker will generate prior samples internally and use these to generate posterior samples. It is sometimes useful to instead create a custom cache of prior samples to use if you, for example, believe that your data or systems of interest have a different eccentricity distribution than that assumed in The Joker. Here, we’ll work through an example of how to generate your own prior samples file and use this with The Joker. In particular, we will assume that the system has no eccentricity (is circular), so we will generate prior samples with \(e = 0\).

[1]:
import os

from astropy.time import Time
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

from thejoker.data import RVData
from thejoker.sampler.io import save_prior_samples
from thejoker.sampler import JokerParams, TheJoker
from thejoker.plot import plot_rv_curves

For this example, I’ve already pre-generated some simulated data that we’re going to use – the true parameters of the main orbit are in the cell below:

[2]:
true_P = 61.1166 * u.day
true_e = 0.0324

The data are provided in-line in the next cell:

[3]:
t = Time([56001.71633991, 56030.37736291, 56033.3716085 , 56035.26497108,
          56037.53465393, 56131.88500321, 56145.72804619, 56215.01481418],
         format='mjd', scale='tcb')
rv = [19.911927, 12.004762, 4.320121, 0.21907635,
      -3.8032671, 36.162901, 28.636109, 8.2394115] * u.km/u.s
rv_err = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2] * u.m/u.s
data = RVData(t, rv, stddev=rv_err)
_ = data.plot()
../_images/examples_custom-prior-sampling_6_0.png

The easiest way to generate your own prior samples is to use the TheJoker.sample_prior() method, then modify the samples that are returned. Here we will use the default period prior (uniform in the logarithm of period between the specified bounds), but will modify the prior samples in eccentricity that are generated:

[4]:
params = JokerParams(P_min=2*u.day, P_max=1024*u.day)
joker = TheJoker(params)

Here we use the .sample_prior() method to generate a million prior samples under the default prior, then we replace the eccentricity, "e", values with zeros:

[5]:
prior_samples = joker.sample_prior(size=1_000_000)

default_e_samples = prior_samples['e']
prior_samples['e'] = np.zeros(len(prior_samples))

To use these samples with The Joker, we have to create a cache file and save the prior samples. We can do that with the save_prior_samples() function in thejoker:

[6]:
prior_cache_file = 'prior_cache.hdf5'

# always overwrite the file:
if os.path.exists(prior_cache_file):
    os.remove(prior_cache_file)

save_prior_samples(prior_cache_file,
                   prior_samples,
                   data.rv.unit);

We now pass this cache filename in to rejection_sample or iterative_rejection_sample:

[7]:
samples = joker.rejection_sample(data, prior_cache_file=prior_cache_file)
INFO: 1 good sample after rejection sampling [thejoker.sampler.sampler]
[8]:
_ = plot_rv_curves(samples, data=data)
../_images/examples_custom-prior-sampling_15_0.png

Notice that the samples (in this case, just 1) that are returned have eccentricity \(e=0\), as we desired:

[9]:
samples
[9]:
JokerSamples([('P', <Quantity [60.99474123] d>),
              ('M0', <Quantity [1.41576253] rad>),
              ('e', <Quantity [0.]>),
              ('omega', <Quantity [0.1012441] rad>),
              ('jitter', <Quantity [0.] km / s>),
              ('K', <Quantity [25.02450126] km / s>),
              ('v0', <Quantity [13.8589499] km / s>)])