GitHub

Fit Relative RV Data

Octofitter includes support for fitting relative radial velocity data. Currently this is only tested with a single companion. Please open an issue if you would like to fit multiple companions simultaneously.

The convention we adopt is that positive relative radial velocity is the velocity of the companion (exoplanets) minus the velocity of the host (star).

To fit relative RV data, start by creating a likelihood object:

using Octofitter
using OctofitterRadialVelocity
using CairoMakie
using Distributions

rel_rv_like = PlanetRelativeRVLikelihood(
        (;epoch=5000.0,  rv=−6.54, σ_rv=1.30),
        (;epoch=5050.1,  rv=−3.33, σ_rv=1.09),
        (;epoch=5100.2,  rv=7.90,  σ_rv=.11),
    )
PlanetRelativeRVLikelihood Table with 4 columns and 3 rows:
     inst_idx  epoch   rv     σ_rv
   ┌──────────────────────────────
 1 │ 1         5000.0  -6.54  1.3
 2 │ 1         5050.1  -3.33  1.09
 3 │ 1         5100.2  7.9    0.11

The columns in the table should be . See the standard radial velocity tutorial for examples on how this data can be loaded from a CSV file.

The relative RV likelihood does not incorporate an instrument-specific RV offset. A jitter parameter called jitter can still be specified in the @planet block, as can parameters for a gaussian process model of stellar noise. Unlike the StarAbsoluteRVLikelihood, only a single instrument jitter parameter is supported. If you need to model relative radial velocities from multiple instruments with different jitters, please open an issue on GitHub.

Next, create a planet and system model, attaching the relative rv likelihood to the planet. Make sure to add a jitter parameter (optionally jitter=0) to the planet.

@planet b Visual{KepOrbit} begin
    a ~ truncated(Normal(10, 4), lower=0, upper=100)
    e ~ Uniform(0.0, 0.5)
    i ~ Sine()
    ω ~ UniformCircular()
    Ω ~ UniformCircular()
    θ ~ UniformCircular()
    tp = θ_at_epoch_to_tperi(system,b,50420)

    jitter ~ LogUniform(1, 1000) # m/s
end rel_rv_like

@system ExampleSystem begin
    M ~ truncated(Normal(1.2, 0.1), lower=0)
    plx ~ truncated(Normal(50.0, 0.02), lower=0)
end b

model = Octofitter.LogDensityModel(ExampleSystem)
LogDensityModel for System ExampleSystem of dimension 12 with fields .ℓπcallback and .∇ℓπcallback

using Random
rng = Random.Xoshiro(1234)

chain = octofit(rng, model)
Chains MCMC chain (1000×29×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 46.58 seconds
Compute duration  = 46.58 seconds
parameters        = M, plx, b_a, b_e, b_i, b_ωy, b_ωx, b_Ωy, b_Ωx, b_θy, b_θx, b_jitter, b_ω, b_Ω, 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      rh ⋯
      Symbol      Float64     Float64    Float64    Float64    Float64   Float ⋯

           M       1.1825      0.1014     0.0042   596.6606   597.7892    0.99 ⋯
         plx      49.9992      0.0209     0.0007   811.1600   534.6795    0.99 ⋯
         b_a      14.6848      3.0433     0.5711    28.7613   128.0578    1.01 ⋯
         b_e       0.2934      0.1413     0.0053   673.1762   504.3042    1.01 ⋯
         b_i       3.0208      0.1692     0.0141   142.6690   133.4957    1.00 ⋯
        b_ωy      -0.0040      0.6925     0.0461   239.4562   438.7579    1.00 ⋯
        b_ωx       0.3412      0.6518     0.1866    16.5009    40.9496    1.05 ⋯
        b_Ωy      -0.0573      0.7022     0.0625   133.6191   389.0335    1.00 ⋯
        b_Ωx       0.0093      0.7178     0.0540   199.9862   602.6893    1.00 ⋯
        b_θy       0.0026      0.7245     0.0513   228.9825   650.1640    1.00 ⋯
        b_θx       0.1111      0.6952     0.0486   197.2977   675.2340    1.00 ⋯
    b_jitter      15.7645     33.6731     5.5142    15.3260    38.5221    1.08 ⋯
         b_ω       0.6414      1.6842     0.2906    55.6819   509.8109    1.03 ⋯
         b_Ω       0.0771      1.8739     0.1393   199.1644   382.5198    1.00 ⋯
         b_θ       0.2198      1.8125     0.1221   254.2465   169.2113    1.00 ⋯
        b_tp   41385.9575   7446.7662   663.5239   147.0845   180.9039    1.00 ⋯
                                                               2 columns omitted

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

           M       0.9971       1.1090       1.1865       1.2522       1.3800
         plx      49.9576      49.9848      49.9989      50.0140      50.0408
         b_a      10.1127      12.2416      14.2243      16.6739      21.1635
         b_e       0.0252       0.1750       0.3150       0.4141       0.4920
         b_i       2.5138       3.0209       3.0736       3.1021       3.1345
        b_ωy      -1.0609      -0.6791       0.0094       0.6411       1.0356
        b_ωx      -1.0071      -0.1054       0.5597       0.8759       1.1107
        b_Ωy      -1.0721      -0.7486      -0.0521       0.5997       1.0377
        b_Ωx      -1.0300      -0.7070       0.0372       0.6985       1.0859
        b_θy      -1.0590      -0.6975      -0.0702       0.7212       1.0851
        b_θx      -1.0768      -0.5300       0.2217       0.7434       1.0876
    b_jitter       1.1585       2.5140       4.7737      12.0969     129.1123
         b_ω      -2.9791      -0.2946       1.0112       1.9226       2.9249
         b_Ω      -3.0228      -1.4773       0.0912       1.7658       3.0142
         b_θ      -2.9845      -1.2890       0.4540       1.7899       3.0394
        b_tp   23685.9849   36214.0091   42971.9825   47972.9902   50297.9010

octoplot(model, chain)
Example block output