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.

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

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

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

# 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
end
gauss_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     = 271.71 seconds
Compute duration  = 271.71 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.8232    0.0190    0.0017   124.5380    76.6284    0.992 ⋯
     rv0_harps   -6702.8758   11.4017    2.5353    20.5993    40.5549    1.091 ⋯
       rv0_pfs      -5.2308   16.3851    5.2647     9.6967    31.4535    1.054 ⋯
  jitter_harps       1.0700    1.1351    0.1847    48.1501    51.5563    0.998 ⋯
    jitter_pfs       5.0222    1.5428    0.1716    79.1528    37.7807    0.999 ⋯
         gp_η₁      27.3924    5.3189    1.6247    11.1771    74.7352    1.098 ⋯
         gp_η₂      13.3158    1.5311    0.1758    79.5433    59.8641    1.033 ⋯
         gp_η₃       9.2864    0.8471    0.1051    73.2295    74.9298    0.992 ⋯
         gp_η₄       0.3861    0.0161    0.0013   143.2572    90.1814    1.000 ⋯
           b_P       0.0010    0.0000    0.0000    94.8622    52.4503    0.994 ⋯
        b_mass       0.0212    0.0046    0.0006    64.0677    55.9475    1.006 ⋯
          b_τy       0.6519    0.1735    0.0180    90.7584    52.4503    0.993 ⋯
          b_τx      -0.7249    0.1554    0.0164    95.9558   116.7774    0.994 ⋯
           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   124.3909    76.6284    0.992 ⋯
           b_τ      -0.1337    0.0342    0.0038    79.1643    78.3393    0.994 ⋯
          b_tp   57781.9506    0.0126    0.0014    79.1643    78.3393    0.994 ⋯
                                                               2 columns omitted

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

             M       0.7870       0.8099       0.8215       0.8355       0.858 ⋯
     rv0_harps   -6722.8730   -6709.0452   -6703.7440   -6696.5815   -6675.521 ⋯
       rv0_pfs     -29.2955     -20.0840      -3.6088       6.7164      23.629 ⋯
  jitter_harps       0.1186       0.2681       0.5651       1.3611       3.767 ⋯
    jitter_pfs       1.8800       4.1409       4.9843       6.0042       8.537 ⋯
         gp_η₁      18.0461      24.1407      27.1589      32.0243      35.990 ⋯
         gp_η₂      10.4429      12.4486      13.1768      14.3107      16.251 ⋯
         gp_η₃       7.7574       8.8969       9.1670       9.6232      11.292 ⋯
         gp_η₄       0.3564       0.3758       0.3862       0.3966       0.420 ⋯
           b_P       0.0010       0.0010       0.0010       0.0010       0.001 ⋯
        b_mass       0.0135       0.0181       0.0214       0.0238       0.030 ⋯
          b_τy       0.3084       0.5240       0.6681       0.7873       0.957 ⋯
          b_τx      -0.9902      -0.8291      -0.7521      -0.6216      -0.393 ⋯
           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.1992      -0.1596      -0.1315      -0.1091      -0.060 ⋯
          b_tp   57781.9264   57781.9411   57781.9514   57781.9597   57781.977 ⋯
                                                                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"
Example block output

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

Corner plot:

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