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.31 seconds
Compute duration = 2.31 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.2006 0.1016 0.0034 902.5940 814.3606 0.9998 ⋯
b_a 1.2772 0.0387 0.0013 907.2830 814.6575 1.0007 ⋯
b_e 0.1355 0.1014 0.0032 874.6595 512.9898 0.9994 ⋯
b_i 1.5741 0.4447 0.0163 801.2321 912.4941 0.9990 ⋯
b_ωy 0.0168 0.7160 0.0426 303.9276 698.3219 0.9995 ⋯
b_ωx -0.0308 0.7114 0.0422 290.5184 661.9088 1.0020 ⋯
b_Ωy 0.0226 0.7138 0.0371 417.1016 815.9685 0.9992 ⋯
b_Ωx 0.0082 0.7193 0.0358 428.5418 577.7946 0.9993 ⋯
b_τy -0.0154 0.7177 0.0421 338.4280 787.3716 1.0007 ⋯
b_τx 0.0208 0.7121 0.0424 286.7535 651.2125 1.0024 ⋯
b_gamma 105.0832 203.8283 6.8552 1095.6684 754.8793 0.9997 ⋯
b_ω -0.0901 1.7954 0.0833 537.5279 627.6239 0.9992 ⋯
b_Ω 0.0035 1.7862 0.0826 532.2335 814.3606 0.9992 ⋯
b_τ 0.0192 0.2914 0.0141 541.7141 790.6797 1.0024 ⋯
b_P 1.3190 0.0226 0.0007 1016.6857 814.6575 1.0020 ⋯
b_tp 6009.2681 140.4231 6.7632 544.2547 788.5790 1.0024 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 1.0082 1.1306 1.1986 1.2662 1.4186
b_a 1.2005 1.2511 1.2770 1.3029 1.3568
b_e 0.0046 0.0529 0.1151 0.1969 0.3717
b_i 0.7916 1.2061 1.6096 1.9366 2.3412
b_ωy -1.0652 -0.6830 0.0373 0.7088 1.0733
b_ωx -1.0611 -0.7145 -0.0758 0.6662 1.0451
b_Ωy -1.0525 -0.6873 0.0218 0.7326 1.0699
b_Ωx -1.0591 -0.6985 -0.0165 0.7072 1.0793
b_τy -1.0767 -0.7217 -0.0178 0.6705 1.0707
b_τx -1.0760 -0.6888 0.0765 0.7167 1.0611
b_gamma 0.1278 1.0534 9.8736 92.6457 782.1337
b_ω -3.0334 -1.6239 -0.1884 1.4287 2.9477
b_Ω -3.0016 -1.5023 -0.0293 1.6277 2.9220
b_τ -0.4666 -0.2280 0.0294 0.2729 0.4810
b_P 1.2725 1.3043 1.3188 1.3340 1.3623
b_tp 5775.1018 5890.5365 6014.0611 6132.2245 6232.2501
octoplot(model, chain, ts=range(4500, 8000, length=200), show_physical_orbit=true)