Fit Radial Velocity
You can use Octofitter as a basic tool for fitting radial velocity data, by itself, or in combination with other kinds of data. Multiple instruments (up to five) are supported.
Radial velocity modelling is supported in Octofitter via the extension package OctofitterRadialVelocity. To install it, run pkg> add OctofitterRadialVelocity
Load the packages we'll need:
using Octofitter, OctofitterRadialVelocity, Distributions, PlanetOrbits, PlotsWe can specify a table of radial velocity data manually by creating a RadialVelocityLikelihood. An example of this is in the other RV tutorial.
We can also directly load in data from the HARPS RVBank dataset:
rvs = OctofitterRadialVelocity.HARPS_rvs("GJ436")RadialVelocityLikelihood Table with 4 columns and 169 rows:
epoch inst_idx rv σ_rv
┌──────────────────────────────────
1 │ 53760.3 1 14.41 1.193
2 │ 53761.3 1 -18.072 1.359
3 │ 53762.3 1 7.407 0.96
4 │ 53763.3 1 9.38 1.146
5 │ 53765.3 1 13.463 1.15
6 │ 53785.3 1 -20.19 1.175
7 │ 53788.3 1 -11.405 1.136
8 │ 54122.3 1 13.836 1.152
9 │ 54135.3 1 8.246 1.163
10 │ 54140.3 1 0.079 1.259
11 │ 54142.3 1 -19.134 1.034
12 │ 54166.3 1 -17.128 1.388
13 │ 54172.2 1 6.696 1.24
14 │ 54194.2 1 13.143 1.488
15 │ 54197.2 1 -1.356 1.267
16 │ 54199.2 1 16.157 1.029
17 │ 54202.2 1 13.661 1.179
⋮ │ ⋮ ⋮ ⋮ ⋮Now, create a planet. Since we're only fitting radial velocity data, we fix some of these parameters
@planet b Visual{KepOrbit} begin
τ ~ UniformCircular(1.0)
mass ~ truncated(Normal(21.3*0.00314558, 8*0.00314558),lower=0)
a ~ truncated(Normal(0.028, 0.02), lower=0)
i = pi/2
e ~ Uniform(0, 0.5)
ω ~ UniformCircular()
Ω = 0
end
# Nor create the system
@system HD82134 begin
M ~ truncated(Normal(0.425,0.009),lower=0)
plx ~ 0.425
# RV zero point of the system for instrument 1
rv0_1 ~ Normal(0,10)
# Jitter term for instrument 1
jitter_1 ~ truncated(Normal(0,10),lower=0)
end rvs bSystem model HD82134
Derived:
plx,
Priors:
M Truncated(Distributions.Normal{Float64}(μ=0.425, σ=0.009); lower=0.0)
rv0_1 Distributions.Normal{Float64}(μ=0.0, σ=10.0)
jitter_1 Truncated(Distributions.Normal{Float64}(μ=0.0, σ=10.0); lower=0.0)
Planet b
Derived:
τ, i, ω, Ω,
Priors:
τy Distributions.Normal{Float64}(μ=0.0, σ=1.0)
τx Distributions.Normal{Float64}(μ=0.0, σ=1.0)
mass Truncated(Distributions.Normal{Float64}(μ=0.067000854, σ=0.02516464); lower=0.0)
a Truncated(Distributions.Normal{Float64}(μ=0.028, σ=0.02); lower=0.0)
e Distributions.Uniform{Float64}(a=0.0, b=0.5)
ωy Distributions.Normal{Float64}(μ=0.0, σ=1.0)
ωx Distributions.Normal{Float64}(μ=0.0, σ=1.0)
Octofitter.UnitLengthPrior{:τy, :τx}: √(τy^2+τx^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:ωy, :ωx}: √(ωy^2+ωx^2) ~ LogNormal(log(1), 0.02)
RadialVelocityLikelihood Table with 4 columns and 169 rows:
epoch inst_idx rv σ_rv
┌──────────────────────────────────
1 │ 53760.3 1 14.41 1.193
2 │ 53761.3 1 -18.072 1.359
3 │ 53762.3 1 7.407 0.96
4 │ 53763.3 1 9.38 1.146
5 │ 53765.3 1 13.463 1.15
6 │ 53785.3 1 -20.19 1.175
7 │ 53788.3 1 -11.405 1.136
8 │ 54122.3 1 13.836 1.152
9 │ 54135.3 1 8.246 1.163
10 │ 54140.3 1 0.079 1.259
11 │ 54142.3 1 -19.134 1.034
12 │ 54166.3 1 -17.128 1.388
13 │ 54172.2 1 6.696 1.24
14 │ 54194.2 1 13.143 1.488
15 │ 54197.2 1 -1.356 1.267
16 │ 54199.2 1 16.157 1.029
17 │ 54202.2 1 13.661 1.179
⋮ │ ⋮ ⋮ ⋮ ⋮
Build model:
model = Octofitter.LogDensityModel(HD82134; autodiff=:ForwardDiff, verbosity=4) # defaults are ForwardDiff, and verbosity=0LogDensityModel for System HD82134 of dimension 10 with fields .ℓπcallback and .∇ℓπcallback
Sample from chains:
results = octofit(
model, 0.8;
adaptation = 1000,
iterations = 1000,
verbosity = 4,
max_depth = 14
)Chains MCMC chain (1000×28×1 Array{Float64, 3}):
Iterations = 1:1:1000
Number of chains = 1
Samples per chain = 1000
Wall duration = 7091.0 seconds
Compute duration = 7091.0 seconds
parameters = M, rv0_1, jitter_1, plx, b_τy, b_τx, b_mass, b_a, b_e, b_ωy, b_ωx, b_τ, b_i, b_ω, b_Ω
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 0.4253 0.0090 0.0003 714.2883 608.3330 1.0017 ⋯
rv0_1 0.1196 0.2160 0.0084 660.5680 403.3236 1.0098 ⋯
jitter_1 2.0699 0.1634 0.0064 653.2460 678.0667 1.0015 ⋯
plx 0.4250 0.0000 0.0000 NaN NaN NaN ⋯
b_τy -0.9539 0.0592 0.0026 639.3066 604.8516 1.0026 ⋯
b_τx -0.2133 0.2048 0.0096 465.4708 405.1876 1.0039 ⋯
b_mass 0.0658 0.0014 0.0000 806.6865 715.0619 1.0040 ⋯
b_a 0.0281 0.0002 0.0000 714.6850 608.3330 1.0017 ⋯
b_e 0.1319 0.0173 0.0005 1229.4010 774.0196 0.9999 ⋯
b_ωy 0.7134 0.0952 0.0034 798.7886 607.7104 0.9996 ⋯
b_ωx -0.6877 0.0981 0.0035 807.4639 689.7141 1.0005 ⋯
b_τ -0.3170 0.3340 0.0153 675.5063 739.5350 1.0009 ⋯
b_i 1.5708 0.0000 0.0000 NaN NaN NaN ⋯
b_ω -0.7671 0.1357 0.0048 793.8533 592.5683 1.0015 ⋯
b_Ω 0.0000 0.0000 NaN NaN NaN NaN ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 0.4078 0.4192 0.4255 0.4315 0.4424
rv0_1 -0.2971 -0.0249 0.1222 0.2696 0.5437
jitter_1 1.7747 1.9581 2.0666 2.1738 2.3912
plx 0.4250 0.4250 0.4250 0.4250 0.4250
b_τy -1.0277 -0.9930 -0.9703 -0.9297 -0.7996
b_τx -0.6036 -0.3595 -0.2076 -0.0806 0.2211
b_mass 0.0633 0.0648 0.0659 0.0668 0.0685
b_a 0.0277 0.0280 0.0281 0.0283 0.0285
b_e 0.0975 0.1202 0.1318 0.1439 0.1646
b_ωy 0.5031 0.6541 0.7216 0.7790 0.8824
b_ωx -0.8670 -0.7561 -0.6941 -0.6223 -0.4746
b_τ -0.4952 -0.4741 -0.4525 -0.4215 0.4959
b_i 1.5708 1.5708 1.5708 1.5708 1.5708
b_ω -1.0388 -0.8561 -0.7664 -0.6724 -0.5018
b_Ω 0.0000 0.0000 0.0000 0.0000 0.0000
using Plots
octoplot(model, results)("HD82134-plot-grid.png", "HD82134-plot-grid-colorbar.png")