Fit RV and Proper Motion Anomaly
In this example, we will fit an orbit model to a combination of radial velocity and Hipparcos-GAIA proper motion anomaly for the star $\epsilon$ Eridani. We will use some of the radial velocity data collated in Mawet et al 2019.
Radial velocity modelling is supported in Octofitter via the extension package OctofitterRadialVelocity. To install it, run pkg> add OctofitterRadialVelocity
Datasets from two different radial velocity insturments are included and modelled together with separate jitters and instrumental offsets.
using Octofitter, OctofitterRadialVelocity, Distributions, PlanetOrbits, CairoMakie
gaia_id = 5164707970261890560
@planet b Visual{KepOrbit} begin
# For speed of example, we are fitting a circular orbit only.s
e = 0
ω=0.0
mass ~ Uniform(0, 3)
a ~ truncated(Normal(3.48,0.1),lower=0)
i ~ Sine()
Ω ~ UniformCircular()
τ ~ UniformCircular(1.0)
P = √(b.a^3/system.M)
tp = b.τ*b.P*365.25 + 58849 # reference epoch for τ. Choose an MJD date near your data.
end # No planet astrometry is included since it has not yet been directly detected
# Convert from JD to MJD
# Data tabulated from Mawet et al
jd(mjd) = mjd - 2400000.5
# # You could also specify it manually like so:
# rvs = StarAbsoluteRVLikelihood(
# (;inst_idx=1, epoch=jd(2455110.97985), rv=−6.54, σ_rv=1.30),
# (;inst_idx=1, epoch=jd(2455171.90825), rv=−3.33, σ_rv=1.09), # units in meters/second
# ...
# )
# We will load in data from two RV
rvharps = OctofitterRadialVelocity.HARPS_RVBank_rvs("HD22049")
rvharps.table.inst_idx .= 1
rvhires = OctofitterRadialVelocity.HIRES_rvs("HD22049")
rvhires.table.inst_idx .= 2
rvs_merged = StarAbsoluteRVLikelihood(
cat(rvhires.table, rvharps.table,dims=1),
instrument_names=["HARPS", "HIRES"]
)
scatter(
rvs_merged.table.epoch[:],
rvs_merged.table.rv[:],
axis=(
xlabel="time [mjd]",
ylabel="RV [m/s]",
)
)
We load the HGCA data for this target:
hgca_like = HGCALikelihood(;gaia_id, N_ave=1)HGCALikelihood Table with 3 columns and 4 rows:
epoch meas inst
┌────────────────────
1 │ 48312.2 ra hip
2 │ 48423.4 dec hip
3 │ 57322.8 ra gaia
4 │ 57278.7 dec gaiaIn the interests of time, we set N_ave=1 to speed up the computation. This parameter controls how the model smears out the simulated Gaia and Hipparcos measurements. For a real target, leave it at the default value once you have completed testing.
@system ϵEri begin
M ~ truncated(Normal(0.78, 0.01),lower=0)
plx ~ gaia_plx(;gaia_id)
pmra ~ Normal(-975, 10)
pmdec ~ Normal(20, 10)
# Radial velocity zero point per instrument
rv0 ~ Product([
Normal(0,10), # m/s
Normal(0,10),
])
# Radial jitter per instrument.
jitter ~ Product([
truncated(Normal(0,10),lower=0), # m/s
truncated(Normal(0,10),lower=0),
])
# You can also set both instruments to the same jitter, eg by putting instead (with = not ~):
# jitter_2 = system.jitter_1
end hgca_like rvs_merged b
# Build model
model = Octofitter.LogDensityModel(ϵEri)LogDensityModel for System ϵEri of dimension 15 and 884 epochs with fields .ℓπcallback and .∇ℓπcallback
Now sample. You could use HMC via octofit or tempered sampling via octofit_pigeons. When using tempered sampling, make sure to start julia with julia --thread=auto. Each additional round doubles the number of posterior samples, so n_rounds=10 gives 1024 samples. You should adjust n_chains to be roughly double the Λ value printed out during sample, and n_chains_variational to be roughly double the Λ_var column.
using Pigeons
results, pt = octofit_pigeons(model, n_rounds=8, n_chains=10, n_chains_variational=10);[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 15
ℓπcallback(θ, args...): 0.000052 seconds (43 allocations: 1.625 KiB)
┌ Info: Tuning autodiff
│ chunk_size = 15
└ t = 3.997e-6
┌ Info: Selected auto-diff chunk size
└ ideal_chunk_size = 15
∇ℓπcallback(θ): 0.000040 seconds (34 allocations: 10.016 KiB)
[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 15
ℓπcallback(θ, args...): 0.000005 seconds (43 allocations: 1.625 KiB)
┌ Info: Tuning autodiff
│ chunk_size = 15
└ t = 4.248e-6
┌ Info: Selected auto-diff chunk size
└ ideal_chunk_size = 15
∇ℓπcallback(θ): 0.000005 seconds (34 allocations: 10.016 KiB)
[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-2883.922757647466, -2873.745439258521)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 5.87 4.21 17.8 8.18e+08 -3.46e+04 0 0.47 0.947 0.973
4 0 4.19 5.13 19.8 5.96e+08 -4.57e+03 0 0.509 0.897 0.938
8 0 4.36 5.49 28.7 8.58e+08 -3.23e+03 0 0.482 0.931 0.944
16 0 5.87 5.63 59 1.75e+09 -3.8e+03 0 0.395 0.924 0.938
32 0 6.72 6.27 118 3.55e+09 -3.08e+03 5.22e-164 0.317 0.925 0.935
64 10 6.52 2.13 252 5.98e+09 -2.9e+03 4.93e-31 0.545 0.923 0.938
128 19 6.86 2.15 485 1.16e+10 -2.9e+03 1.12e-29 0.526 0.92 0.935
256 30 7.08 2.26 970 2.33e+10 -2.9e+03 0.000135 0.509 0.926 0.936
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────We can now plot the results with a multi-panel plot:
octoplot(model, results[1:200,:,:], show_mass=true)
We can also plot just the RV curve from the maximum a-posteriori fit:
fig = OctofitterRadialVelocity.rvpostplot(model, results)
We can see what the visual orbit looks like for the maximum a-posteriori orbit:
i_max = argmax(results[:logpost][:])
fig = octoplot(
model,
results[i_max,:,:],
# change the colour map a bit:
colormap=Makie.cgrad([Makie.wong_colors()[1], "#FAFAFA"]),
show_astrom=true,
show_astrom_time=false,
show_rv=false,
show_hgca=false,
mark_epochs_mjd=[
mjd("2030")
]
)
Label(fig[0,1], "Maximum a-posteriori orbit sample")
fig
And a corner plot:
using CairoMakie, PairPlots
octocorner(model, results, small=true)