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, rv=-24022.74287804528, σ_rv=15000.0),
(epoch=5100, rv=-18571.333891168735, σ_rv=15000.0),
(epoch=5200, rv= 14221.562142944855, σ_rv=15000.0),
(epoch=5300, rv= 26076.885281031347, σ_rv=15000.0),
(epoch=5400, rv= -459.2622916989299, σ_rv=15000.0),
(epoch=5500, rv=-26319.264894263993, σ_rv=15000.0),
(epoch=5600, rv=-13430.95547916007, σ_rv=15000.0),
(epoch=5700, rv= 19230.962951723584, σ_rv=15000.0),
(epoch=5800, rv= 23580.261108170227, σ_rv=15000.0),
(epoch=5900, rv= -6786.277919597756, σ_rv=15000.0),
(epoch=6000, rv=-27161.777481651112, σ_rv=15000.0),
(epoch=6100, rv= -7548.583094927461, σ_rv=15000.0),
(epoch=6200, rv= 23177.948014103342, σ_rv=15000.0),
(epoch=6300, rv= 19780.94394128632, σ_rv=15000.0),
(epoch=6400, rv=-12738.38520102873, σ_rv=15000.0),
(epoch=6500, rv=-26503.73597982596, σ_rv=15000.0),
(epoch=6600, rv= -1249.188767321913, σ_rv=15000.0),
(epoch=6700, rv= 25844.465894406647, σ_rv=15000.0),
(epoch=6800, rv= 14888.827293969505, σ_rv=15000.0),
(epoch=6900, rv=-17986.75959839915, σ_rv=15000.0),
(epoch=7000, rv=-24381.49393255423, σ_rv=15000.0),
(epoch=7100, rv= 5119.21707156116, σ_rv=15000.0),
(epoch=7200, rv= 27083.2046462065, σ_rv=15000.0),
(epoch=7300, rv= 9174.176455190982, σ_rv=15000.0),
(epoch=7400, rv=-22241.45434114139, σ_rv=15000.0),
instrument_name="simulated data",
jitter=:gamma # name of jitter variable to use
)PlanetRelativeRVLikelihood Table with 3 columns and 25 rows:
epoch rv σ_rv
┌─────────────────────────
1 │ 5000 -24022.7 15000.0
2 │ 5100 -18571.3 15000.0
3 │ 5200 14221.6 15000.0
4 │ 5300 26076.9 15000.0
5 │ 5400 -459.262 15000.0
6 │ 5500 -26319.3 15000.0
7 │ 5600 -13431.0 15000.0
8 │ 5700 19231.0 15000.0
9 │ 5800 23580.3 15000.0
10 │ 5900 -6786.28 15000.0
11 │ 6000 -27161.8 15000.0
12 │ 6100 -7548.58 15000.0
13 │ 6200 23177.9 15000.0
14 │ 6300 19780.9 15000.0
15 │ 6400 -12738.4 15000.0
16 │ 6500 -26503.7 15000.0
17 │ 6600 -1249.19 15000.0
⋮ │ ⋮ ⋮ ⋮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 can still be specified in the @planet block, as can parameters for a gaussian process model of stellar noise. Currently 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 KepOrbit begin
a ~ Uniform(0,10)
e ~ Uniform(0.0, 0.5)
i ~ Sine()
ω ~ UniformCircular()
Ω ~ UniformCircular()
τ ~ UniformCircular(1.0)
P = √(b.a^3/system.M)
tp = b.τ*b.P*365.25 + 6000 # reference epoch for τ. Choose an MJD date near your data.
gamma ~ LogUniform(0.1, 1000) # m/s
end rel_rv_like
@system ExampleSystem begin
M ~ truncated(Normal(1.2, 0.1), lower=0.1)
end b
model = Octofitter.LogDensityModel(ExampleSystem)LogDensityModel for System ExampleSystem of dimension 11 and 28 epochs with fields .ℓπcallback and .∇ℓπcallback
using Random
rng = Random.Xoshiro(123)
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 = 2.93 seconds
Compute duration = 2.93 seconds
parameters = M, b_a, b_e, b_i, b_ωy, b_ωx, b_Ωy, b_Ωx, b_τy, b_τx, b_gamma, b_ω, b_Ω, b_τ, b_P, 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 1.1994 0.1007 0.0035 830.9348 584.8309 1.0045 ⋯
b_a 1.2784 0.0389 0.0013 885.0113 592.2101 0.9995 ⋯
b_e 0.1338 0.1015 0.0034 795.2039 636.1553 0.9994 ⋯
b_i 1.6051 0.4340 0.0180 637.3436 724.6882 1.0006 ⋯
b_ωy 0.0790 0.7125 0.0395 351.1947 576.0815 1.0071 ⋯
b_ωx 0.0180 0.7069 0.0391 353.7675 780.5245 1.0149 ⋯
b_Ωy 0.0269 0.7073 0.0363 438.2557 941.3886 1.0013 ⋯
b_Ωx 0.0381 0.7179 0.0395 339.7477 697.0767 0.9997 ⋯
b_τy -0.0852 0.7133 0.0395 337.2047 587.8279 1.0072 ⋯
b_τx -0.0069 0.7082 0.0385 376.2813 699.6500 1.0113 ⋯
b_gamma 112.1491 208.1277 7.4822 731.2408 663.0044 0.9997 ⋯
b_ω 0.0646 1.7258 0.0855 467.3656 836.8318 1.0162 ⋯
b_Ω 0.0391 1.7720 0.0793 586.5977 732.3138 0.9993 ⋯
b_τ -0.0100 0.3036 0.0138 612.8632 639.6115 1.0070 ⋯
b_P 1.3215 0.0239 0.0009 726.1005 765.2791 0.9996 ⋯
b_tp 5994.9986 146.7359 6.6761 616.0199 596.0926 1.0067 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 1.0066 1.1312 1.1984 1.2687 1.3995
b_a 1.2027 1.2517 1.2793 1.3047 1.3528
b_e 0.0056 0.0527 0.1143 0.1889 0.3786
b_i 0.8401 1.2624 1.6211 1.9570 2.3869
b_ωy -1.0610 -0.6131 0.1622 0.7737 1.0680
b_ωx -1.0494 -0.6815 0.0525 0.6797 1.0915
b_Ωy -1.0411 -0.6586 0.0346 0.7198 1.0650
b_Ωx -1.0689 -0.6522 0.0688 0.7365 1.0666
b_τy -1.0843 -0.7607 -0.1769 0.6138 1.0593
b_τx -1.0536 -0.6769 -0.0818 0.6540 1.0764
b_gamma 0.1166 0.9492 10.1992 116.1508 808.4711
b_ω -2.9951 -1.3516 0.0910 1.4710 2.9661
b_Ω -2.9393 -1.5239 0.1159 1.5415 2.9174
b_τ -0.4802 -0.2870 -0.0239 0.2673 0.4745
b_P 1.2762 1.3054 1.3209 1.3375 1.3682
b_tp 5767.9186 5863.0236 5988.7513 6128.7993 6228.5936
octoplot(model, chain, ts=range(4500, 8000, length=200), show_physical_orbit=true)