Fit RV and Astrometry

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.

Note

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
    e = 0
    ω=0.0
    mass ~ Uniform(0, 3)
    a~Uniform(1, 5)
    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),
#     ...
# )

# 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[:])


@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)

    rv0_1 ~ Normal(0,10)
    rv0_2 ~ Normal(0,10)
    jitter_1 ~ truncated(Normal(0,10),lower=0)
    jitter_2 ~ truncated(Normal(0,10),lower=0)
end HGCALikelihood(;gaia_id) rvs_merged b

# Build model
model = Octofitter.LogDensityModel(ϵEri; autodiff=:ForwardDiff, verbosity=4) # defaults are ForwardDiff, and verbosity=0
LogDensityModel for System ϵEri of dimension 15 with fields .ℓπcallback and .∇ℓπcallback

Now sample:

using Random
rng = Random.Xoshiro(0)

results = octofit(rng, model)
Chains MCMC chain (1000×34×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 117.89 seconds
Compute duration  = 117.89 seconds
parameters        = M, plx, pmra, pmdec, rv0_1, rv0_2, jitter_1, jitter_2, b_mass, b_a, b_i, b_Ωy, b_Ωx, b_τy, b_τx, b_e, 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       0.7802    0.0103    0.0003   1291.9677   675.7607    1.0001 ⋯
         plx     310.5781    0.1383    0.0040   1174.8248   760.0771    1.0069 ⋯
        pmra    -975.0702    0.0202    0.0005   1646.4091   840.7347    1.0061 ⋯
       pmdec      20.0405    0.0165    0.0004   1424.7227   688.9930    1.0008 ⋯
       rv0_1       3.6867    0.5305    0.0138   1474.6597   656.6717    0.9997 ⋯
       rv0_2       0.8096    1.0489    0.0292   1301.9413   602.5053    1.0019 ⋯
    jitter_1       5.5702    0.1467    0.0038   1441.8613   734.4811    1.0047 ⋯
    jitter_2      11.1887    0.7843    0.0225   1305.7107   734.6116    0.9991 ⋯
      b_mass       0.5794    0.1090    0.0031   1284.9795   874.2977    1.0013 ⋯
         b_a       1.7317    0.0089    0.0002   1408.2346   665.9359    0.9995 ⋯
         b_i       2.7171    0.1080    0.0031   1399.6216   842.7230    1.0004 ⋯
        b_Ωy       0.9824    0.1119    0.0032   1242.0733   660.1217    1.0005 ⋯
        b_Ωx       0.0438    0.2265    0.0068   1121.0058   642.2759    0.9994 ⋯
        b_τy       0.4786    0.1052    0.0030   1232.7972   641.2952    1.0024 ⋯
        b_τx       0.8790    0.0967    0.0025   1467.6900   699.6500    0.9993 ⋯
         b_e       0.0000    0.0000       NaN         NaN        NaN       NaN ⋯
         b_ω       0.0000    0.0000       NaN         NaN        NaN       NaN ⋯
      ⋮            ⋮           ⋮         ⋮          ⋮          ⋮          ⋮    ⋱
                                                     1 column and 4 rows omitted

Quantiles
  parameters         2.5%        25.0%        50.0%        75.0%        97.5%
      Symbol      Float64      Float64      Float64      Float64      Float64

           M       0.7603       0.7734       0.7800       0.7867       0.8018
         plx     310.3184     310.4770     310.5778     310.6707     310.8458
        pmra    -975.1110    -975.0826    -975.0699    -975.0576    -975.0294
       pmdec      20.0092      20.0286      20.0409      20.0515      20.0738
       rv0_1       2.6451       3.3340       3.7007       4.0592       4.6622
       rv0_2      -1.2968       0.1269       0.8243       1.5293       2.8630
    jitter_1       5.2813       5.4731       5.5650       5.6597       5.8659
    jitter_2       9.7894      10.6481      11.1342      11.6672      12.8108
      b_mass       0.3660       0.5021       0.5798       0.6565       0.7836
         b_a       1.7145       1.7255       1.7317       1.7378       1.7487
         b_i       2.4254       2.6625       2.7364       2.7920       2.8675
        b_Ωy       0.7790       0.9053       0.9729       1.0549       1.2026
        b_Ωx      -0.4179      -0.1080       0.0440       0.2025       0.4817
        b_τy       0.2671       0.4100       0.4792       0.5503       0.6762
        b_τx       0.6975       0.8148       0.8772       0.9410       1.0727
         b_e       0.0000       0.0000       0.0000       0.0000       0.0000
         b_ω       0.0000       0.0000       0.0000       0.0000       0.0000
      ⋮            ⋮            ⋮            ⋮            ⋮            ⋮
                                                                 4 rows omitted

We can now plot the results with a multi-panel plot:

## Save results plot
octoplot(model, results)
Example block output

We can also plot just the RV fit:

fig = OctofitterRadialVelocity.rvpostplot(model, results)
Example block output

And a corner plot:

using CairoMakie, PairPlots
octocorner(model, results, small=true)
Example block output