Fit Radial Velocity

You can use Octofitter as a basic tool for fitting radial velocity data, by itself, or in combination with other kinds of data. Multiple instruments (up to five) are supported.

Tip

Radial velocity modelling is supported in Octofitter via the extension package OctofitterRadialVelocity. To install it, run pkg> add OctofitterRadialVelocity

Load the packages we'll need:

using Octofitter, OctofitterRadialVelocity, Distributions, PlanetOrbits, Plots

We can specify a table of radial velocity data manually by creating a RadialVelocityLikelihood. An example of this is in the other RV tutorial.

We can also directly load in data from the HARPS RVBank dataset:

rvs = OctofitterRadialVelocity.HARPS_rvs("GJ436")
RadialVelocityLikelihood Table with 4 columns and 169 rows:
      epoch    inst_idx  rv       σ_rv
    ┌──────────────────────────────────
 1  │ 53760.3  1         14.41    1.193
 2  │ 53761.3  1         -18.072  1.359
 3  │ 53762.3  1         7.407    0.96
 4  │ 53763.3  1         9.38     1.146
 5  │ 53765.3  1         13.463   1.15
 6  │ 53785.3  1         -20.19   1.175
 7  │ 53788.3  1         -11.405  1.136
 8  │ 54122.3  1         13.836   1.152
 9  │ 54135.3  1         8.246    1.163
 10 │ 54140.3  1         0.079    1.259
 11 │ 54142.3  1         -19.134  1.034
 12 │ 54166.3  1         -17.128  1.388
 13 │ 54172.2  1         6.696    1.24
 14 │ 54194.2  1         13.143   1.488
 15 │ 54197.2  1         -1.356   1.267
 16 │ 54199.2  1         16.157   1.029
 17 │ 54202.2  1         13.661   1.179
 ⋮  │    ⋮        ⋮         ⋮       ⋮

Now, create a planet. Since we're only fitting radial velocity data, we fix some of these parameters

@planet b Visual{KepOrbit} begin
    τ ~ UniformCircular(1.0)
    mass ~ truncated(Normal(21.3*0.00314558, 8*0.00314558),lower=0)
    a ~ truncated(Normal(0.028, 0.02), lower=0)
    i = pi/2
    e ~ Uniform(0, 0.5)
    ω ~ UniformCircular()
    Ω = 0
end

# Nor create the system
@system HD82134 begin
    M ~ truncated(Normal(0.425,0.009),lower=0)
    plx ~ 0.425

    # RV zero point of the system for instrument 1
    rv0_1 ~ Normal(0,10)
    # Jitter term for instrument 1
    jitter_1 ~ truncated(Normal(0,10),lower=0)
end rvs b
System model HD82134
Derived:
  plx, 
Priors:
  M	Truncated(Distributions.Normal{Float64}(μ=0.425, σ=0.009); lower=0.0)
  rv0_1	Distributions.Normal{Float64}(μ=0.0, σ=10.0)
  jitter_1	Truncated(Distributions.Normal{Float64}(μ=0.0, σ=10.0); lower=0.0)
Planet b
Derived:
  τ, i, ω, Ω, 
Priors:
  τy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  τx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  mass	Truncated(Distributions.Normal{Float64}(μ=0.067000854, σ=0.02516464); lower=0.0)
  a	Truncated(Distributions.Normal{Float64}(μ=0.028, σ=0.02); lower=0.0)
  e	Distributions.Uniform{Float64}(a=0.0, b=0.5)
  ωy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  ωx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
Octofitter.UnitLengthPrior{:τy, :τx}: √(τy^2+τx^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:ωy, :ωx}: √(ωy^2+ωx^2) ~ LogNormal(log(1), 0.02)


RadialVelocityLikelihood Table with 4 columns and 169 rows:
      epoch    inst_idx  rv       σ_rv
    ┌──────────────────────────────────
 1  │ 53760.3  1         14.41    1.193
 2  │ 53761.3  1         -18.072  1.359
 3  │ 53762.3  1         7.407    0.96
 4  │ 53763.3  1         9.38     1.146
 5  │ 53765.3  1         13.463   1.15
 6  │ 53785.3  1         -20.19   1.175
 7  │ 53788.3  1         -11.405  1.136
 8  │ 54122.3  1         13.836   1.152
 9  │ 54135.3  1         8.246    1.163
 10 │ 54140.3  1         0.079    1.259
 11 │ 54142.3  1         -19.134  1.034
 12 │ 54166.3  1         -17.128  1.388
 13 │ 54172.2  1         6.696    1.24
 14 │ 54194.2  1         13.143   1.488
 15 │ 54197.2  1         -1.356   1.267
 16 │ 54199.2  1         16.157   1.029
 17 │ 54202.2  1         13.661   1.179
 ⋮  │    ⋮        ⋮         ⋮       ⋮

Build model:

model = Octofitter.LogDensityModel(HD82134; autodiff=:ForwardDiff, verbosity=4) # defaults are ForwardDiff, and verbosity=0
LogDensityModel for System HD82134 of dimension 10 with fields .ℓπcallback and .∇ℓπcallback

Sample from chains:

results = octofit(
    model, 0.8;
    adaptation =  1000,
    iterations =  1000,
    verbosity = 4,
    max_depth = 14
)
Chains MCMC chain (1000×28×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 7091.0 seconds
Compute duration  = 7091.0 seconds
parameters        = M, rv0_1, jitter_1, plx, b_τy, b_τx, b_mass, b_a, b_e, b_ωy, b_ωx, b_τ, b_i, b_ω, b_Ω
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.4253    0.0090    0.0003    714.2883   608.3330    1.0017    ⋯
       rv0_1    0.1196    0.2160    0.0084    660.5680   403.3236    1.0098    ⋯
    jitter_1    2.0699    0.1634    0.0064    653.2460   678.0667    1.0015    ⋯
         plx    0.4250    0.0000    0.0000         NaN        NaN       NaN    ⋯
        b_τy   -0.9539    0.0592    0.0026    639.3066   604.8516    1.0026    ⋯
        b_τx   -0.2133    0.2048    0.0096    465.4708   405.1876    1.0039    ⋯
      b_mass    0.0658    0.0014    0.0000    806.6865   715.0619    1.0040    ⋯
         b_a    0.0281    0.0002    0.0000    714.6850   608.3330    1.0017    ⋯
         b_e    0.1319    0.0173    0.0005   1229.4010   774.0196    0.9999    ⋯
        b_ωy    0.7134    0.0952    0.0034    798.7886   607.7104    0.9996    ⋯
        b_ωx   -0.6877    0.0981    0.0035    807.4639   689.7141    1.0005    ⋯
         b_τ   -0.3170    0.3340    0.0153    675.5063   739.5350    1.0009    ⋯
         b_i    1.5708    0.0000    0.0000         NaN        NaN       NaN    ⋯
         b_ω   -0.7671    0.1357    0.0048    793.8533   592.5683    1.0015    ⋯
         b_Ω    0.0000    0.0000       NaN         NaN        NaN       NaN    ⋯
                                                                1 column omitted

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

           M    0.4078    0.4192    0.4255    0.4315    0.4424
       rv0_1   -0.2971   -0.0249    0.1222    0.2696    0.5437
    jitter_1    1.7747    1.9581    2.0666    2.1738    2.3912
         plx    0.4250    0.4250    0.4250    0.4250    0.4250
        b_τy   -1.0277   -0.9930   -0.9703   -0.9297   -0.7996
        b_τx   -0.6036   -0.3595   -0.2076   -0.0806    0.2211
      b_mass    0.0633    0.0648    0.0659    0.0668    0.0685
         b_a    0.0277    0.0280    0.0281    0.0283    0.0285
         b_e    0.0975    0.1202    0.1318    0.1439    0.1646
        b_ωy    0.5031    0.6541    0.7216    0.7790    0.8824
        b_ωx   -0.8670   -0.7561   -0.6941   -0.6223   -0.4746
         b_τ   -0.4952   -0.4741   -0.4525   -0.4215    0.4959
         b_i    1.5708    1.5708    1.5708    1.5708    1.5708
         b_ω   -1.0388   -0.8561   -0.7664   -0.6724   -0.5018
         b_Ω    0.0000    0.0000    0.0000    0.0000    0.0000
using Plots
octoplot(model, results)
("HD82134-plot-grid.png", "HD82134-plot-grid-colorbar.png")