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..

Note

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 Distributions

We 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:
        epoch    rv        σ_rv   inst_idx
      ┌───────────────────────────────────
 1,1  │ 57782.2  -6682.78  3.71   1
 2,1  │ 57783.1  -6710.55  5.95   1
 3,1  │ 57783.2  -6698.96  8.76   1
 4,1  │ 57812.1  -6672.32  4.0    1
 5,1  │ 57812.2  -6672.0   4.22   1
 6,1  │ 57812.2  -6690.31  4.63   1
 7,1  │ 57813.0  -6701.27  4.12   1
 8,1  │ 57813.1  -6700.71  5.34   1
 9,1  │ 57813.1  -6695.97  3.67   1
 10,1 │ 57813.1  -6705.87  4.23   1
 11,1 │ 57813.1  -6694.27  5.12   1
 12,1 │ 57813.2  -6699.99  5.0    1
 13,1 │ 57813.2  -6696.17  10.43  1
 14,1 │ 57836.0  -6675.4   7.03   1
 15,1 │ 57836.0  -6665.92  7.48   1
 16,1 │ 57836.0  -6661.9   6.03   1
 17,1 │ 57836.1  -6657.92  5.08   1
  ⋮   │    ⋮        ⋮        ⋮       ⋮

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).

    # HARPS-N
    rv0_1 ~ Normal(0,10000) # m/s
    jitter_1 ~ LogUniform(0.01,100) # m/s

    # FPS
    rv0_2 ~ Normal(0,10000) # m/s
    jitter_2 ~ LogUniform(0.01,100) # m/s
end rvlike b
System model k2_132
Derived:
  
Priors:
  M	Truncated(Distributions.Normal{Float64}(μ=0.82, σ=0.02); lower=0.0)
  rv0_1	Distributions.Normal{Float64}(μ=0.0, σ=10000.0)
  jitter_1	Distributions.LogUniform{Float64}(a=0.01, b=100.0)
  rv0_2	Distributions.Normal{Float64}(μ=0.0, σ=10000.0)
  jitter_2	Distributions.LogUniform{Float64}(a=0.01, 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:
        epoch    rv        σ_rv   inst_idx
      ┌───────────────────────────────────
 1,1  │ 57782.2  -6682.78  3.71   1
 2,1  │ 57783.1  -6710.55  5.95   1
 3,1  │ 57783.2  -6698.96  8.76   1
 4,1  │ 57812.1  -6672.32  4.0    1
 5,1  │ 57812.2  -6672.0   4.22   1
 6,1  │ 57812.2  -6690.31  4.63   1
 7,1  │ 57813.0  -6701.27  4.12   1
 8,1  │ 57813.1  -6700.71  5.34   1
 9,1  │ 57813.1  -6695.97  3.67   1
 10,1 │ 57813.1  -6705.87  4.23   1
 11,1 │ 57813.1  -6694.27  5.12   1
 12,1 │ 57813.2  -6699.99  5.0    1
 13,1 │ 57813.2  -6696.17  10.43  1
 14,1 │ 57836.0  -6675.4   7.03   1
 15,1 │ 57836.0  -6665.92  7.48   1
 16,1 │ 57836.0  -6661.9   6.03   1
 17,1 │ 57836.1  -6657.92  5.08   1
  ⋮   │    ⋮        ⋮        ⋮       ⋮

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     = 26.3 seconds
Compute duration  = 26.3 seconds
parameters        = M, rv0_1, jitter_1, rv0_2, 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.8201    0.0193    0.0005   1274.3451   500.3212    1.0000 ⋯
       rv0_1    6693.6090    2.2242    0.0562   1534.9792   852.6950    1.0006 ⋯
    jitter_1      13.7415    1.9722    0.0596   1192.0104   732.3881    1.0003 ⋯
       rv0_2      -1.5417    3.3008    0.1225    738.1237   550.6909    1.0010 ⋯
    jitter_2      17.7521    2.7468    0.0823   1188.7967   647.0574    0.9994 ⋯
         b_P       0.0010    0.0000    0.0000   1088.1344   786.9206    0.9993 ⋯
        b_τy       0.7586    0.1594    0.0052   1012.0497   623.2499    1.0045 ⋯
        b_τx      -0.6143    0.1913    0.0068    769.9153   752.7169    1.0018 ⋯
      b_mass       0.0419    0.0088    0.0003   1024.7643   583.5152    0.9999 ⋯
         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   1274.5119   500.3212    1.0000 ⋯
         b_τ      -0.1081    0.0370    0.0013    819.7533   717.2317    1.0026 ⋯
        b_tp   57781.9601    0.0137    0.0005    819.7857   717.2317    1.0025 ⋯
                                                                1 column omitted

Quantiles
  parameters         2.5%        25.0%        50.0%        75.0%        97.5%
      Symbol      Float64      Float64      Float64      Float64      Float64

           M       0.7828       0.8074       0.8211       0.8334       0.8575
       rv0_1    6689.4474    6692.0466    6693.5929    6695.1320    6698.0439
    jitter_1      10.5928      12.2987      13.5502      14.9908      17.9181
       rv0_2      -8.1279      -3.6562      -1.3882       0.7486       4.8593
    jitter_2      13.1972      15.8133      17.4410      19.3698      23.6145
         b_P       0.0010       0.0010       0.0010       0.0010       0.0010
        b_τy       0.4055       0.6661       0.7655       0.8628       1.0557
        b_τx      -0.9728      -0.7473      -0.6247      -0.4918      -0.2062
      b_mass       0.0250       0.0365       0.0417       0.0476       0.0600
         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.1819      -0.1325      -0.1078      -0.0847      -0.0326
        b_tp   57781.9328   57781.9511   57781.9602   57781.9687   57781.9880

Excellent! Let's plot the maximum likelihood orbit:

using CairoMakie: Makie
fig = OctofitterRadialVelocity.rvpostplot(model, chain) # saved to "k2_132-rvpostplot.png"
Example block output

Create a corner plot:

using PairPlots
using CairoMakie: Makie
octocorner(model, chain, small=true)
Example block output

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
end
gauss_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:
        epoch    rv        σ_rv   inst_idx
      ┌───────────────────────────────────
 1,1  │ 57782.2  -6682.78  3.71   1
 2,1  │ 57783.1  -6710.55  5.95   1
 3,1  │ 57783.2  -6698.96  8.76   1
 4,1  │ 57812.1  -6672.32  4.0    1
 5,1  │ 57812.2  -6672.0   4.22   1
 6,1  │ 57812.2  -6690.31  4.63   1
 7,1  │ 57813.0  -6701.27  4.12   1
 8,1  │ 57813.1  -6700.71  5.34   1
 9,1  │ 57813.1  -6695.97  3.67   1
 10,1 │ 57813.1  -6705.87  4.23   1
 11,1 │ 57813.1  -6694.27  5.12   1
 12,1 │ 57813.2  -6699.99  5.0    1
 13,1 │ 57813.2  -6696.17  10.43  1
 14,1 │ 57836.0  -6675.4   7.03   1
 15,1 │ 57836.0  -6665.92  7.48   1
 16,1 │ 57836.0  -6661.9   6.03   1
 17,1 │ 57836.1  -6657.92  5.08   1
  ⋮   │    ⋮        ⋮        ⋮       ⋮

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_1 ~ Normal(0,10000)
    jitter_1 ~ LogUniform(0.01,10)
    # FPS
    rv0_2 ~ Normal(0,10000)
    jitter_2 ~ LogUniform(0.01,10)

    # Add priors on GP kernel hyper-parameters.
    gp_η₁ ~ truncated(Normal(25,10),lower=0)
    gp_η₂ ~ truncated(Normal(gp_explength_mean,gp_explength_unc),lower=0)
    gp_η₃ ~ truncated(Normal(gp_per_mean,1),lower=0)
    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     = 174.24 seconds
Compute duration  = 174.24 seconds
parameters        = M, rv0_1, jitter_1, rv0_2, 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.8228    0.0179    0.0019    87.6161   112.1358    1.0184  ⋯
       rv0_1    6694.3794    8.8497    1.0865    73.8616    54.8565    1.0118  ⋯
    jitter_1       0.5429    0.6571    0.1146    28.7168    40.3486    1.0051  ⋯
       rv0_2      11.7672    7.8891    2.2005    14.5372    18.5440    1.1540  ⋯
    jitter_2       4.8548    0.8809    0.0698   152.5735   111.3156    0.9927  ⋯
       gp_η₁      24.8050    4.3166    0.5418    66.4700    52.4503    0.9969  ⋯
       gp_η₂      13.3590    1.2691    0.1347    89.7389    92.5822    0.9943  ⋯
       gp_η₃       8.9361    0.9710    0.1398    49.0221    34.3628    1.0039  ⋯
       gp_η₄       0.3825    0.0175    0.0020    74.8528    85.3949    0.9931  ⋯
         b_P       0.0010    0.0000    0.0000    38.8740    78.3393    1.0479  ⋯
      b_mass       0.0211    0.0037    0.0004    91.6137    31.7840    0.9989  ⋯
        b_τy       0.6571    0.1108    0.0094   145.1000   108.1275    0.9902  ⋯
        b_τx      -0.7384    0.1100    0.0115    87.7388    76.8221    1.0060  ⋯
         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    87.6420   112.1358    1.0194  ⋯
         b_τ      -0.1342    0.0211    0.0020   100.8556    78.3393    1.0023  ⋯
        b_tp   57781.9504    0.0078    0.0007   100.8556    78.3393    1.0024  ⋯
                                                                1 column omitted

Quantiles
  parameters         2.5%        25.0%        50.0%        75.0%        97.5%
      Symbol      Float64      Float64      Float64      Float64      Float64

           M       0.7933       0.8094       0.8225       0.8353       0.8631
       rv0_1    6680.6710    6688.5026    6693.6245    6699.2126    6712.5899
    jitter_1       0.0105       0.0491       0.2612       0.8191       2.3502
       rv0_2       0.8882       5.4391      10.5244      17.0441      27.6617
    jitter_2       3.4228       4.2927       4.9417       5.3917       6.4695
       gp_η₁      18.4239      22.0537      24.1592      26.3464      34.9753
       gp_η₂      11.0168      12.4598      13.4879      14.3559      15.5855
       gp_η₃       7.0705       8.3382       9.0853       9.3825      11.3554
       gp_η₄       0.3551       0.3711       0.3818       0.3931       0.4150
         b_P       0.0010       0.0010       0.0010       0.0010       0.0010
      b_mass       0.0141       0.0187       0.0215       0.0233       0.0287
        b_τy       0.4507       0.5769       0.6568       0.7386       0.8578
        b_τx      -0.9482      -0.8095      -0.7392      -0.6666      -0.5651
         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.1724      -0.1486      -0.1347      -0.1171      -0.1007
        b_tp   57781.9363   57781.9451   57781.9502   57781.9568   57781.9628

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"
Example block output

Corner plot:

octocorner(model, chain, small=true) # saved to "k2_132-pairplot-small.png"
Example block output

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=10
)
display(chain)