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.

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]",
    )
)

@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 HGCALikelihood(;gaia_id) rvs_merged b

# Build model
model = Octofitter.LogDensityModel(ϵEri)

LogDensityModel for System ϵEri of dimension 15 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);

┌ Warning:
│ This model's log density function is not type stable, but all of its components are.
│ This may indicate a performance bug in Octofitter; please consider filing an issue on GitHub.
└ @ Octofitter ~/work/Octofitter.jl/Octofitter.jl/src/logdensitymodel.jl:294
┌ Warning:
│ This model's log density function is not type stable, but all of its components are.
│ This may indicate a performance bug in Octofitter; please consider filing an issue on GitHub.
└ @ Octofitter ~/work/Octofitter.jl/Octofitter.jl/src/logdensitymodel.jl:294
[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└   initial_logpost_range = (-2896.1593155258165, -2882.438484925881)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  scans     restarts      Λ        Λ_var      time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
        2          0        3.4       3.85         14   4.06e+09  -2.06e+04          0      0.618      0.964      0.979
        4          0       4.55       4.72       15.2   5.98e+09  -8.39e+03          0      0.512      0.913      0.943
        8          0       5.44       4.99       21.6   8.22e+09  -4.64e+03          0      0.451      0.919      0.944
       16          0       5.83        5.3         47   1.82e+10  -3.53e+03          0      0.415      0.924      0.942
       32          0       6.22       5.92       92.4   3.58e+10  -2.91e+03   2.34e-77      0.361       0.93       0.94
       64          8       6.72       2.23        188    7.1e+10   -2.9e+03   2.08e-40      0.529      0.925      0.939
      128         15       6.72       2.06        379   1.43e+11   -2.9e+03   1.46e-38      0.538      0.929      0.938
      256         26       6.88       2.23        753   2.85e+11   -2.9e+03    8.1e-05      0.521      0.925      0.935
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

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)