Fit Gaussian Process
This example shows how to fit a Gaussian process to model stellar activity in RV data. It continues from Basic RV Fit.
Radial velocity modelling is supported in Octofitter via the extension package OctofitterRadialVelocity. To install it, run pkg> add OctofitterRadialVelocity
For this example, we will fit the orbit of the planet K2-131 to perform the same fit as in the RadVel Gaussian Process Fitting tutorial.
We will use the following packages:
using Octofitter
using OctofitterRadialVelocity
using PlanetOrbits
using CairoMakie
using PairPlots
using CSV
using DataFrames
using DistributionsWe will start by downloading and preparing a table of radial velocity measurements, and create a StarAbsoluteRVLikelihood object to hold them.
The following functions allow you to directly load data from various public RV databases:
HARPS_DR1_rvs("star-name")HARPS_RVBank_observations("star-name")Lick_rvs("star-name")HIRES_rvs("star-name")
If you would like to manually specify RV data, use the following format:
rv_like = StarAbsoluteRVLikelihood(
# epoch is in units of MJD. `jd` is a helper function to convert.
# you can also put `years2mjd(2016.1231)`.
# rv and σ_rv are in units of meters/second
(epoch=jd2mjd(2455110.97985), rv=−6.54, σ_rv=1.30),
(epoch=jd2mjd(2455171.90825), rv=−3.33, σ_rv=1.09),
# ... etc.
instrument_name="insert name here",
offset=:rv0_harps,
jitter=:jitter_harps,
)StarAbsoluteRVLikelihood Table with 3 columns and 2 rows:
epoch rv σ_rv
┌─────────────────────
1 │ 55110.5 -6.54 1.3
2 │ 55171.4 -3.33 1.09For this example, to replicate the results of RadVel, we will download their example data for K2-131 and format it for Octofitter:
rv_file = download("https://raw.githubusercontent.com/California-Planet-Search/radvel/master/example_data/k2-131.txt")
rv_dat_raw = CSV.read(rv_file, DataFrame, delim=' ')
rv_dat = DataFrame();
rv_dat.epoch = jd2mjd.(rv_dat_raw.time)
rv_dat.rv = rv_dat_raw.mnvel
rv_dat.σ_rv = rv_dat_raw.errvel
tels = sort(unique(rv_dat_raw.tel))
# This table includes data from two insturments. We create a separate
# likelihood object for each:
rvlike_harps = StarAbsoluteRVLikelihood(
rv_dat[rv_dat_raw.tel .== "harps-n",:],
instrument_name="harps-n",
offset=:rv0_harps,
jitter=:jitter_harps,
)
rvlike_pfs = StarAbsoluteRVLikelihood(
rv_dat[rv_dat_raw.tel .== "pfs",:],
instrument_name="pfs",
offset=:rv0_pfs,
jitter=:jitter_pfs,
)StarAbsoluteRVLikelihood Table with 3 columns and 31 rows:
epoch rv σ_rv
┌──────────────────────
1 │ 57828.4 -45.17 3.88
2 │ 57828.4 -36.47 4.03
3 │ 57829.2 -20.52 3.55
4 │ 57829.2 -4.58 6.07
5 │ 57830.1 17.71 3.56
6 │ 57830.1 23.34 3.48
7 │ 57830.2 21.34 3.51
8 │ 57830.2 19.92 3.26
9 │ 57830.2 -0.74 3.44
10 │ 57830.3 16.9 3.5
11 │ 57830.3 8.25 4.64
12 │ 57830.4 -4.06 4.91
13 │ 57832.1 -29.51 4.96
14 │ 57832.1 -21.56 4.11
15 │ 57832.2 -25.76 5.59
16 │ 57833.1 2.1 3.91
17 │ 57833.1 4.38 3.73
⋮ │ ⋮ ⋮ ⋮Gaussian Process Fit
Let us now add a Gaussian process to model stellar activity. This should improve the fit.
We start by writing a function that creates a Gaussian process kernel from a set of system parameters. We will create a quasi-periodic kernel.
using AbstractGPs
function gauss_proc_quasi_periodic(θ_system)
# θ_system is a named tuple of variable values for this posterior draw.
# Note: for multiple instruments, you can supply different gaussian process
# functions that e.g. use different hyper-parameter names.
η₁ = θ_system.gp_η₁ # h
η₂ = θ_system.gp_η₂ # λ
η₃ = θ_system.gp_η₃ # θ
η₄ = θ_system.gp_η₄ # ω
kernel = η₁^2 *
(SqExponentialKernel() ∘ ScaleTransform(1/(η₂))) *
(PeriodicKernel(r=[η₄]) ∘ ScaleTransform(1/(η₃)))
gp = GP(kernel)
return gp
endgauss_proc_quasi_periodic (generic function with 1 method)Now provide this function when we construct our likelihood object:
rvlike_harps = StarAbsoluteRVLikelihood(
rv_dat[rv_dat_raw.tel .== "harps-n",:],
instrument_name="harps-n",
offset=:rv0_harps,
jitter=:jitter_harps,
# NEW:
gaussian_process = gauss_proc_quasi_periodic
)
rvlike_pfs = StarAbsoluteRVLikelihood(
rv_dat[rv_dat_raw.tel .== "pfs",:],
instrument_name="pfs",
offset=:rv0_pfs,
jitter=:jitter_pfs,
# NEW:
gaussian_process = gauss_proc_quasi_periodic
)StarAbsoluteRVLikelihood Table with 3 columns and 31 rows:
epoch rv σ_rv
┌──────────────────────
1 │ 57828.4 -45.17 3.88
2 │ 57828.4 -36.47 4.03
3 │ 57829.2 -20.52 3.55
4 │ 57829.2 -4.58 6.07
5 │ 57830.1 17.71 3.56
6 │ 57830.1 23.34 3.48
7 │ 57830.2 21.34 3.51
8 │ 57830.2 19.92 3.26
9 │ 57830.2 -0.74 3.44
10 │ 57830.3 16.9 3.5
11 │ 57830.3 8.25 4.64
12 │ 57830.4 -4.06 4.91
13 │ 57832.1 -29.51 4.96
14 │ 57832.1 -21.56 4.11
15 │ 57832.2 -25.76 5.59
16 │ 57833.1 2.1 3.91
17 │ 57833.1 4.38 3.73
⋮ │ ⋮ ⋮ ⋮Note that the two instruments do not need to use the same Gaussian process kernels, nor the same hyper parameter names.
We create a new system model with added parameters gp_η₁ etc for the Gaussian process fit.
gp_explength_mean = 9.5*sqrt(2.) # sqrt(2)*tau in Dai+ 2017 [days]
gp_explength_unc = 1.0*sqrt(2.)
gp_perlength_mean = 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
# No change to planet model
@planet b RadialVelocityOrbit begin
e = 0
ω = 0.0
P ~ truncated(Normal(0.3693038/365.25, 0.0000091/365.25),lower=0.00001)
a = cbrt(system.M * b.P^2)
mass ~ LogUniform(0.001, 10)
τ ~ UniformCircular(1.0)
tp = b.τ*b.P*365.25 + 57782 # reference epoch for τ. Choose an MJD date near your data.
end
@system k2_132_gp begin
M ~ truncated(Normal(0.82, 0.02),lower=0.1) # (Baines & Armstrong 2011).
rv0_harps ~ Normal(-6693,100) # m/s
rv0_pfs ~ Normal(0,100) # m/s
jitter_harps ~ LogUniform(0.1,100) # m/s
jitter_pfs ~ LogUniform(0.1,100) # m/s
# Add priors on GP kernel hyper-parameters.
gp_η₁ ~ truncated(Normal(25,10),lower=0.1,upper=100)
# Important: ensure the period and exponential length scales
# have physically plausible lower and upper limits to avoid poor numerical conditioning
gp_η₂ ~ truncated(Normal(gp_explength_mean,gp_explength_unc),lower=5,upper=100)
gp_η₃ ~ truncated(Normal(gp_per_mean,1),lower=2, upper=100)
gp_η₄ ~ truncated(Normal(gp_perlength_mean,gp_perlength_unc),lower=0.2, upper=10)
end rvlike_harps rvlike_pfs b
model = Octofitter.LogDensityModel(k2_132_gp)LogDensityModel for System k2_132_gp of dimension 13 and 71 epochs with fields .ℓπcallback and .∇ℓπcallback
Sample from the model as before:
# Seed the random number generator
using Random
rng = Random.Xoshiro(0)
chain = octofit(
rng, model,
adaptation = 100,
iterations = 100,
)Chains MCMC chain (100×31×1 Array{Float64, 3}):
Iterations = 1:1:100
Number of chains = 1
Samples per chain = 100
Wall duration = 14.79 seconds
Compute duration = 14.79 seconds
parameters = M, rv0_harps, rv0_pfs, jitter_harps, jitter_pfs, gp_η₁, gp_η₂, gp_η₃, gp_η₄, b_P, b_mass, b_τy, b_τx, b_e, b_ω, b_a, b_τ, b_tp
internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rha ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float6 ⋯
M 0.8168 0.0184 0.0023 64.6491 101.7705 1.001 ⋯
rv0_harps -6695.6863 13.4018 1.4315 92.1456 52.4503 1.024 ⋯
rv0_pfs -6.9228 13.8633 1.5789 79.7984 58.9391 1.015 ⋯
jitter_harps 0.8189 0.7936 0.0916 36.7594 13.0721 1.001 ⋯
jitter_pfs 4.8293 1.2241 0.1089 138.5005 58.9391 0.997 ⋯
gp_η₁ 26.7013 5.1117 0.4364 128.5726 113.4378 1.034 ⋯
gp_η₂ 13.4843 1.1682 0.1892 46.8912 40.5549 1.029 ⋯
gp_η₃ 9.2630 0.8327 0.1035 67.2155 49.8569 0.994 ⋯
gp_η₄ 0.3844 0.0181 0.0039 27.9357 37.3971 1.018 ⋯
b_P 0.0010 0.0000 0.0000 91.8962 61.5167 1.023 ⋯
b_mass 0.0207 0.0040 0.0006 46.5259 78.3393 1.050 ⋯
b_τy 0.6285 0.1624 0.0146 122.8012 87.2420 1.002 ⋯
b_τx -0.7769 0.1446 0.0140 103.3247 78.3393 0.995 ⋯
b_e 0.0000 0.0000 NaN NaN NaN Na ⋯
b_ω 0.0000 0.0000 NaN NaN NaN Na ⋯
b_a 0.0094 0.0001 0.0000 64.5699 101.7705 1.001 ⋯
b_τ -0.1417 0.0288 0.0026 133.1477 76.6284 1.010 ⋯
b_tp 57781.9477 0.0106 0.0009 133.1477 76.6284 1.010 ⋯
2 columns omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5 ⋯
Symbol Float64 Float64 Float64 Float64 Float6 ⋯
M 0.7838 0.8040 0.8151 0.8272 0.857 ⋯
rv0_harps -6718.4882 -6703.8768 -6696.9319 -6686.7020 -6670.425 ⋯
rv0_pfs -46.2721 -13.2023 -6.2129 2.8742 14.868 ⋯
jitter_harps 0.1041 0.2180 0.4668 1.3510 2.533 ⋯
jitter_pfs 2.6906 4.0536 4.6648 5.4337 7.561 ⋯
gp_η₁ 17.9144 22.9258 26.3059 30.0721 37.090 ⋯
gp_η₂ 11.4765 12.6244 13.3660 14.1361 16.207 ⋯
gp_η₃ 7.5370 8.7687 9.3406 9.7008 10.971 ⋯
gp_η₄ 0.3494 0.3739 0.3850 0.3970 0.415 ⋯
b_P 0.0010 0.0010 0.0010 0.0010 0.001 ⋯
b_mass 0.0135 0.0180 0.0203 0.0239 0.027 ⋯
b_τy 0.3301 0.5299 0.6177 0.7357 0.948 ⋯
b_τx -1.0802 -0.8702 -0.7704 -0.6704 -0.536 ⋯
b_e 0.0000 0.0000 0.0000 0.0000 0.000 ⋯
b_ω 0.0000 0.0000 0.0000 0.0000 0.000 ⋯
b_a 0.0093 0.0094 0.0094 0.0095 0.009 ⋯
b_τ -0.1910 -0.1608 -0.1407 -0.1206 -0.088 ⋯
b_tp 57781.9295 57781.9406 57781.9480 57781.9555 57781.967 ⋯
1 column omitted
For real data, we would want to increase the adaptation and iterations to about 1000 each.
Plot one sample from the results:
fig = Octofitter.rvpostplot(model, chain) # saved to "k2_132-rvpostplot.png"
Plot many samples from the results:
fig = octoplot(
model,
chain,
# Some optional tweaks to the appearance:
N=50, # only plot 50 samples
figscale=1.5, # make it larger
alpha=0.05, # make each sample more transparent
colormap="#0072b2",
) # saved to "k2_132-plot-grid.png"
Corner plot:
octocorner(model, chain, small=true) # saved to "k2_132-pairplot-small.png"