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.54 seconds
Compute duration = 2.54 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.2028 0.0971 0.0029 1096.4454 786.2712 1.0004 ⋯
b_a 1.2786 0.0376 0.0011 1101.2986 839.8138 1.0011 ⋯
b_e 0.1270 0.1013 0.0031 644.2779 323.2730 1.0109 ⋯
b_i 1.5535 0.4389 0.0159 829.7299 761.7017 0.9992 ⋯
b_ωy 0.0496 0.7014 0.0424 295.0209 504.0216 1.0000 ⋯
b_ωx -0.0043 0.7251 0.0387 395.9884 762.3949 0.9997 ⋯
b_Ωy -0.0121 0.7318 0.0411 337.7874 738.2036 1.0080 ⋯
b_Ωx 0.0522 0.6906 0.0385 350.4210 682.5667 1.0026 ⋯
b_τy -0.0392 0.6925 0.0408 319.2962 766.6262 0.9999 ⋯
b_τx 0.0032 0.7266 0.0386 402.5417 814.3606 0.9999 ⋯
b_gamma 103.4582 199.9965 7.2593 938.7737 714.3181 0.9997 ⋯
b_ω -0.0133 1.7416 0.0849 490.7034 636.2739 0.9994 ⋯
b_Ω 0.1407 1.8422 0.0972 425.9464 564.3586 1.0115 ⋯
b_τ -0.0133 0.2951 0.0134 582.9975 702.3693 0.9996 ⋯
b_P 1.3197 0.0236 0.0010 626.8262 789.8283 0.9996 ⋯
b_tp 5993.6657 142.3333 6.4460 597.9629 646.4712 0.9997 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 1.0195 1.1333 1.2013 1.2680 1.3945
b_a 1.2034 1.2551 1.2794 1.3037 1.3526
b_e 0.0039 0.0479 0.1020 0.1833 0.3776
b_i 0.8044 1.1991 1.5354 1.9288 2.3287
b_ωy -1.0566 -0.6182 0.0961 0.7268 1.0709
b_ωx -1.0727 -0.7604 0.0516 0.6956 1.0581
b_Ωy -1.0961 -0.7369 -0.0008 0.7291 1.0508
b_Ωx -1.0591 -0.6123 0.1298 0.7005 1.0670
b_τy -1.0614 -0.7197 -0.0463 0.6078 1.0650
b_τx -1.0492 -0.7123 -0.0318 0.7674 1.0978
b_gamma 0.1252 0.9471 11.1641 85.9736 778.2238
b_ω -2.9235 -1.5302 0.0766 1.4141 2.9271
b_Ω -3.0176 -1.4105 0.2520 1.7934 3.0086
b_τ -0.4786 -0.2804 -0.0070 0.2377 0.4763
b_P 1.2757 1.3039 1.3187 1.3356 1.3684
b_tp 5768.9990 5865.7238 5996.6471 6114.8212 6230.3110
octoplot(model, chain, ts=range(4500, 8000, length=200), show_physical_orbit=true)