Astrometry, PMA, and RV
Octofitter.jl supports fitting orbit models to astrometric motion in the form of GAIA-Hipparcos proper motion anomaly (HGCA; https://arxiv.org/abs/2105.11662). These data points are calculated by finding the difference between a long term proper motion of a star between the Hipparcos and GAIA catalogs, and their proper motion calculated within the windows of each catalog. This gives four data points that can constrain the dynamical mass & orbits of planetary companions (assuming we subtract out the net trend).
If your star of interest is in the HGCA, all you need is it's GAIA DR3 ID number. You can find this number by searching for your target on SIMBAD.
For this tutorial, we will examine the star and companion HD 91312 A & B discovered by SCExAO. We will use their published astrometry and proper motion anomaly extracted from the HGCA.
We will also perform a model comparison: we will fit the same model to four different subsets of data to see how each dataset are impacting the final constraints. This is an important consistency check, especially with proper motion / absolute astrometry data which be susceptible to systematic errors.
The first step is to find the GAIA source ID for your object. For HD 91312, SIMBAD tells us the GAIA DR3 ID is 756291174721509376.
Work through the the Proper Motion Anomaly tutorial before attempting a joint fit.
Model: PMA Only
Initial setup:
using Octofitter, Distributions, Random
using CairoMakie
using PairPlots
using PigeonsPrecompiling packages...
14785.6 ms ? IntervalArithmetic → IntervalArithmeticForwardDiffExt
7771.1 ms ? Unitful → ForwardDiffExt
22995.8 ms ? PlanetOrbits → PlanetOrbitsMakieExt
Info Given OctofitterPairPlotsExt was explicitly requested, output will be shown live
┌ Warning: Module Octofitter with build ID fafbfcfd-0ef1-fd64-1fb5-e2fe57ab2a86 is missing from the cache.
│ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2709
27816.4 ms ? Octofitter → OctofitterMakieExt
27517.2 ms ? Octofitter → OctofitterPairPlotsExt
┌ Warning: Module Octofitter with build ID fafbfcfd-0ef1-fd64-1fb5-e2fe57ab2a86 is missing from the cache.
│ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2709
Precompiling packages...
Info Given OctofitterPigeonsExt was explicitly requested, output will be shown live
┌ Warning: Module Octofitter with build ID fafbfcfd-0ef1-fd64-1fb5-e2fe57ab2a86 is missing from the cache.
│ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2709
55893.9 ms ? Octofitter → OctofitterPigeonsExt
118205.3 ms ? Pigeons → PigeonsForwardDiffExt
┌ Warning: Module Octofitter with build ID fafbfcfd-0ef1-fd64-1fb5-e2fe57ab2a86 is missing from the cache.
│ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2709
Precompiling packages...
Info Given PigeonsForwardDiffExt was explicitly requested, output will be shown live
┌ Warning: Module ForwardDiff with build ID fafbfcfd-ce51-eccc-9fc8-6019f98c1385 is missing from the cache.
│ This may mean ForwardDiff [f6369f11-7733-5829-9624-2563aa707210] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2709
3998.9 ms ? Pigeons → PigeonsForwardDiffExt
┌ Warning: Module ForwardDiff with build ID fafbfcfd-ce51-eccc-9fc8-6019f98c1385 is missing from the cache.
│ This may mean ForwardDiff [f6369f11-7733-5829-9624-2563aa707210] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2709We begin by finding orbits that are consistent with the astrometric motion. Later, we will add in relative astrometry to the fit from direct imaging to further constrain the planet's orbit and mass.
Compared to previous tutorials, we now add a prior on the mass of the companion, called mass. All masses are in solar masses — mjup is a plain multiplicative constant, so a Jupiter-mass prior reads LogUniform(0.5mjup, 1000mjup). There is also nothing to do to "place a prior on the host mass rather than the system total mass": every body carries its own mass, and an orbit's total mass comes from the hierarchy.
The bodies
A = Body(
name="A",
variables=@variables begin
mass ~ truncated(Normal(1.61, 0.1), lower=0.1) # Msol
flux_G = 1.0 # the host sets the contrast scale
flux_Hp = 1.0
end
)
planet_b = Body(
name="b",
about=A,
variables=@variables begin
mass ~ LogUniform(0.5mjup, 1000mjup) # Msol
a ~ LogUniform(0.1, 20)
e ~ Uniform(0, 0.999)
ω ~ Uniform(0, 2pi)
i ~ Sine() # the Sine() distribution is defined by Octofitter
Ω ~ Uniform(0, 2pi)
# `θ` (position angle) + `epoch` is a phase parametrization the orbit
# constructor understands directly.
θ ~ Uniform(0, 2pi)
epoch = 57423.0 # epoch of the GAIA measurement
flux_G = 0.0 # dark to Gaia…
flux_Hp = 0.0 # …and to Hipparcos
end
)The two flux_* lines are contrast ratios against the host — set the host to 1.0 and every other body's flux is a ratio — and they are how the model places the Gaia photocentre and modulates the Hipparcos abscissa. Give a luminous, unresolved companion a real prior here (flux_G ~ Uniform(0, 1)); omitting all four lines is legal and equivalent to 0.0.
Retrieving the HGCA
hgca_obs = HGCAObs(
gaia_id = 756291174721509376,
target = A,
blends = (planet_b,),
ref = Barycentre,
freeze_epochs = true, # fast but approximate; drop for a production fit
)
hgca_obs.table.kind6-element Vector{Symbol}:
:ra_hip
:dec_hip
:ra_hg
:dec_hg
:ra_dr3
:dec_dr3System Model & Specifying Proper Motion Anomaly
Now that we have our bodies, we create a system model to contain them. Observations are no longer attached to a planet — they are a flat list on the system, and each one names what it is a measurement of (target/host) and what it is measured against (ref).
We specify priors on plx as usual, but here we use the gaia_plx helper function to read the parallax and uncertainty directly from the Gaia DR3 catalog using its source ID.
We also add parameters for the star's long term proper motion. This is usually close to the long term trend between the Hipparcos and GAIA measurements.
When fitting HGCA data, use wide, uninformative priors for pmra and pmdec, such as Normal(0, 1000) (0 ± 1000 mas/yr). Do not use Gaia DR3 proper motion values as informative priors—this would double-count the information since the HGCA already incorporates Gaia astrometry and will constrain the system's proper motion through the likelihood. The pmra/pmdec parameters represent the center-of-mass proper motion, which the HGCA measurements help determine.
ra_deg = 158.30707896392835
dec_deg = 40.42555422701387
sys = System(
name="HD91312_pma",
bodies=[A, planet_b],
observations=[hgca_obs],
variables=@variables begin
plx ~ gaia_plx(gaia_id=756291174721509376)
# Priors on the center of mass proper motion
pmra ~ Normal(-137, 10)
pmdec ~ Normal(2, 10)
# The rest of the absolute frame
ra = $ra_deg
dec = $dec_deg
rv = 0.0 # m/s
ref_epoch = $(Octofitter.meta_gaia_DR3.ref_epoch_mjd)
end
)
model_pma = Octofitter.LogDensityModel(sys)LogDensityModel for System HD91312_pma of dimension 11 and 6 epochs with fields .ℓπcallback and .∇ℓπcallback
Sampling from the posterior (PMA only)
Because proper motion anomaly data is quite sparse, it can often produce multi-modal posteriors. If your orbit already has several relative astrometry or RV data points, this is less of an issue.
chain_pma, pt = octofit_pigeons(model_pma, n_rounds=8)
display(chain_pma)[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: [HD91312_pma] observing_geometry = true (user)
[ Info: [HD91312_pma] barycentric_lighttime = true (user)
[ Info: [HD91312_pma] observing_geometry = true (user)
[ Info: [HD91312_pma] barycentric_lighttime = true (user)
┌ Info: Starting values not provided for all parameters! Guessing starting point using global optimization:
│ num_params = 11
└ num_fixed = 0
┌ Warning: Verbosity toggle: unrecognized_stop_reason
│ Unrecognized stop reason: Too many steps (101) without any function evaluations (probably search has converged). Defaulting to ReturnCode.Default.
└ @ OptimizationBase ~/.julia/packages/OptimizationBase/Jfw5O/src/utils.jl:170
┌ Warning: Pathfinder error: MethodError(copy, (SplittableRandoms.SplittableRandom(0xc82883f44ae5e0cb, 0xaabf13ae126f9159),), 0x0000000000009f30)
└ @ Octofitter ~/octo-maintenance/runner/actions-runner/_work/Octofitter.jl/Octofitter.jl/src/initialization.jl:944
┌ Warning: Using BBO result as fallback
└ @ Octofitter ~/octo-maintenance/runner/actions-runner/_work/Octofitter.jl/Octofitter.jl/src/initialization.jl:945
┌ Info: Found sample of initial positions
│ logpost_range = (-9.803133125692725, -9.803133125692725)
└ mean_logpost = -9.803133125692739
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 2.87 3.11 2.39 2.03e+09 -3.37e+04 0 0.807 1 1
4 0 4.34 4.22 1.35 3.6e+09 -1.03e+03 0 0.724 1 1
8 0 4.23 4.75 2.19 6.53e+09 -258 2.67e-218 0.71 1 1
16 0 4.81 5.6 4.13 1.19e+10 -59.3 6.33e-44 0.664 0.996 1
32 0 6.62 6.18 7.99 2.26e+10 -16.5 0.0581 0.587 1 1
64 0 7.27 6.17 16.9 4.73e+10 -17.4 0.000271 0.566 1 1
128 2 7.44 5.85 31.5 9.05e+10 -21.2 0.0543 0.571 1 1
256 8 7.83 6.13 59 1.76e+11 -22.2 0.213 0.55 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────Pair Plot
If we wish to examine the covariance between parameters in more detail, we can construct a pair-plot (aka. corner plot).
# Create a corner plot / pair plot.
# We can access any property from the chain specified in Variables
octocorner(model_pma, chain_pma, small=true)
Model: PMA & Relative Astrometry
The first orbit fit to only Hipparcos/GAIA data was very unconstrained. We will now add six epochs of relative astrometry (measured from direct images) gathered from the discovery paper.
astrom_dat = Table(;
epoch = [mjd("2016-12-15"), mjd("2017-03-12"), mjd("2017-03-13"), mjd("2018-02-08"), mjd("2018-11-28"), mjd("2018-12-15")],
ra = [133., 126., 127., 083., 058., 056.],
dec = [-174., -176., -172., -133., -122., -104.],
σ_ra = [07.0, 04.0, 04.0, 10.0, 10.0, 08.0],
σ_dec = [07.0, 04.0, 04.0, 10.0, 20.0, 08.0],
cor = [0.2, 0.3, 0.1, 0.4, 0.3, 0.2]
)
astrom_obs = RelAstromObs(
astrom_dat,
target = planet_b, # the body the offsets are of
ref = A, # …and the body they are measured against
name = "SCExAO",
variables = @variables begin
# Fixed values for this example - could be free variables:
jitter = 0 # mas [could use: jitter ~ Uniform(0, 10)]
northangle = 0 # radians [could use: northangle ~ Normal(0, deg2rad(1))]
platescale = 1 # relative [could use: platescale ~ truncated(Normal(1, 0.01), lower=0)]
end
)
scatter(astrom_obs.table.ra, astrom_obs.table.dec)
We use the same bodies as before, but now condition the model on the astrometry by adding it to the system's observations list. We wrap it in ObsPriorONeil2019 — the observable-based prior of O'Neil et al. (2019), which reweights the orbit prior by the Jacobian of the observables.
List ObsPriorONeil2019(astrom_obs) in observations=, not the wrapper and astrom_obs. Listing both counts the astrometry twice.
The wrapper needs to know which orbit the prior applies to. It defaults to the wrapped observation's target, which is correct here and for relative RVs. A stellar-reflex RadialVelocityObs(…; target=A, ref=Barycentre) has no orbit of its own, so wrapping one requires ObsPriorONeil2019(rvs; orbit=planet_b) explicitly (orbit=(b, c) sums over several orbits).
sys_astrom = System(
name="HD91312_pma_astrom",
bodies=[A, planet_b],
observations=[hgca_obs, ObsPriorONeil2019(astrom_obs)],
variables=@variables begin
plx ~ gaia_plx(gaia_id=756291174721509376)
pmra ~ Normal(-137, 10)
pmdec ~ Normal(2, 10)
ra = $ra_deg
dec = $dec_deg
rv = 0.0
ref_epoch = $(Octofitter.meta_gaia_DR3.ref_epoch_mjd)
end
)
model_pma_astrom = Octofitter.LogDensityModel(sys_astrom, verbosity=4)
chain_pma_astrom, pt = octofit_pigeons(model_pma_astrom, n_rounds=7, explorer=SliceSampler())[ Info: [HD91312_pma_astrom] observing_geometry = true (auto): "HGCA" does not report comparable predictions — over 300 prior draws (seed 0xc70f177e5000001)
[ Info: [HD91312_pma_astrom] barycentric_lighttime = true (auto): "HGCA" does not report comparable predictions — over 300 prior draws (seed 0xc70f177e5000001)
[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 11
┌ Info: Determined number type
└ T = Float64
ℓπcallback(θ): 0.000082 seconds (42 allocations: 73.719 KiB)
∇ℓπcallback(θ): 1.020920 seconds (2.24 M allocations: 97.461 MiB, 99.94% compilation time)
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: [HD91312_pma_astrom] observing_geometry = true (user)
[ Info: [HD91312_pma_astrom] barycentric_lighttime = true (user)
[ Info: [HD91312_pma_astrom] observing_geometry = true (user)
[ Info: [HD91312_pma_astrom] barycentric_lighttime = true (user)
┌ Info: Starting values not provided for all parameters! Guessing starting point using global optimization:
│ num_params = 11
└ num_fixed = 0
┌ Warning: Verbosity toggle: unrecognized_stop_reason
│ Unrecognized stop reason: Too many steps (101) without any function evaluations (probably search has converged). Defaulting to ReturnCode.Default.
└ @ OptimizationBase ~/.julia/packages/OptimizationBase/Jfw5O/src/utils.jl:170
┌ Warning: Pathfinder error: MethodError(copy, (SplittableRandoms.SplittableRandom(0xc82883f44ae5e0cb, 0xaabf13ae126f9159),), 0x0000000000009f59)
└ @ Octofitter ~/octo-maintenance/runner/actions-runner/_work/Octofitter.jl/Octofitter.jl/src/initialization.jl:944
┌ Warning: Using BBO result as fallback
└ @ Octofitter ~/octo-maintenance/runner/actions-runner/_work/Octofitter.jl/Octofitter.jl/src/initialization.jl:945
┌ Info: Found sample of initial positions
│ logpost_range = (-95.38883728836572, -95.38883728836572)
└ mean_logpost = -95.38883728836574
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 3.28 3.61 2.04 1.89e+09 -3.39e+04 0 0.778 1 1
4 0 4.5 6.65 1.28 3.44e+09 -159 0 0.64 1 1
8 0 5.97 7.24 2.22 6.27e+09 -109 4.4e-41 0.574 0.996 1
16 0 6.83 8.38 4.41 1.17e+10 -98 7.15e-15 0.51 1 1
32 0 8.84 9.13 8.49 2.25e+10 -83.6 0.00977 0.42 1 1
64 0 9.02 5.34 18.4 4.78e+10 -91.6 0.0006 0.537 1 1
128 0 9.24 6.85 35.8 9.95e+10 -80 0.0505 0.481 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────octoplot(model_pma_astrom, chain_pma_astrom)
octoplot draws one time-series panel per plottable data channel: RelAstromObs contributes separation and position angle (or Δα⋆/Δδ), and G23HObs contributes pmra and pmdec — the five catalog proper motions against the model's reflex curve, with each point's mission averaging window drawn as a horizontal bar.
Wrapping an observation changes nothing about how it is drawn: ObsPriorONeil2019 reweights the prior, and the astrometry inside it still contributes its points to the sky panel and its own time-series panels, calibrated by the fitted platescale/northangle exactly as an unwrapped dataset is. The panels are labelled with the wrapper's name (obspri_SCExAO rather than SCExAO), which is also the name its calibration parameters take in the chain.
Model: PMA & Relative Astrometry & RVs
We now add in three additional epochs of stellar RVs.
rv_dat_abs = Table(;
epoch = [mjd("2008-05-01"), mjd("2010-02-15"), mjd("2016-03-01")],
rv = [1300, 700, -2700],
σ_rv = [150, 150, 150]
)
rvlike = RadialVelocityObs(
rv_dat_abs,
target = A, # the star's own reflex motion…
ref = Barycentre, # …against the system barycentre
name = "SOPHIE",
# Per dataset, and the default: with a full absolute frame like this one,
# the secular (perspective) acceleration term is non-zero and modelled.
# Pass `:data_corrected` for a series whose pipeline already removed it.
secular_acceleration = :model,
variables = @variables begin
jitter ~ truncated(Normal(10, 5), lower=0) # m/s
offset ~ Normal(0, 1000) # m/s
end
)RadialVelocityObs lives in core Octofitter and covers both cases: target=A, ref=Barycentre is the star's reflex motion against the barycentre, and target=b, ref=A is a companion's velocity relative to its host. No using OctofitterRadialVelocity is needed for either — that package is for Celerite, MarginalizedRVObs, and the archive loaders.
sys_rv_astrom = System(
name="HD91312_pma_rv_astrom",
bodies=[A, planet_b],
observations=[hgca_obs, rvlike, ObsPriorONeil2019(astrom_obs)],
variables=@variables begin
plx ~ gaia_plx(gaia_id=756291174721509376)
# Priors on the centre of mass proper motion
pmra ~ Normal(-137, 10)
pmdec ~ Normal(2, 10)
ra = $ra_deg
dec = $dec_deg
rv = 0.0 # m/s
ref_epoch = $(Octofitter.meta_gaia_DR3.ref_epoch_mjd)
end
)
model_pma_rv_astrom = Octofitter.LogDensityModel(sys_rv_astrom, verbosity=4)
chain_pma_rv_astrom, pt = octofit_pigeons(model_pma_rv_astrom, n_rounds=9, explorer=SliceSampler())
display(chain_pma_rv_astrom)[ Info: [HD91312_pma_rv_astrom] observing_geometry = true (auto): "HGCA" does not report comparable predictions — over 300 prior draws (seed 0xc70f177e5000001)
[ Info: [HD91312_pma_rv_astrom] barycentric_lighttime = true (auto): "HGCA" does not report comparable predictions — over 300 prior draws (seed 0xc70f177e5000001)
[ Info: [HD91312_pma_rv_astrom] "SOPHIE": Einstein term, peak-to-peak [m/s] ≈ 0.0244 (0.000163× its tightest σ = 150.0)
[ Info: [HD91312_pma_rv_astrom] "SOPHIE": secular-acceleration drift, peak-to-peak [m/s] ≈ 0.117 (0.000777× its tightest σ = 150.0)
[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 13
┌ Info: Determined number type
└ T = Float64
ℓπcallback(θ): 0.000077 seconds (42 allocations: 73.719 KiB)
∇ℓπcallback(θ): 1.408093 seconds (2.32 M allocations: 100.188 MiB, 18.07% gc time, 99.96% compilation time)
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: [HD91312_pma_rv_astrom] observing_geometry = true (user)
[ Info: [HD91312_pma_rv_astrom] barycentric_lighttime = true (user)
[ Info: [HD91312_pma_rv_astrom] observing_geometry = true (user)
[ Info: [HD91312_pma_rv_astrom] barycentric_lighttime = true (user)
┌ Info: Starting values not provided for all parameters! Guessing starting point using global optimization:
│ num_params = 13
└ num_fixed = 0
┌ Warning: Pathfinder error: MethodError(copy, (SplittableRandoms.SplittableRandom(0x01e2d3f16bc10c5f, 0xaabf13ae126f9159),), 0x0000000000009f82)
└ @ Octofitter ~/octo-maintenance/runner/actions-runner/_work/Octofitter.jl/Octofitter.jl/src/initialization.jl:944
┌ Warning: Using BBO result as fallback
└ @ Octofitter ~/octo-maintenance/runner/actions-runner/_work/Octofitter.jl/Octofitter.jl/src/initialization.jl:945
┌ Info: Found sample of initial positions
│ logpost_range = (-247.97027917704654, -247.97027917704654)
└ mean_logpost = -247.97027917704673
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 3.3 3.79 2.49 2.63e+09 -8.88e+03 0 0.771 1 1
4 0 5.64 4.21 1.83 4.83e+09 -801 0 0.682 1 1
8 0 7.12 5.63 3.08 8.79e+09 -414 6.45e-154 0.589 1 1
16 0 9.41 6.92 6.11 1.66e+10 -322 1.35e-66 0.473 1 1
32 0 9.65 8.38 11.9 3.18e+10 -243 0.0156 0.419 1 1
64 0 9.42 7.86 24.6 6.54e+10 -258 0.0132 0.442 1 1
128 0 10.7 9.29 45.4 1.29e+11 -182 0.0308 0.354 1 1
256 4 11.2 6.61 88.3 2.55e+11 -266 0.0489 0.426 0.999 1
512 19 10.4 5.52 179 5.16e+11 -92 0.154 0.486 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────The mass vs. semi-major axis posterior is now much more constrained:
pairplot(
(; a=chain_pma_rv_astrom["b_a"][:], mass=chain_pma_rv_astrom["b_mass"][:] ./ mjup) =>
(
PairPlots.Scatter(color=:red, markersize=5),
PairPlots.MarginHist(),
PairPlots.MarginQuantileText()
),
labels=Dict(:mass=>"mass [Mⱼᵤₚ]", :a=>"sma. [au]"),
)
It is now useful to display the orbits projected onto the plane of the sky using octoplot. This produces a sky panel plus one time-series panel per data channel — here, separation and position angle from the astrometry, and radial velocity (with a phase-folded panel) from the RVs.
octoplot(model_pma_rv_astrom, chain_pma_rv_astrom)
Model: Relative Astrometry & RVs (no PMA)
There is a final model we should consider: one using the RV and astrometry data, but not the proper motion anomaly. With no absolute astrometry there is no absolute frame to declare, so the system block shrinks to plx alone:
sys_final = System(
name="HD91312_rv_astrom",
bodies=[A, planet_b],
observations=[rvlike, ObsPriorONeil2019(astrom_obs)],
variables=@variables begin
plx ~ gaia_plx(gaia_id=756291174721509376)
end
)
model_rv_astrom = Octofitter.LogDensityModel(sys_final, verbosity=4)
chain_rv_astrom, pt = octofit_pigeons(model_rv_astrom, n_rounds=11)[ Info: [HD91312_rv_astrom] observing_geometry = false (auto): worst accumulated bias 0.038σ, 2.63× inside the 0.1σ limit — over 300 prior draws (seed 0xc70f177e5000001)
[ Info: [HD91312_rv_astrom] barycentric_lighttime = false (auto): changes no prediction at all — over 300 prior draws (seed 0xc70f177e5000001)
[ Info: [HD91312_rv_astrom] "SOPHIE": Einstein term, peak-to-peak [m/s] ≈ 0.0177 (0.000118× its tightest σ = 150.0)
[ Info: [HD91312_rv_astrom] "SOPHIE": secular-acceleration drift, peak-to-peak [m/s] ≈ 0.0 (0.0× its tightest σ = 150.0)
[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 11
┌ Info: Determined number type
└ T = Float64
ℓπcallback(θ): 0.000015 seconds
∇ℓπcallback(θ): 0.000021 seconds (1 allocation: 32 bytes)
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: [HD91312_rv_astrom] observing_geometry = false (user)
[ Info: [HD91312_rv_astrom] barycentric_lighttime = false (user)
[ Info: [HD91312_rv_astrom] observing_geometry = false (user)
[ Info: [HD91312_rv_astrom] barycentric_lighttime = false (user)
┌ Info: Starting values not provided for all parameters! Guessing starting point using global optimization:
│ num_params = 11
└ num_fixed = 0
┌ Warning: Pathfinder error: MethodError(copy, (SplittableRandoms.SplittableRandom(0x4612a50ef2b9dec0, 0xaabf13ae126f9159),), 0x0000000000009f97)
└ @ Octofitter ~/octo-maintenance/runner/actions-runner/_work/Octofitter.jl/Octofitter.jl/src/initialization.jl:944
┌ Warning: Using BBO result as fallback
└ @ Octofitter ~/octo-maintenance/runner/actions-runner/_work/Octofitter.jl/Octofitter.jl/src/initialization.jl:945
┌ Info: Found sample of initial positions
│ logpost_range = (-130.42610000303077, -130.42610000303077)
└ mean_logpost = -130.4261000030308
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 3.52 2.28 1.09 6.09e+06 -321 2.21e-187 0.813 1 1
4 0 3.19 4.88 0.157 2.11e+05 -127 2.18e-22 0.74 1 1
8 0 6.06 5.77 0.295 3.36e+05 -122 0.00354 0.619 1 1
16 0 4.56 5.7 0.59 6.9e+05 -117 0.145 0.669 1 1
32 0 6.94 7.16 0.915 1.03e+06 -90.4 0.0693 0.545 1 1
64 1 7.18 5.07 2.24 5.42e+07 -61.6 0.226 0.605 1 1
128 5 7.08 5.17 2.7 8.94e+07 -75.1 0.351 0.605 1 1
256 16 7.67 3.36 5.24 1.79e+08 -74.9 0.348 0.644 1 1
512 45 7.34 3.23 10.6 3.61e+08 -75 0.413 0.659 1 1
1.02e+03 94 7.51 3.61 21.3 7.23e+08 -75 0.449 0.641 1 1
2.05e+03 212 7.41 3.5 42.5 1.44e+09 -74.9 0.479 0.648 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────octoplot(model_rv_astrom, chain_rv_astrom)
Model Comparison
Let's now display the constraints provided by each data set in a single corner plot
octocorner(
model_pma,
chain_pma,
chain_pma_astrom,
chain_rv_astrom,
chain_pma_rv_astrom,
small=true,
axis=(;
b_a = (;lims=(low=0, high=25))
),
viz=(
PairPlots.MarginDensity(),
PairPlots.Scatter()
)
)
We see that the constraints provided by the PMA, the astrometry, and the radial velocity data all individually overlap, and agree with the joint model constraint. This is means that none of the datasets are in tension with each other, which might suggest an issue with the data or with the modelling assumptions (e.g. single planet).