Changing the priors

By default, kima uses a set of priors for the model parameters (see below). If you don't set new ones in the kima_setup.cpp file, the default ones will be used and the model will run.

But sometimes (most times?) you will want to consider different priors for some parameters.

To change specific priors, we just need to re-define some variables inside the RVmodel constructor.
The rule is: priors for the planet orbital parameters require the conditional object, those for other parameters don't.

For example, we can define two new priors

RVmodel::RVmodel():fix(true),npmax(1)
{
    // define new prior for the extra white noise parameter, s
    Jprior = make_prior<Uniform>(1., 20.); // m/s

    // get the conditional prior object
    auto conditional = planets.get_conditional_prior();
    // define new prior for the orbital period(s)
    conditional->Pprior = make_prior<Gaussian>(15, 0.1); // days
}

Note the small difference in syntax when defining the two priors.

Important:
The priors for the orbital parameters are the same for all planets. This is intentional!
If you're looking for a workaround, see this page.

The arguments between < > in the calls to make_prior (called template arguments) specify the prior distribution and the arguments inside parentheses set the parameters of each distribution. kima (and DNest4) currently implement the following distributions:

distribution	arguments	link
Uniform	lower, upper
LogUniform	lower, upper	scipy
ModifiedLogUniform	knee, upper
Rayleigh	scale
Triangular	lower, center, upper
Laplace	center, width	wiki
Kumaraswamy	alpha, beta	wiki
Gaussian	mean, stdev
Exponential	mean
Cauchy	location, scale	wiki
Fixed	value

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.

The kima-checkpriors script is provided with the pykima package to make histograms of the samples of each parameter. Simply call it with the column of the parameter you want to check (look inside sample.txt to see which one is which)

For example

kima-checkpriors 1

might print

Histogram of column 1: extra_sigma
number of samples: 10000
max value: 1999.744937
min value: 0.000128

and display the histogram of the prior samples

Default priors

// orbital periods and semi-amplitudes

// if hyperpriors=true
Pprior = make_shared<Laplace>(exp(log_muP), wP);
Kprior = make_shared<Exponential>(exp(log_muK));

// if hyperpriors=false
Pprior = make_shared<LogUniform>(1.0, 1e5);
Kprior = make_shared<ModifiedLogUniform>(1.0, 1e3);


// eccentricities
eprior = make_shared<Uniform>(0, 1);
    // another good option is
    // eprior = make_shared<Kumaraswamy>(0.867, 3.03);

// orbital phase
phiprior = make_shared<Uniform>(0, 2*M_PI);
// argument of periastron
wprior = make_shared<Uniform>(0, 2*M_PI);


// systemic velocity; m/s
Cprior = make_prior<Uniform>( data.get_RV_min(), data.get_RV_max() );
                          // so between the minimum and maximum RV

// additional white noise; m/s
Jprior = make_prior<ModifiedLogUniform>(
            min(1.0, 0.1*data.get_max_RV_span()), data.get_max_RV_span() );

// coefficients of the trend (if trend=true);
// units are m/s/day, m/s/day², m/s/day³
slope_prior = make_prior<Gaussian>( 0.0, pow(10, data.get_trend_magnitude(1)) );
quadr_prior = make_prior<Gaussian>( 0.0, pow(10, data.get_trend_magnitude(2)) );
cubic_prior = make_prior<Gaussian>( 0.0, pow(10, data.get_trend_magnitude(3)) );

// between-instrument offsets; m/s
offsets_prior = make_prior<Uniform>( -data.get_RV_span(), data.get_RV_span() );

// GP hyperparameters
log_eta1_prior = make_prior<Uniform>(-5, 5);
eta2_prior = make_prior<LogUniform>(1, 100);
eta3_prior = make_prior<Uniform>(10, 40);
log_eta4_prior = make_prior<Uniform>(-1, 1);

// degrees of freedom for Student-t likelihood
nu_prior = make_prior<LogUniform>(2, 1000);

// moving average parameters
sigmaMA_prior = make_prior<Uniform>(-1, 1);
tauMA_prior = make_prior<LogUniform>(1, 100);

// hyperparameters for orbital period and semi-amplitude hierachical priors
// only if hyperpriors=true
log_muP_prior = make_shared<TruncatedCauchy>(log(365), 1.0, log(365)-21, log(365)+21);
wP_prior = make_shared<Uniform>(0.1, 3.0);
log_muK_prior = make_shared<TruncatedCauchy>(0.0, 1.0, -21.0, 21.0);

// correlation coefficients with activity indicators
betaprior = make_prior<Gaussian>(0, 1);

In the example kima/examples/default_priors, kima will sample from (some of) these default priors.