<jupyter_start><jupyter_text>Gaussian Process Fitting by <NAME> PrerequisitesThis tutorial assumes knowledge of the basic `radvel` API for $\chi^2$ likelihood fitting. As such, please complete the following before beginning this tutorial:- radvel/docs/tutorials/164922_Fitting+MCMC.ipynb- radvel/docs/tutorials/K2-24_Fitting+MCMC.ipynbThis tutorial also assumes knowledge of Gaussian Processes (GPs) as applied to radial velocity (RV) timeseries modeling. Grunblatt et al. (2015) and Rajpaul et al. (2015) contain excellent introductions to this topic. Also check out "Gaussian Processes for Machine Learning," by <NAME>, a free online textbook hosted at gaussianprocesses.org. ObjectivesUsing the K2-131 (EPIC-228732031) dataset published in Dai et al. (2017), I will show how to:- perform a maximum a posteriori (MAP) fit using a quasi-periodic kernel GP regression to model stellar activity (with data from multiple telescopes)- do an MCMC exploration of the corresponding parameter space (with data from multiple telescopes) Tutorial Do some preliminary imports:<jupyter_code>import numpy as np
import pandas as pd
import os
import radvel
import radvel.likelihood
from radvel.plot import orbit_plots, mcmc_plots
from scipy import optimize

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Read in RV data from Dai et al. (2017):<jupyter_code>data = pd.read_csv(os.path.join(radvel.DATADIR,'k2-131.txt'), sep=' ')

t = np.array(data.time)
vel = np.array(data.mnvel)
errvel = np.array(data.errvel)
tel = np.array(data.tel)

telgrps = data.groupby('tel').groups
instnames = telgrps.keys()<jupyter_output><empty_output><jupyter_text>We'll use a quasi-periodic covariance kernel in this fit. An element in the covariance matrix, $C_{ij}$ is defined as follows:$$C_{ij} = \eta_1^2 exp[-\frac{|t_i-t_j|^2}{\eta_2^2} -\frac{sin^2(\pi|t_i-t_j|/\eta_3)}{2\eta_4^2}] $$Several other kernels are implemented in `radvel`. The code for all kernels lives in radvel/gp.py. Check out that file if you'd like to implement a new kernel. Side Note: to see a list of all implemented kernels and examples of possible names for their associated hyperparameters...<jupyter_code>print(radvel.gp.KERNELS)<jupyter_output>{'SqExp': ['gp_length', 'gp_amp'], 'Per': ['gp_per', 'gp_length', 'gp_amp'], 'QuasiPer': ['gp_per', 'gp_perlength', 'gp_explength', 'gp_amp'], 'Celerite': ['gp_B', 'gp_C', 'gp_L', 'gp_Prot']}<jupyter_text>Define the GP hyperparameters we will use in our fit:<jupyter_code>hnames = [
    'gp_amp', # eta_1; GP variability amplitude
    'gp_explength', # eta_2; GP non-periodic characteristic length
    'gp_per', # eta_3; GP variability period
    'gp_perlength', # eta_4; GP periodic characteristic length
]<jupyter_output><empty_output><jupyter_text>Define some numbers (derived from photometry) that we will use in our priors on the GP hyperparameters:<jupyter_code>gp_explength_mean = 9.5*np.sqrt(2.) # sqrt(2)*tau in Dai+ 2017 [days]
gp_explength_unc = 1.0*np.sqrt(2.)
gp_perlength_mean = np.sqrt(1./(2.*3.32)) # sqrt(1/(2*gamma)) in Dai+ 2017
gp_perlength_unc = 0.019
gp_per_mean = 9.64 # T_bar in Dai+ 2017 [days]
gp_per_unc = 0.12

Porb = 0.3693038 # orbital period [days]
Porb_unc = 0.0000091 
Tc = 2457582.9360 # [BJD]
Tc_unc = 0.0011<jupyter_output><empty_output><jupyter_text>Dai et al. (2017) derive the above from photometry (see sect 7.2.1). I'm currently working on implementing joint modeling of RVs & photometry and RVs & activity indicators in `radvel`, so stay tuned if you'd like to use those features! Initialize `radvel.Parameters` object:<jupyter_code>nplanets=1
params = radvel.Parameters(nplanets,basis='per tc secosw sesinw k')<jupyter_output><empty_output><jupyter_text>Set initial guesses for each fitting parameter:<jupyter_code>params['per1'] = radvel.Parameter(value=Porb)
params['tc1'] = radvel.Parameter(value=Tc)
params['sesinw1'] = radvel.Parameter(value=0.,vary=False) # fix eccentricity = 0
params['secosw1'] = radvel.Parameter(value=0.,vary=False)
params['k1'] = radvel.Parameter(value=6.55)
params['dvdt'] = radvel.Parameter(value=0.,vary=False)
params['curv'] = radvel.Parameter(value=0.,vary=False)<jupyter_output><empty_output><jupyter_text>Set initial guesses for GP hyperparameters:<jupyter_code>params['gp_amp'] = radvel.Parameter(value=25.0)
params['gp_explength'] = radvel.Parameter(value=gp_explength_mean) 
params['gp_per'] = radvel.Parameter(value=gp_per_mean) 
params['gp_perlength'] = radvel.Parameter(value=gp_perlength_mean)<jupyter_output><empty_output><jupyter_text>Instantiate a `radvel.model.RVmodel` object, with `radvel.Parameters` object as attribute:<jupyter_code>gpmodel = radvel.model.RVModel(params)<jupyter_output><empty_output><jupyter_text>Initialize `radvel.likelihood.GPLikelihood` objects (one for each telescope):<jupyter_code>jit_guesses = {'harps-n':0.5, 'pfs':5.0}

likes = []
def initialize(tel_suffix):

    # Instantiate a separate likelihood object for each instrument.
    # Each likelihood must use the same radvel.RVModel object.
    indices = telgrps[tel_suffix]
    like = radvel.likelihood.GPLikelihood(gpmodel, t[indices], vel[indices], 
                                          errvel[indices], hnames, suffix='_'+tel_suffix, 
                                          kernel_name="QuasiPer"
                                          )
    # Add in instrument parameters
    like.params['gamma_'+tel_suffix] = radvel.Parameter(value=np.mean(vel[indices]), vary=False, linear=True)
    like.params['jit_'+tel_suffix] = radvel.Parameter(value=jit_guesses[tel_suffix], vary=True)
    likes.append(like)


for tel in instnames:
    initialize(tel)<jupyter_output><empty_output><jupyter_text>Instantiate a `radvel.likelihood.CompositeLikelihood` object that has both GP likelihoods as attributes:<jupyter_code>gplike = radvel.likelihood.CompositeLikelihood(likes)<jupyter_output><empty_output><jupyter_text>Instantiate a `radvel.Posterior` object:<jupyter_code>gppost = radvel.posterior.Posterior(gplike)<jupyter_output><empty_output><jupyter_text>Add in priors (see Dai et al. 2017 section 7.2):<jupyter_code>gppost.priors += [radvel.prior.Gaussian('per1', Porb, Porb_unc)]
gppost.priors += [radvel.prior.Gaussian('tc1', Tc, Tc_unc)]
gppost.priors += [radvel.prior.Jeffreys('k1', 0.01, 10.)] # min and max for Jeffrey's priors estimated by Sarah
gppost.priors += [radvel.prior.Jeffreys('gp_amp', 0.01, 100.)]
gppost.priors += [radvel.prior.Jeffreys('jit_pfs', 0.01, 10.)]
gppost.priors += [radvel.prior.Jeffreys('jit_harps-n', 0.01,10.)]
gppost.priors += [radvel.prior.Gaussian('gp_explength', gp_explength_mean, gp_explength_unc)]
gppost.priors += [radvel.prior.Gaussian('gp_per', gp_per_mean, gp_per_unc)]
gppost.priors += [radvel.prior.Gaussian('gp_perlength', gp_perlength_mean, gp_perlength_unc)]<jupyter_output><empty_output><jupyter_text>Note: our prior on `'gp_perlength'` isn't equivalent to the one Dai et al. (2017) use because our formulations of the quasi-periodic kernel are slightly different. The results aren't really affected. Do a MAP fit:<jupyter_code>res = optimize.minimize(
    gppost.neglogprob_array, gppost.get_vary_params(), method='Nelder-Mead',
    options=dict(maxiter=200, maxfev=100000, xatol=1e-8)
)

print(gppost)<jupyter_output>parameter                     value      vary
per1                       0.369302       True
tc1                     2.45758e+06       True
secosw1                           0      False
sesinw1                           0      False
k1                          6.71302       True
dvdt                              0      False
curv                              0      False
gp_amp                      23.3362       True
gp_explength                13.3721       True
gp_per                      9.62823       True
gp_perlength               0.381962       True
gamma_harps-n              -6695.13      False
jit_harps-n                     0.5      False
gamma_pfs                   1.66129      False
jit_pfs                           5      False

Priors
------
Gaussian prior on per1, mu=0.3693038, sigma=9.1e-06
Gaussian prior on tc1, mu=2457582.936, sigma=0.0011
Jeffrey's prior on k1, min=0.01, max=10.0
Jeffrey's prior on gp_amp, min=0.01, max=100.0
Jeffrey's prior on jit_pfs, min=0.01, max[...]<jupyter_text>Explore the parameter space with MCMC:<jupyter_code>chains = radvel.mcmc(gppost,nrun=100,ensembles=3,savename='rawchains.h5')<jupyter_output>15000/15000 (100.0%) steps complete; Running 1231.70 steps/s; Mean acceptance rate = 45.7%; Min Auto Factor =  12; Max Auto Relative-Change = 0.978; Min Tz = 248.2; Max G-R = 1.012
Discarding burn-in now that the chains are marginally well-mixed



MCMC: WARNING: chains did not pass convergence tests. They are likely not well-mixed.<jupyter_text>Note: for reliable results, run MCMC until the chains have converged. For this example, nrun=10000 should do the trick, but that would take a minute or two, and I won't presume to take up that much of your time with this tutorial.Make some nice plots:<jupyter_code># try switching some of these (optional) keywords to "True" to see what they do!
GPPlot = orbit_plots.GPMultipanelPlot(
    gppost,
    subtract_gp_mean_model=False,
    plot_likelihoods_separately=False, 
    subtract_orbit_model=False
) 
GPPlot.plot_multipanel()
Corner = mcmc_plots.CornerPlot(gppost, chains) # posterior distributions
Corner.plot()
quants = chains.quantile([0.159, 0.5, 0.841]) # median & 1sigma limits of posterior distributions

for par in gppost.params.keys():
    if gppost.params[par].vary:
        med = quants[par][0.5]
        high = quants[par][0.841] - med
        low = med - quants[par][0.159]
        err = np.mean([high,low])
        err = radvel.utils.round_sig(err)
        med, err, errhigh = radvel.utils.sigfig(med, err)
        print('{} : {} +/- {}'.format(par, med, err))<jupyter_output>per1 : 0.3693031 +/- 2.1e-06
tc1 : 2457582.9366 +/- 0.0036
k1 : 6.94 +/- 0.83
gp_amp : 24.7 +/- 3.1
gp_explength : 13.9 +/- 1.1
gp_per : 9.71 +/- 0.24
gp_perlength : 0.39 +/- 0.02
{'fraction_non_alphanumeric': 0.11724275326227712, 'fraction_numerical': 0.04392867815519474, 'mean_word_length': 4.075834175935288, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 63, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '23632', 'license': 'MIT', 'n_tokens_mistral': 3807, 'n_tokens_neox': 3390, 'n_words': 859, 'path': 'docs/tutorials/GaussianProcess-tutorial.ipynb', 'repo_name': 'zdebeurs/radvel', 'size': '10058'}