Fit Radial Velocity
You can use Octofitter to fit radial velocity data, either alone or in combination with other kinds of data. Multiple instruments (up to 10) are supported, as are gaussian processes to model stellar activity..
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")
Basic Fit
Start by fitting a 1-planet Keplerian model.
Calling those functions with the name of a star will return a StarAbsoluteRVLikelihood table.
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
(;inst_idx=1, epoch=jd(2455110.97985), rv=−6.54, σ_rv=1.30),
(;inst_idx=1, epoch=jd(2455171.90825), rv=−3.33, σ_rv=1.09),
# ... etc.
)For 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))
rv_dat.inst_idx = map(rv_dat_raw.tel) do tel
findfirst(==(tel), tels)
end
rvlike = StarAbsoluteRVLikelihood(
rv_dat,
instrument_names=["harps-n", "psf"],
)StarAbsoluteRVLikelihood Table with 4 columns and 70×1 rows:
inst_idx epoch rv σ_rv
┌───────────────────────────────────
1,1 │ 1 57782.2 -6682.78 3.71
2,1 │ 1 57783.1 -6710.55 5.95
3,1 │ 1 57783.2 -6698.96 8.76
4,1 │ 1 57812.1 -6672.32 4.0
5,1 │ 1 57812.2 -6672.0 4.22
6,1 │ 1 57812.2 -6690.31 4.63
7,1 │ 1 57813.0 -6701.27 4.12
8,1 │ 1 57813.1 -6700.71 5.34
9,1 │ 1 57813.1 -6695.97 3.67
10,1 │ 1 57813.1 -6705.87 4.23
11,1 │ 1 57813.1 -6694.27 5.12
12,1 │ 1 57813.2 -6699.99 5.0
13,1 │ 1 57813.2 -6696.17 10.43
14,1 │ 1 57836.0 -6675.4 7.03
15,1 │ 1 57836.0 -6665.92 7.48
16,1 │ 1 57836.0 -6661.9 6.03
17,1 │ 1 57836.1 -6657.92 5.08
⋮ │ ⋮ ⋮ ⋮ ⋮Now, create a planet. We can use the RadialVelocityOrbit type from PlanetOrbits.jl that requires fewer parameters (eg no inclination or longitude of ascending node). We could instead use a Visual{KepOrbit} or similar if we wanted to include these parameters and visualize the orbit in the plane of the sky.
@planet b RadialVelocityOrbit 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.00001)
a = cbrt(system.M * b.P^2) # note the equals sign.
τ ~ UniformCircular(1.0)
tp = b.τ*b.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
@system k2_132 begin
# total mass [solar masses]
M ~ truncated(Normal(0.82, 0.02),lower=0) # (Baines & Armstrong 2011).
rv0 ~ Product([
# HARPS-N
Normal(-6693,100), # m/s
# FPS
Normal(0,100), # m/s
])
jitter ~ Product([
LogUniform(0.1,100), # m/s
LogUniform(0.1,100), # m/s
])
end rvlike bSystem model k2_132
Derived:
Priors:
M Truncated(Distributions.Normal{Float64}(μ=0.82, σ=0.02); lower=0.0)
rv0 Distributions.Product{Distributions.Continuous, Distributions.Normal{Float64}, Vector{Distributions.Normal{Float64}}}(v=Distributions.Normal{Float64}[Distributions.Normal{Float64}(μ=-6693.0, σ=100.0), Distributions.Normal{Float64}(μ=0.0, σ=100.0)])
jitter Distributions.Product{Distributions.Continuous, Distributions.LogUniform{Float64}, Vector{Distributions.LogUniform{Float64}}}(v=Distributions.LogUniform{Float64}[Distributions.LogUniform{Float64}(a=0.1, b=100.0), Distributions.LogUniform{Float64}(a=0.1, b=100.0)])
Planet b
Derived:
e, ω, a, τ, tp,
Priors:
P Truncated(Distributions.Normal{Float64}(μ=0.001011081092683449, σ=2.4914008313533153e-8); lower=1.0e-5)
τy Distributions.Normal{Float64}(μ=0.0, σ=1.0)
τx Distributions.Normal{Float64}(μ=0.0, σ=1.0)
mass Distributions.LogUniform{Float64}(a=0.001, b=10.0)
Octofitter.UnitLengthPrior{:τy, :τx}: √(τy^2+τx^2) ~ LogNormal(log(1), 0.02)
StarAbsoluteRVLikelihood Table with 4 columns and 70×1 rows:
inst_idx epoch rv σ_rv
┌───────────────────────────────────
1,1 │ 1 57782.2 -6682.78 3.71
2,1 │ 1 57783.1 -6710.55 5.95
3,1 │ 1 57783.2 -6698.96 8.76
4,1 │ 1 57812.1 -6672.32 4.0
5,1 │ 1 57812.2 -6672.0 4.22
6,1 │ 1 57812.2 -6690.31 4.63
7,1 │ 1 57813.0 -6701.27 4.12
8,1 │ 1 57813.1 -6700.71 5.34
9,1 │ 1 57813.1 -6695.97 3.67
10,1 │ 1 57813.1 -6705.87 4.23
11,1 │ 1 57813.1 -6694.27 5.12
12,1 │ 1 57813.2 -6699.99 5.0
13,1 │ 1 57813.2 -6696.17 10.43
14,1 │ 1 57836.0 -6675.4 7.03
15,1 │ 1 57836.0 -6665.92 7.48
16,1 │ 1 57836.0 -6661.9 6.03
17,1 │ 1 57836.1 -6657.92 5.08
⋮ │ ⋮ ⋮ ⋮ ⋮
Note how the rvlike object was attached to the k2_132 system instead of the planet. This is because the observed radial velocity is of the star, and is caused by any/all orbiting planets.
The rv0 and jitter parameters specify priors for the instrument-specific offset and white noise jitter standard deviation. The _i index matches the inst_idx used to create the observation table.
Note also here that the mass variable is really msini, or the minimum mass of the planet.
We can now prepare our model for sampling.
model = Octofitter.LogDensityModel(k2_132)LogDensityModel for System k2_132 of dimension 9 with fields .ℓπcallback and .∇ℓπcallback
Sample:
using Random
rng = Random.Xoshiro(0)
chain = octofit(rng, model)Chains MCMC chain (1000×27×1 Array{Float64, 3}):
Iterations = 1:1:1000
Number of chains = 1
Samples per chain = 1000
Wall duration = 37.09 seconds
Compute duration = 37.09 seconds
parameters = M, rv0_1, rv0_2, jitter_1, jitter_2, b_P, b_τy, b_τx, b_mass, 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 rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
M 0.8193 0.0195 0.0006 1285.5216 740.0623 1.0001 ⋯
rv0_1 -6693.6339 2.5011 0.0763 1064.9995 777.6055 1.0074 ⋯
rv0_2 1.5283 3.3251 0.1098 921.3913 636.5256 1.0007 ⋯
jitter_1 14.0045 2.0882 0.0672 1079.8779 636.8365 1.0014 ⋯
jitter_2 17.6918 2.6934 0.0910 884.5163 743.7114 0.9996 ⋯
b_P 0.0010 0.0000 0.0000 1060.7258 535.9508 1.0013 ⋯
b_τy 0.7292 0.1759 0.0063 875.6661 490.2580 1.0025 ⋯
b_τx -0.6545 0.1836 0.0070 738.0436 567.3051 1.0010 ⋯
b_mass 0.0412 0.0091 0.0003 1015.0473 787.3716 0.9990 ⋯
b_e 0.0000 0.0000 NaN NaN NaN NaN ⋯
b_ω 0.0000 0.0000 NaN NaN NaN NaN ⋯
b_a 0.0094 0.0001 0.0000 1285.4258 740.0623 1.0000 ⋯
b_τ -0.1165 0.0374 0.0015 650.5123 403.0202 1.0095 ⋯
b_tp 57781.9570 0.0138 0.0006 650.5322 403.0202 1.0096 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 0.7816 0.8061 0.8190 0.8321 0.8596
rv0_1 -6698.6984 -6695.2325 -6693.4983 -6691.8760 -6689.1283
rv0_2 -4.7429 -0.6659 1.5020 3.7206 8.3080
jitter_1 10.7305 12.5435 13.8325 15.1425 18.7317
jitter_2 13.1840 15.8172 17.4754 19.1940 23.8477
b_P 0.0010 0.0010 0.0010 0.0010 0.0010
b_τy 0.3444 0.6170 0.7449 0.8426 1.0368
b_τx -0.9755 -0.7764 -0.6609 -0.5417 -0.2472
b_mass 0.0246 0.0349 0.0416 0.0477 0.0586
b_e 0.0000 0.0000 0.0000 0.0000 0.0000
b_ω 0.0000 0.0000 0.0000 0.0000 0.0000
b_a 0.0093 0.0094 0.0094 0.0095 0.0096
b_τ -0.1908 -0.1407 -0.1155 -0.0934 -0.0414
b_tp 57781.9295 57781.9481 57781.9573 57781.9655 57781.9847
Excellent! Let's plot the maximum likelihood orbit:
using CairoMakie: Makie
fig = OctofitterRadialVelocity.rvpostplot(model, chain) # saved to "k2_132-rvpostplot.png"
We can also plot a number of samples from the posterior:
using CairoMakie: Makie
octoplot(model, chain)
Create a corner plot:
using PairPlots
using CairoMakie: Makie
octocorner(model, chain, small=true)
Gaussian Process Fit
Let us now add a Gaussian process to model stellar activity. This should improve the fit.
We start by writing a functino 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.
η₁ = θ_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 = StarAbsoluteRVLikelihood(rv_dat,
instrument_names=["harps-n", "psf"],
gaussian_process = gauss_proc_quasi_periodic
)StarAbsoluteRVLikelihood Table with 4 columns and 70×1 rows:
inst_idx epoch rv σ_rv
┌───────────────────────────────────
1,1 │ 1 57782.2 -6682.78 3.71
2,1 │ 1 57783.1 -6710.55 5.95
3,1 │ 1 57783.2 -6698.96 8.76
4,1 │ 1 57812.1 -6672.32 4.0
5,1 │ 1 57812.2 -6672.0 4.22
6,1 │ 1 57812.2 -6690.31 4.63
7,1 │ 1 57813.0 -6701.27 4.12
8,1 │ 1 57813.1 -6700.71 5.34
9,1 │ 1 57813.1 -6695.97 3.67
10,1 │ 1 57813.1 -6705.87 4.23
11,1 │ 1 57813.1 -6694.27 5.12
12,1 │ 1 57813.2 -6699.99 5.0
13,1 │ 1 57813.2 -6696.17 10.43
14,1 │ 1 57836.0 -6675.4 7.03
15,1 │ 1 57836.0 -6665.92 7.48
16,1 │ 1 57836.0 -6661.9 6.03
17,1 │ 1 57836.1 -6657.92 5.08
⋮ │ ⋮ ⋮ ⋮ ⋮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 begin
M ~ truncated(Normal(0.82, 0.02),lower=0) # (Baines & Armstrong 2011).
# HARPS-N
rv0 ~ Product([
# HARPS-N
Normal(0,10000), # m/s
# FPS
Normal(0,10000), # m/s
])
jitter ~ Product([
LogUniform(0.01,100), # m/s
LogUniform(0.01,100), # m/s
])
# Add priors on GP kernel hyper-parameters.
gp_η₁ ~ truncated(Normal(25,10),lower=0.1)
gp_η₂ ~ truncated(Normal(gp_explength_mean,gp_explength_unc),lower=0.1)
gp_η₃ ~ truncated(Normal(gp_per_mean,1),lower=0.1)
gp_η₄ ~ truncated(Normal(gp_perlength_mean,gp_perlength_unc),lower=0)
end rvlike b
model = Octofitter.LogDensityModel(k2_132; autodiff=:ForwardDiff)LogDensityModel for System k2_132 of dimension 13 with fields .ℓπcallback and .∇ℓπcallback
Sample from the model as before:
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 = 1357.7 seconds
Compute duration = 1357.7 seconds
parameters = M, rv0_1, rv0_2, jitter_1, jitter_2, 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 rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
M 0.8200 0.0190 0.0022 70.6228 78.3393 1.0055 ⋯
rv0_1 -6696.1039 9.0837 0.8283 119.6675 112.1358 1.0368 ⋯
rv0_2 -9.9709 9.5439 1.1587 71.4481 41.0413 1.1092 ⋯
jitter_1 0.3715 0.5176 0.0543 68.7626 78.3393 0.9943 ⋯
jitter_2 4.8103 1.0415 0.1222 77.3964 17.2397 1.0297 ⋯
gp_η₁ 24.8792 3.5212 0.3283 109.8340 69.6192 0.9925 ⋯
gp_η₂ 13.6153 1.2569 0.1441 77.7819 58.9391 1.0091 ⋯
gp_η₃ 8.8543 0.8211 0.1673 25.7362 31.8676 1.0008 ⋯
gp_η₄ 0.3820 0.0156 0.0027 36.3375 33.2304 1.0178 ⋯
b_P 0.0010 0.0000 0.0000 71.1714 76.6284 1.0656 ⋯
b_mass 0.0213 0.0030 0.0003 98.3435 78.3393 1.0091 ⋯
b_τy 0.6208 0.1026 0.0105 92.0745 74.7352 1.0332 ⋯
b_τx -0.7958 0.1043 0.0215 25.2635 35.6350 1.0352 ⋯
b_e 0.0000 0.0000 NaN NaN NaN NaN ⋯
b_ω 0.0000 0.0000 NaN NaN NaN NaN ⋯
b_a 0.0094 0.0001 0.0000 70.6421 78.3393 1.0070 ⋯
b_τ -0.1445 0.0165 0.0017 86.6636 112.1358 1.0134 ⋯
b_tp 57781.9466 0.0061 0.0006 86.6636 112.1358 1.0134 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 0.7877 0.8059 0.8197 0.8341 0.8564
rv0_1 -6715.4600 -6701.4539 -6695.7258 -6689.8379 -6680.0674
rv0_2 -29.9679 -16.0618 -8.6449 -4.3451 6.4307
jitter_1 0.0143 0.0332 0.1321 0.5145 1.7714
jitter_2 2.2528 4.1947 4.8123 5.3132 6.7735
gp_η₁ 19.2658 21.8813 24.3038 27.9899 30.9067
gp_η₂ 10.9504 12.8017 13.6752 14.3220 16.1688
gp_η₃ 7.3791 8.0685 9.0398 9.3980 10.2704
gp_η₄ 0.3539 0.3715 0.3809 0.3920 0.4128
b_P 0.0010 0.0010 0.0010 0.0010 0.0010
b_mass 0.0163 0.0193 0.0213 0.0233 0.0278
b_τy 0.4529 0.5356 0.6205 0.6881 0.8416
b_τx -1.0003 -0.8569 -0.7822 -0.7341 -0.6070
b_e 0.0000 0.0000 0.0000 0.0000 0.0000
b_ω 0.0000 0.0000 0.0000 0.0000 0.0000
b_a 0.0093 0.0094 0.0094 0.0095 0.0096
b_τ -0.1754 -0.1528 -0.1468 -0.1336 -0.1124
b_tp 57781.9352 57781.9436 57781.9458 57781.9507 57781.9585
For real data, we would want to increase the adaptation and iterations to about 1000 each.
Plot the results:
fig = OctofitterRadialVelocity.rvpostplot(model, chain) # saved to "k2_132-rvpostplot.png"
Corner plot:
octocorner(model, chain, small=true) # saved to "k2_132-pairplot-small.png"
It is expensive (per iteration) to fit Guassian processes. If you have many cores (or a cluster) available you might also try starting julia with multiple threads (julia --threads=auto) and using octofit_pigeons. Pigeons is less efficient with single-core, but has better parallelization scaling:
using Pigeons
chain, pt = octofit_pigeons(
model,
n_rounds=9
)
display(chain)