For models with multiple instruments, this function "bursts" the
computed RV into n_instruments individual arrays. This is mostly
useful for plotting the RV model together with the original data.
Evaluate the deterministic part of the model at one posterior sample.
If t is None, use the observed times. Instrument offsets are only
added if t is None, but the systemic velocity is always added.
To evaluate at all posterior samples, consider using
np.apply_along_axis(self.eval_model, 1, self.posterior_sample)
Note: this function does not evaluate the GP component of the model.
Arguments
sample : ndarray
One posterior sample, with shape (npar,)
t : ndarray
Times at which to evaluate the model, or None to use observed times
include_planets : bool, default True
Whether to include the contribution from the planets
include_known_object : bool, default True
Whether to include the contribution from the known object planets
include_trend : bool, default True
Whether to include the contribution from the trend
single_planet : int
Index of a single planet to include in the model, starting at 1.
Use positive values (1, 2, …) for the Np planets and negative
values (-1, -2, …) for the known object planets.
except_planet : int or list
Index (or list of indices) of a single planet to exclude from
the model, starting at 1. Use positive values (1, 2, …) for the
Np planets and negative
values (-1, -2, …) for the known object planets.
def from_prior(self, n=1)
Generate n samples from the priors for all parameters.
def full_model(self, sample, t=None, **kwargs)
Evaluate the full model at one posterior sample, including the GP.
If t is None, use the observed times. Instrument offsets are only
added if t is None, but the systemic velocity is always added.
To evaluate at all posterior samples, consider using
np.apply_along_axis(self.full_model, 1, self.posterior_sample)
Arguments
sample : ndarray
One posterior sample, with shape (npar,)
t : ndarray (optional)
Times at which to evaluate the model, or None to use observed times
**kwargs
Keyword arguments passed directly to eval_model
def get_marginals(self)
Get the marginal posteriors from the posterior_sample matrix.
They go into self.T, self.A, self.E, etc
def get_medians(self)
return the median values of all the parameters
def get_sorted_planet_samples(self, full=True)
def hist_MA(self)
Plot the histogram of the posterior for the MA parameters
def hist_correlations(self)
Plot the histogram of the posterior for the activity correlations
def hist_jitter(self, show_prior=False, **kwargs)
Plot the histogram of the posterior for the additional white noise
def hist_nu(self, show_prior=False, **kwargs)
Plot the histogram of the posterior for the Student-t degrees of freedom
Plot the histogram of the posterior for the systemic velocity and for
the between-instrument offsets (if show_offsets is True and the model
has multiple instruments). If specific is not None, it should be a
tuple with the name of the datafiles for two instruments (matching
self.data_file). In that case, this function works out the RV offset
between the specific[0] and specific[1] instruments.
Plot the histogram (or the kde) of the posterior for the orbital period(s).
Optionally provide the number of histogram bins, the bins themselves, the
limits in orbital period, or the kde bandwidth. If both kde and
show_peaks are true, the routine attempts to plot the most prominent
peaks in the posterior.
Plot the 2d histograms of the posteriors for semi-amplitude and orbital
period and eccentricity and orbital period. If points is True, plot
each posterior sample, else plot hexbins
def make_plot4(self, *args, **kwargs)
Plot histograms for the GP hyperparameters. If Np is not None,
highlight the samples with Np Keplerians.
Get the maximum likelihood sample. By default, this is the highest
likelihood sample found by DNest4. If from_posterior is True, this
returns instead the highest likelihood sample from those that represent
the posterior. The latter may change, due to random choices, between
different calls to "showresults". If Np is given, select only samples
with that number of planets.
def median_sample(self, Np=None, printit=True)
Get the median posterior sample. If Np is given, select only from
posterior samples with that number of planets.
Evaluate the planet part of the model at one posterior sample. If t
is None, use the observed times.
To evaluate at all posterior samples, consider using
np.apply_along_axis(self.planet_model, 1, self.posterior_sample)
Note: this function does not evaluate the GP component of the model
nor the systemic velocity and instrument offsets.
Arguments
sample : ndarray
One posterior sample, with shape (npar,)
t : ndarray
Times at which to evaluate the model, or None to use observed times
include_known_object : bool, default True
Whether to include the contribution from the known object planets
single_planet : int
Index of a single planet to include in the model, starting at 1.
Use positive values (1, 2, …) for the Np planets and negative
values (-1, -2, …) for the known object planets.
except_planet : int or list
Index (or list of indices) of a single planet to exclude from
the model, starting at 1. Use positive values (1, 2, …) for the
Np planets and negative
values (-1, -2, …) for the known object planets.
Plot the histogram (or the kde) of the posterior for the orbital period(s).
Optionally provide the number of histogram bins, the bins themselves, the
limits in orbital period, or the kde bandwidth. If both kde and
show_peaks are true, the routine attempts to plot the most prominent
peaks in the posterior.
Plot the 2d histograms of the posteriors for semi-amplitude and orbital
period and eccentricity and orbital period. If points is True, plot
each posterior sample, else plot hexbins
def plot4(self, *args, **kwargs)
Plot histograms for the GP hyperparameters. If Np is not None,
highlight the samples with Np Keplerians.
Display the RV data together with curves from the posterior predictive.
A total of ncurves random samples are chosen, and the Keplerian
curves are calculated covering 100 + over% of the data timespan.
If the model has a GP component, the prediction is calculated using the
GP hyperparameters for each of the random samples.
def plot_posterior_np(self, **kwargs)
Plot the histogram of the posterior for Np
def plot_random_planets(self, **kwargs)
Display the RV data together with curves from the posterior predictive.
A total of ncurves random samples are chosen, and the Keplerian
curves are calculated covering 100 + over% of the data timespan.
If the model has a GP component, the prediction is calculated using the
GP hyperparameters for each of the random samples.
Evaluate the stochastic part of the model (GP) at one posterior sample.
If t is None, use the observed times. Instrument offsets are only
added if t is None, but the systemic velocity is always added.
To evaluate at all posterior samples, consider using
np.apply_along_axis(self.stochastic_model, 1, self.posterior_sample)
Arguments
sample : ndarray
One posterior sample, with shape (npar,)
t : ndarray (optional)
Times at which to evaluate the model, or None to use observed times
return_std : bool (optional, default False)
Whether to return the st.d. of the predictive
derivative : bool (optional, default False)
Return the first time derivative of the GP prediction instead
Module pykima.analysis
Functions
def column_dynamic_ranges(results)
Return the range of each column in the posterior file
def columns_with_dynamic_range(results)
Return the columns in the posterior file which vary
Calculate detection limits using samples with more than Np planets. By
default, this function uses the value of Np which passes the posterior
probability threshold.
Arguments
star_mass : float or tuple
Stellar mass and optionally its uncertainty [in solar masses].
Np : int
Consider only posterior samples with more than Np planets.
bins : int
Number of bins at which to calculate the detection limits. The period
ranges from the minimum to the maximum orbital period in the posterior.
plot : bool (default True)
Whether to plot the detection limits
semi_amplitude : bool (default False)
Show the detection limits for semi-amplitude, instead of planet mass
semi_major_axis : bool (default False)
Show semi-major axis in the x axis, instead of orbital period
logX : bool (default True)
Should X-axis bins be logarithmic ?
sorted_samples : bool
undoc
return_mask : bool
undoc
remove_nan : bool
remove bins with no counts (no posterior samples)
add_jitter …
show_prior …
show_prior : bool (default False)
If true, color the posterior samples according to orbital eccentricity
smooth : bool (default False)
Smooth the binned maximum with a polynomial
smooth_degree : int (default 3)
Degree of the polynomial used for smoothing
Returns
P, K, E, M : ndarray
Orbital periods, semi-amplitudes, eccentricities, and planet masses used
in the calculation of the detection limits. These correspond to all the
posterior samples with more than Np planets
s : DLresult, namedtuple
Detection limits result, with attributes max and bins. The max
array is in units of Earth masses, bins is in days.
Estimate which observations are outliers, for a model with a Student t
likelihood. This function first calculates the residuals given the
parameters in sample. Then it computes the relative probability of each
residual point given a Student-t (Td) and a Gaussian (Gd) likelihoods. If
the probability Td is larger than Gd (by a factor of threshold), the point
is flagged as an outlier. The function returns an "outlier mask".
A simple context manager to switch directories temporarily
def clipped_mean(arr, min, max)
Mean of arr between min and max
def clipped_std(arr, min, max)
Standard deviation of arr between min and max
def date_to_jd(year, month, day)
Convert a date (year, month, day) to Julian Day.
Algorithm from 'Practical Astronomy with your Calculator or Spreadsheet',
4th ed., Duffet-Smith and Zwart, 2011.
Parameters
year : int
Year as integer. Years preceding 1 A.D. should be 0 or negative.
The year before 1 A.D. is 0, 10 B.C. is year -9.
month : int
Month as integer, Jan = 1, Feb. = 2, etc.
day : float
Day, may contain fractional part.
Returns
jd : float
Julian Day
Examples
Convert 6 a.m., February 17, 1985 to Julian Day
>>> date_to_jd(1985, 2, 17.25)
2446113.75
def datetime_to_jd(date)
Convert a datetime.datetime object to Julian Day.
Parameters
date : datetime.datetime instance
Returns
jd : float
Julian day.
Examples
>>> d = datetime.datetime(1985, 2, 17, 6)
>>> d
datetime.datetime(1985, 2, 17, 6, 0)
>>> jdutil.datetime_to_jd(d)
2446113.75
def days_to_hmsm(days)
Convert fractional days to hours, minutes, seconds, and microseconds.
Precision beyond microseconds is rounded to the nearest microsecond.
Parameters
days : float
A fractional number of days. Must be less than 1.
Returns
hour : int
Hour number.
min : int
Minute number.
sec : int
Second number.
micro : int
Microsecond number.
Raises
ValueError
If days is >= 1.
Examples
>>> days_to_hmsm(0.1)
(2, 24, 0, 0)
def find_prior_limits(prior)
Find lower and upper limits of a prior from the kima_model_setup.txt file.
def find_prior_parameters(prior)
Find parameters of a prior from the kima_model_setup.txt file.
def get_instrument_name(data_file)
Find instrument name
def get_planet_mass(P, K, e, star_mass=1.0, full_output=False, verbose=False)
Calculate the planet (minimum) mass Msini given orbital period P,
semi-amplitude K, eccentricity e, and stellar mass. If star_mass is a
tuple with (estimate, uncertainty), this (Gaussian) uncertainty will be
taken into account in the calculation.
mean Msini [Mjup], std Msini [Mjup], Msini [Mjup] (array)
else: mean Msini [Mjup], std Msini [Mjup], mean Msini [Mearth], std Msini [Mearth]
def get_planet_mass_and_semimajor_axis(P, K, e, star_mass=1.0, full_output=False, verbose=False)
Calculate the planet (minimum) mass Msini and the semi-major axis given
orbital period P, semi-amplitude K, eccentricity e, and stellar mass.
If star_mass is a tuple with (estimate, uncertainty), this (Gaussian)
uncertainty will be taken into account in the calculation.
Units
P [days]
K [m/s]
e []
star_mass [Msun]
Returns
(M, A) where
M is the output of get_planet_mass
A is the output of get_planet_semimajor_axis
def get_planet_mass_latex(P, K, e, star_mass=1.0, earth=False, **kargs)
Tstar = stellar effective temperature in kelvins (K)
Rstar = stellar radius in solar radii
a = distance from the star in AU (change to solRad)
A = bond albedo
f = redistribution factor
Returns equilibrium temperature in Kelvins
def get_prior(prior)
Return a scipt.stats-like prior from a kima_model_setup.txt prior
def get_star_name(data_file)
Find star name (usually works for approx. standard filenames)
Read data from datafile for multiple instruments.
Can be str, in which case the 4th column is assumed to contain an integer
identifier of the instrument.
Or list, in which case each element will be one different filename
containing three columns each.
def read_datafile_rvfwhm(datafile, skip)
Read data from datafile for multiple instruments and RV-FWHM data.
Can be str, in which case the 4th column is assumed to contain an integer
identifier of the instrument.
Or list, in which case each element will be one different filename
containing three columns each.