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 = 3.23 seconds
Compute duration = 3.23 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.2054 0.0993 0.0032 930.9540 786.8536 1.0002 ⋯
b_a 1.2797 0.0389 0.0013 942.2525 720.1504 0.9996 ⋯
b_e 0.1278 0.0984 0.0033 719.0320 539.8307 1.0008 ⋯
b_i 1.5872 0.4484 0.0181 696.9473 795.5730 1.0053 ⋯
b_ωy 0.0584 0.7241 0.0407 336.4811 589.1057 1.0010 ⋯
b_ωx 0.0682 0.7001 0.0391 346.6255 704.3724 1.0030 ⋯
b_Ωy 0.0638 0.7269 0.0418 323.1707 649.3709 1.0074 ⋯
b_Ωx 0.0195 0.7001 0.0346 431.5938 774.1110 0.9992 ⋯
b_τy -0.0510 0.7142 0.0398 339.0060 678.9875 1.0007 ⋯
b_τx -0.0719 0.6998 0.0400 318.3490 664.4983 1.0043 ⋯
b_gamma 109.6636 204.7931 6.5825 851.8009 697.5131 1.0061 ⋯
b_ω 0.1453 1.7585 0.0886 447.8169 656.2286 1.0035 ⋯
b_Ω -0.0077 1.7481 0.0772 588.2423 695.6823 1.0076 ⋯
b_τ -0.0210 0.2972 0.0139 556.4650 767.7530 1.0014 ⋯
b_P 1.3201 0.0248 0.0009 689.0794 743.9013 0.9995 ⋯
b_tp 5990.0115 143.4618 6.7154 560.8385 745.7481 1.0013 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 1.0170 1.1399 1.2028 1.2705 1.4110
b_a 1.2032 1.2536 1.2792 1.3065 1.3570
b_e 0.0055 0.0492 0.1063 0.1902 0.3641
b_i 0.7930 1.2263 1.6147 1.9515 2.3415
b_ωy -1.0515 -0.6841 0.1485 0.7544 1.0945
b_ωx -1.0589 -0.5791 0.1200 0.7329 1.0663
b_Ωy -1.0697 -0.6594 0.1217 0.7738 1.0746
b_Ωx -1.0617 -0.6606 0.0173 0.7259 1.0629
b_τy -1.0725 -0.7500 -0.1402 0.6631 1.0345
b_τx -1.0692 -0.7466 -0.1428 0.5826 1.0541
b_gamma 0.1210 1.1153 10.9876 100.1666 812.2740
b_ω -2.9957 -1.3158 0.2504 1.5598 3.0213
b_Ω -2.9922 -1.4667 0.0193 1.4279 2.9221
b_τ -0.4738 -0.2937 -0.0431 0.2486 0.4776
b_P 1.2731 1.3033 1.3190 1.3363 1.3698
b_tp 5770.7702 5858.5148 5979.3610 6120.6100 6232.7554
octoplot(model, chain, ts=range(4500, 8000, length=200), show_physical_orbit=true)