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 radial velocity instruments are modelled together with separate jitters and instrumental offsets.

using Octofitter, OctofitterRadialVelocity, Distributions, PlanetOrbits, CairoMakie
using Pigeons

gaia_id = 5164707970261890560

5164707970261890560

We build the bodies first. The host is an ordinary body, with a mass and-0-because Gaia and Hipparcos see a photocentre–- a flux in each of their bands. Setting the host's flux to 1 makes every other body's flux a contrast ratio. Set dark companions to 0.

A = Body(
    name="A",
    variables=@variables begin
        mass ~ truncated(Normal(0.82, 0.02), lower=0.5, upper=1.5) # M⊙ (Baines & Armstrong 2011)
        flux_G  = 1.0   # Gaia G
        flux_Hp = 1.0   # Hipparcos Hp
    end
)

b = Body(
    name="b",
    about=A,
    # No relative astrometry is included since the planet has not yet been
    # directly detected.
    variables=@variables begin
        # For speed of example, we are fitting a circular orbit only.
        e = 0
        ω = 0.0
        # Masses are solar masses; `mjup` is a plain multiplicative constant.
        mass ~ Uniform(0, 3mjup)      # M⊙
        a ~ Uniform(3, 10)            # AU
        i ~ Sine()
        Ω ~ Uniform(0, 2pi)
        # Phase: the mean anomaly at a reference epoch.
        M0 ~ Uniform(0, 2pi)
        epoch = 58849.0
        flux_G  = 0.0   # dark to Gaia and Hipparcos
        flux_Hp = 0.0
    end
)

We load in data from one RV instrument. We use MarginalizedRVObs instead of RadialVelocityObs to analytically marginalize out the radial velocity zero point of each instrument, saving one parameter.

hires_data = OctofitterRadialVelocity.HIRES_rvs("HD22049")
rvlike_hires = MarginalizedRVObs(
    hires_data;
    target=A, ref=Barycentre,   # the star's reflex motion against the barycentre
    name="HIRES",
    variables=@variables begin
        jitter ~ LogUniform(0.1, 100) # m/s
    end
)

We load the G23H data for this target. target= and blends= declare which bodies this source is modelled from:

pma = G23HObs(; gaia_id, target=A, blends=(b,), ref=Barycentre, freeze_epochs=true,
)

┌ Info: Count of missed or rejected transits:
└   dr3 = 1
┌ Warning: Gaia DR2 matched-transit count exceeds the geometric DR2-window pool; the excess must be doubly-downlinked transits and is modelled as repeated epochs.
│   n_pool = 22
│   n_dr2_total = 28
└ @ Octofitter ~/octo-maintenance/runner/actions-runner/_work/Octofitter.jl/Octofitter.jl/src/likelihoods/g23h.jl:1219
┌ Info: DR2/DR3 epoch selection
│   n2_win = 21
│   n_tail = 9
│   n_dr2_total = 28
└   n_dr2_distinct_range = (14, 22)

freeze_epochs=true fixes which of Gaia's forecast scans were actually used, instead of marginalizing over that selection. It is the fast-but-approximate setting (the same one the G23H tutorial recommends for exploration). Without it, the model gains one free parameter per forecast scan — 31 for this target. Set freeze_epochs=false (the default) for production fits.

The catalog row it fetched is available as pma.catalog, which is convenient for centring the frame priors below:

cat = pma.catalog
(cat.parallax, cat.pmra_dr3, cat.pmdec_dr3)

(310.5772928005821, -974.758145249517, 20.875840089774208)

Now the system. Absolute astrometry needs a full absolute frame: plx, ra, dec, pmra, pmdec, rv and ref_epoch must all be declared together in the system block (a partial frame is rejected at model-build time). pmra/pmdec are the reference point's proper motion, which is what G23HObs compares its catalog proper motions against.

sys = System(
    name="ϵEri",
    bodies=[A, b],
    observations=[pma, rvlike_hires],
    variables=@variables begin
        plx ~ gaia_plx(; gaia_id)
        ra  = $(cat.ra)
        dec = $(cat.dec)
        pmra  ~ Normal(cat.pmra_dr3, 10)
        pmdec ~ Normal(cat.pmdec_dr3, 10)
        rv = 0.0             # systemic radial velocity [km/s]
        ref_epoch = 57388.5  # Gaia DR3 reference epoch, J2016.0 (MJD)
    end
)
# Build model
model = Octofitter.LogDensityModel(sys)

LogDensityModel for System ϵEri of dimension 19 and 129 epochs with fields .ℓπcallback and .∇ℓπcallback

Find good starting points and visualize the starting position + data:

init_chain = initialize!(model)
octoplot(model, init_chain)

using Pigeons
results, pt = octofit_pigeons(model, n_rounds=10, n_chains=10, n_chains_variational=0, explorer=SliceSampler());

[ Info: Sampler running with multiple threads     : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: [ϵEri] observing_geometry = true (user)
[ Info: [ϵEri] barycentric_lighttime = true (user)
─────────────────────────────────────────────────────────────────────────────────────────────────────────────
  scans     restarts      Λ        time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
        2          0       3.25      0.581   8.45e+08  -1.41e+05          0      0.639          1          1
        4          0       5.32      0.868   1.55e+09  -4.86e+03          0      0.409          1          1
        8          0       6.02       1.51   3.06e+09  -6.32e+03          0      0.331          1          1
       16          0       6.71       2.74   5.93e+09  -3.11e+03          0      0.254          1          1
       32          0       6.45       5.11   1.15e+10  -1.41e+03          0      0.284          1          1
       64          4       2.46       12.1    2.6e+10       -256      0.384      0.726          1          1
      128          9       2.61       24.5   5.26e+10       -601       0.51       0.71          1          1
      256         28        2.7       49.5   1.06e+11       -601      0.566        0.7          1          1
      512         49       2.67         94    2.1e+11       -601      0.643      0.703          1          1
 1.02e+03        101       2.73        180   4.21e+11       -601      0.632      0.696          1          1
─────────────────────────────────────────────────────────────────────────────────────────────────────────────

We can now plot the results with a multi-panel plot. octoplot derives its panels from the model, so the sky orbit, the RV time series and its phase-folded panel all appear without being asked for:

octoplot(model, results)

rvplot is the complementary view. Where octoplot draws many samples from the posterior — one RV curve per draw, each carrying its own instrument offsets, with the measurements left exactly as reported — rvplot renders a single draw, which is what lets it put every instrument back on one axis:

rvplot(model, results)

We can see what the orbit looks like for the maximum a-posteriori sample (note, we would need to run an optimizer to get the true MAP value; this is just the MCMC sample with highest posterior density). Slicing the chain is how you plot a single draw:

i_max = argmax(vec(results[:logpost]))
res = octoplot(model, results[i_max:i_max, :, :])
Label(res.figure[0, 1], "Maximum a-posteriori orbit sample")
Makie.resize_to_layout!(res.figure)
res.figure

And a corner plot:

using PairPlots
octocorner(model, results, small=true)