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
There are two different GP packages supported by OctofitterRadialVelocity: AbstractGPs, and Celerite. Important note: Celerite.jl does not support Julia 1.0+, so we currently bundle a fork that has been patched to work. When / if Celerite.jl is updated we will switch back to the public package.
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 pick up from our tutorial Basic RV Fit with the data already downloaded and available as a table called rv_dat:
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))2-element Vector{InlineStrings.String7}:
"harps-n"
"pfs"Gaussian Process Fit with AbstractGPs
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. We provide this function as an arugment gaussian_process to the likelihood constructor:
using AbstractGPs
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
rvlike_harps = StarAbsoluteRVLikelihood(
rv_dat[rv_dat_raw.tel .== "harps-n",:],
name="harps-n",
variables=(@variables begin
offset ~ Normal(-6693,100) # m/s
jitter ~ LogUniform(0.1,100) # m/s
# Add priors on GP kernel hyper-parameters.
η_1 ~ 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
η_2 ~ truncated(Normal(gp_explength_mean,gp_explength_unc),lower=5,upper=100)
η_3 ~ truncated(Normal(gp_per_mean,1),lower=2, upper=100)
η_4 ~ truncated(Normal(gp_perlength_mean,gp_perlength_unc),lower=0.2, upper=10)
end),
gaussian_process = θ_obs -> GP(
θ_obs.η_1^2 *
(SqExponentialKernel() ∘ ScaleTransform(1/(θ_obs.η_2))) *
(PeriodicKernel(r=[θ_obs.η_4]) ∘ ScaleTransform(1/(θ_obs.η_3)))
)
)
rvlike_pfs = StarAbsoluteRVLikelihood(
rv_dat[rv_dat_raw.tel .== "pfs",:],
name="pfs",
variables=(@variables begin
offset ~ Normal(0,100) # m/s
jitter ~ LogUniform(0.1,100) # m/s
# Add priors on GP kernel hyper-parameters.
η_1 ~ 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
η_2 ~ truncated(Normal(gp_explength_mean,gp_explength_unc),lower=5,upper=100)
η_3 ~ truncated(Normal(gp_per_mean,1),lower=2, upper=100)
η_4 ~ truncated(Normal(gp_perlength_mean,gp_perlength_unc),lower=0.2, upper=10)
end),
gaussian_process = θ_obs -> GP(
θ_obs.η_1^2 *
(SqExponentialKernel() ∘ ScaleTransform(1/(θ_obs.η_2))) *
(PeriodicKernel(r=[θ_obs.η_4]) ∘ ScaleTransform(1/(θ_obs.η_3)))
)
)
## No change to the rest of the model
planet_1 = Planet(
name="b",
basis=RadialVelocityOrbit,
likelihoods=[],
variables=@variables begin
e = 0
ω = 0.0
# To match RadVel, we set a prior on Period and calculate semi-major axis from it
P ~ truncated(
Normal(0.3693038/365.256360417, 0.0000091/365.256360417),
lower=0.0001
)
M = system.M
a = cbrt(M * P^2) # note the equals sign.
τ ~ UniformCircular(1.0)
tp = τ*P*365.256360417 + 57782 # reference epoch for τ. Choose an MJD date near your data.
# minimum planet mass [jupiter masses]. really m*sin(i)
mass ~ LogUniform(0.001, 10)
end
)
sys = System(
name = "k2_132",
companions=[planet_1],
likelihoods=[rvlike_harps, rvlike_pfs],
variables=@variables begin
M ~ truncated(Normal(0.82, 0.02),lower=0.1) # (Baines & Armstrong 2011).
end
)
model = Octofitter.LogDensityModel(sys)LogDensityModel for System k2_132 of dimension 17 and 71 epochs with fields .ℓπcallback and .∇ℓπcallback
Note that the two instruments do not need to use the same Gaussian process kernels, nor the same hyper parameter names.
Tip: If you want the instruments to share the Gaussian process kernel hyper parameters, move the variables up to the system's @variables block, and forward them to the observation variables block e.g. η₁ = system.η₁, η₂ = system.η₂.
Initialize the starting points, and confirm the data are entered correcly:
init_chain = initialize!(model)
fig = Octofitter.rvpostplot(model, init_chain)
Sample from the model using MCMC (the no U-turn sampler)
# Seed the random number generator
using Random
rng = Random.Xoshiro(0)
chain = octofit(
rng, model,
adaptation = 100,
iterations = 100,
)Chains MCMC chain (100×36×1 Array{Float64, 3}):
Iterations = 1:1:100
Number of chains = 1
Samples per chain = 100
Wall duration = 447.28 seconds
Compute duration = 447.28 seconds
parameters = M, harps_n_offset, harps_n_jitter, harps_n_η_1, harps_n_η_2, harps_n_η_3, harps_n_η_4, pfs_offset, pfs_jitter, pfs_η_1, pfs_η_2, pfs_η_3, pfs_η_4, b_P, b_τx, b_τy, b_mass, b_τ, b_e, b_ω, b_M, b_a, 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 r ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Floa ⋯
M 0.8200 0.0195 0.0020 80.6685 73.5692 0.9 ⋯
harps_n_offset -6689.2968 9.0992 4.1649 4.4022 13.2103 1.1 ⋯
harps_n_jitter 0.9073 0.8489 0.1052 56.5264 77.3213 1.0 ⋯
harps_n_η_1 25.8828 5.8889 0.6111 86.7155 42.4389 1.0 ⋯
harps_n_η_2 13.3750 1.3068 0.1345 96.2553 76.6284 0.9 ⋯
harps_n_η_3 9.4095 0.9327 0.0821 123.4904 95.2535 1.0 ⋯
harps_n_η_4 0.3863 0.0180 0.0017 116.5122 101.7705 0.9 ⋯
pfs_offset 11.8753 16.0809 8.4468 3.8872 17.3473 1.2 ⋯
pfs_jitter 4.1119 1.8930 0.5793 14.2941 30.2536 1.0 ⋯
pfs_η_1 31.3284 6.7905 1.3463 28.9818 77.3213 1.0 ⋯
pfs_η_2 13.3553 1.4184 0.1614 76.8907 52.4503 0.9 ⋯
pfs_η_3 9.5102 0.9957 0.0711 200.0000 80.9764 0.9 ⋯
pfs_η_4 0.3847 0.0206 0.0017 136.8723 100.0590 0.9 ⋯
b_P 0.0010 0.0000 0.0000 168.5997 92.5822 1.0 ⋯
b_τx 0.6422 0.1540 0.0184 69.9936 76.6284 1.0 ⋯
b_τy -0.7364 0.1516 0.0190 63.0404 76.6284 1.0 ⋯
b_mass 0.0203 0.0046 0.0005 85.5129 76.6284 0.9 ⋯
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱
2 columns and 6 rows omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97 ⋯
Symbol Float64 Float64 Float64 Float64 Floa ⋯
M 0.7839 0.8050 0.8214 0.8327 0.8 ⋯
harps_n_offset -6702.9539 -6698.0365 -6688.7003 -6680.4262 -6676.2 ⋯
harps_n_jitter 0.1637 0.2750 0.4756 1.3022 2.8 ⋯
harps_n_η_1 15.8536 21.2890 25.2474 29.2565 38.7 ⋯
harps_n_η_2 10.5919 12.5289 13.3706 14.0864 15.9 ⋯
harps_n_η_3 7.6756 8.8516 9.4949 10.0039 11.1 ⋯
harps_n_η_4 0.3514 0.3747 0.3869 0.3991 0.4 ⋯
pfs_offset -16.0804 1.5391 10.9045 21.4932 42.1 ⋯
pfs_jitter 0.1145 3.4169 4.3895 5.2835 7.0 ⋯
pfs_η_1 20.2538 26.1074 31.0630 35.9610 44.6 ⋯
pfs_η_2 10.7434 12.5543 13.3358 14.2621 16.2 ⋯
pfs_η_3 7.6528 8.7906 9.5909 10.3066 11.2 ⋯
pfs_η_4 0.3449 0.3706 0.3854 0.3980 0.4 ⋯
b_P 0.0010 0.0010 0.0010 0.0010 0.0 ⋯
b_τx 0.3814 0.5235 0.6280 0.7225 0.9 ⋯
b_τy -1.0000 -0.8462 -0.7318 -0.6213 -0.4 ⋯
b_mass 0.0120 0.0172 0.0209 0.0229 0.0 ⋯
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱
1 column and 6 rows 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"
Gaussian Process Fit with Celerite
We now demonstrate an approximate quasi-static kernel implemented using Celerite. For the class of kernels supported by Celerite, the performance scales much better with the number of data points. This makes it a good choice for modelling large RV datasets.
Make sure that you type using OctofitterRadialVelocity.Celerite and not using Celerite. Celerite.jl does not support Julia 1.0+, so we currently bundle a fork that has been patched to work. When / if Celerite.jl is updated we will switch back to the public package.
using OctofitterRadialVelocity.Celerite
rvlike_harps = StarAbsoluteRVLikelihood(
rv_dat[rv_dat_raw.tel .== "harps-n",:],
name="harps-n",
variables=(@variables begin
offset ~ Normal(-6693,100) # m/s
jitter ~ LogUniform(0.1,100) # m/s
# Add priors on GP kernel hyper-parameters.
B ~ Uniform(0.00001, 2000000)
C ~ Uniform(0.00001, 200)
L ~ Uniform(2, 200)
Prot ~ Uniform(8.5, 20)#Uniform(0, 20)
end),
gaussian_process = θ_obs -> Celerite.CeleriteGP(
Celerite.RealTerm(
#=log_a=# log(θ_obs.B*(1+θ_obs.C)/(2+θ_obs.C)),
#=log_c=# log(1/θ_obs.L)
) + Celerite.ComplexTerm(
#=log_a=# log(θ_obs.B/(2+θ_obs.C)),
#=log_b=# -Inf,
#=log_c=# log(1/θ_obs.L),
#=log_d=# log(2pi/θ_obs.Prot)
)
)
)
rvlike_pfs = StarAbsoluteRVLikelihood(
rv_dat[rv_dat_raw.tel .== "pfs",:],
name="pfs",
variables=(@variables begin
offset ~ Normal(0,100) # m/s
jitter ~ LogUniform(0.1,100) # m/s
# Add priors on GP kernel hyper-parameters.
B ~ Uniform(0.00001, 2000000)
C ~ Uniform(0.00001, 200)
L ~ Uniform(2, 200)
Prot ~ Uniform(8.5, 20)#Uniform(0, 20)
end),
gaussian_process = θ_obs -> Celerite.CeleriteGP(
Celerite.RealTerm(
#=log_a=# log(θ_obs.B*(1+θ_obs.C)/(2+θ_obs.C)),
#=log_c=# log(1/θ_obs.L)
) + Celerite.ComplexTerm(
#=log_a=# log(θ_obs.B/(2+θ_obs.C)),
#=log_b=# -Inf,
#=log_c=# log(1/θ_obs.L),
#=log_d=# log(2pi/θ_obs.Prot)
)
)
)
## No change to the rest of the model
planet_1 = Planet(
name="b",
basis=RadialVelocityOrbit,
likelihoods=[],
variables=@variables begin
e = 0
ω = 0.0
# To match RadVel, we set a prior on Period and calculate semi-major axis from it
P ~ truncated(
Normal(0.3693038/365.256360417, 0.0000091/365.256360417),
lower=0.0001
)
M = system.M
a = cbrt(M * P^2) # note the equals sign.
τ ~ UniformCircular(1.0)
tp = τ*P*365.256360417 + 57782 # reference epoch for τ. Choose an MJD date near your data.
# minimum planet mass [jupiter masses]. really m*sin(i)
mass ~ LogUniform(0.001, 10)
end
)
sys = System(
name = "k2_132",
companions=[planet_1],
likelihoods=[rvlike_harps, rvlike_pfs],
variables=@variables begin
M ~ truncated(Normal(0.82, 0.02),lower=0.1) # (Baines & Armstrong 2011).
end
)
using DifferentiationInterface
using FiniteDiff
model = Octofitter.LogDensityModel(sys, autodiff=AutoFiniteDiff())LogDensityModel for System k2_132 of dimension 17 and 71 epochs with fields .ℓπcallback and .∇ℓπcallback
The Celerite implementation doesn't support our default autodiff-backend (ForwardDiff.jl), so we disable autodiff by setting it to finite differences, and then using the Pigeons slice sampler which doesn't require gradients or (B) use Enzyme autodiff,
Initialize the starting points, and confirm the data are entered correcly:
init_chain = initialize!(model)
fig = Octofitter.rvpostplot(model, init_chain)
using Pigeons
chain, pt = octofit_pigeons(model, n_rounds=7)
fig = Octofitter.rvpostplot(model, chain)