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](cd8462b0.png)