Fitting Gaia DR4 IAD
This tutorial shows how you can use Octofitter to fit preliminary RV and absolute astrometry data from DR4, using various already published data, including:
- Gaia-BH3 black hole (the only real DR4 data released so far)
- Data from https://dace.unige.ch/openData/?record=10.82180/dace-gaia-ohp
The final format of the Gaia IAD data may change, in which case this tutorial will be updated.
using CairoMakie
using Octofitter
using Distributions
using CSV, DataFramesOHP Gaia Splinter Session
Download and plot the data:
fname = download("https://dace.unige.ch/downloads/open_data/dace-gaia-ohp/files/target_1.csv")
df = CSV.read(fname, DataFrame);df = CSV.read(joinpath(@__DIR__, "target_1.csv"), DataFrame)We can now construct a likelihood object for this data. We must also supply the Gaia ID, which will be used to query the full Gaia solution for this object (for now, using DR3):
ref_epoch_mjd = 57936.375
gaia_dr4_obs = GaiaDR4AstromObs(
df,
# For plotting reasons, you need to supply a Gaia ID to know e.g. the absolute Ra and Dec
# For these ~simulated examples, you can pick one at random
gaia_id=4373465352415301632,
variables=@variables begin
astrometric_jitter ~ LogUniform(0.00001, 10) # mas
ra_offset_mas ~ Normal(0, 10000)
dec_offset_mas ~ Normal(0, 10000)
pmra ~ Uniform(-1000, 1000) # mas/yr
pmdec ~ Uniform(-1000, 1000) # mas/yr
# ra_offset_mas = -0.00154
# dec_offset_mas = 0.0019354
# pmra = 5.421
# pmdec = -24.121
plx = system.plx
ref_epoch = $ref_epoch_mjd
end
)GaiaDR4AstromObs Table with 8 columns and 102 rows:
obs_time_tcb centroid_pos_al centroid_pos_error_… parallax_factor_al ⋯
┌───────────────────────────────────────────────────────────────────────────
1 │ 2.45704e6 -48.1529 0.044486 0.694068 ⋯
2 │ 2.45706e6 -40.9114 0.0447412 -0.583294 ⋯
3 │ 2.45716e6 36.8968 0.0441832 0.36178 ⋯
4 │ 2.45716e6 36.6379 0.0488587 0.358281 ⋯
5 │ 2.45717e6 -20.7433 0.061433 -0.522976 ⋯
6 │ 2.45717e6 -20.9825 0.0470352 -0.526464 ⋯
7 │ 2.45721e6 60.4545 0.0409268 0.7256 ⋯
8 │ 2.45721e6 60.3531 0.0423611 0.725443 ⋯
9 │ 2.45724e6 14.2189 0.0413549 -0.611316 ⋯
10 │ 2.45724e6 14.3917 0.0447546 -0.609065 ⋯
11 │ 2.45727e6 45.7711 0.0384871 0.385425 ⋯
12 │ 2.45727e6 45.7168 0.0406302 0.387544 ⋯
13 │ 2.45733e6 9.56681 0.0398049 -0.162231 ⋯
14 │ 2.45733e6 9.15272 0.0394856 -0.156233 ⋯
15 │ 2.45733e6 9.0219 0.0425311 -0.153719 ⋯
16 │ 2.45733e6 8.60603 0.0402893 -0.147658 ⋯
17 │ 2.45733e6 8.43312 0.0416049 -0.145114 ⋯
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋱mjup2msol = Octofitter.mjup2msol
ref_epoch_mjd = Octofitter.meta_gaia_DR3.ref_epoch_mjd
orbit_ref_epoch = mean(gaia_dr4_obs.table.epoch)
b = Planet(
name="b",
basis=Visual{KepOrbit},
observations=[],
variables=@variables begin
a ~ LogUniform(0.01, 100)
e ~ Uniform(0, 0.99)
ω ~ Uniform(0,2pi)
i ~ Sine()
Ω ~ Uniform(0,2pi)
θ ~ Uniform(0,2pi)
tp = θ_at_epoch_to_tperi(θ, $orbit_ref_epoch; M=system.M, e, a, i, ω, Ω)
mass ~ LogUniform(0.01, 1000)
end
)
sys = System(
name="target_1",
companions=[b],
observations=[gaia_dr4_obs,],
variables=@variables begin
M = 1.0
# Note: keep these physically plausible to prevent numerical errors
plx ~ Uniform(0.01,100) # mas
end
)
model = Octofitter.LogDensityModel(sys, verbosity=4)LogDensityModel for System target_1 of dimension 13 and 102 epochs with fields .ℓπcallback and .∇ℓπcallback
Find a starting point via global optimization and variational approximation, and plot the initial points against the data:
init_chain = initialize!(model)
octoplot(model, init_chain)Sample using parallel tempering (could also use HMC for these unimodel distributions):
using Pigeons
chain, pt = octofit_pigeons(model, n_rounds=10)(chain = MCMC chain (1024×21×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))Plot the results:
octoplot(model, chain)If we pick an individual draw, we can also plot the orbit against the Gaia data more directly. Like RV, this only works with individual draws because the Gaia points are "detrended" in a sense from the values of a particular draw (they move around the plot depending on the draw).
# Picks MAP sample
# gaiastarplot(model, chain)
# Pick specific arbitrary sample
idx = rand(1:size(chain,1)) # pick an integer randomly
Octofitter.gaiastarplot(model, chain, idx)
A good practice is to plot a few different values from the posterior, or e.g. plot draws from 5th, 50th, and 95th percentile in a key orbit parameter like
- semi-major axis
- eccentricity
- inclination
Here, we show semi-major axis/period
fig = Figure(size=(920,345))
percentile_positions = round.(Int, [0.05, 0.50, 0.95] .* size(chain,1))
indices = [partialsortperm(chain["b_a"][:], k) for k in percentile_positions]
# hint! Try `"b_e", "b_a", and "b_i"
ax1 = Octofitter.gaiastarplot!(fig[1,1], model, chain, indices[1])
ax2 = Octofitter.gaiastarplot!(fig[1,2], model, chain, indices[2])
ax3 = Octofitter.gaiastarplot!(fig[1,3], model, chain, indices[3])
ax1.title = "5th percentile period"
ax2.title = "50th percentile period"
ax3.title = "95th percentile period"
hideydecorations!(ax2)
hideydecorations!(ax3)
linkaxes!(ax1,ax2,ax3)
figWe can also plot a sky track, reconstructed against the official gaia solution (queried automatically when constructing the observation object).
# specific sample
idx = 42
# optional parameter keplerian_mult: exagerates the keplerian component for visualization
Octofitter.skytrackplot(model, chain, idx, keplerian_mult=10)
And of course, you can make a corner plot as usual:
using PairPlots
octocorner(model, chain, small=true)
Cross-validation
You can use the full suite of tools for construting models that subset different amounts of data in different ways. See "Cross Validataion".
The most powerful is exhaustive leave-one-out cross validataion plus calculation of expected log pointwise density (ELPD) to score the fit.
likelihood_mat, epochs = Octofitter.pointwise_like(model, chain)
# `likelihood_mat` is now a N_sample x N_data matrix.
using ParetoSmooth
result = psis_loo(
collect(likelihood_mat'),
chain_index=ones(Int,size(chain,1))
)Simulate data from a posterior draw, and re-fit with or without noise
Optional consistency checks–-could be used in a loop as part of e.g. simulation based calibration.
template = Octofitter.mcmcchain2result(model,chain,1)
sim_system = Octofitter.generate_from_params(model.system, template; add_noise=true)
sim_model = Octofitter.LogDensityModel(sim_system)
# Optional initialization speed up hack:
sim_model.starting_points = model.starting_points;
# then re-fit...
sim_chain, pt = octofit_pigeons(sim_model, n_rounds=5)
Octofitter.gaiastarplot(sim_model, sim_chain)Gaia BH 3
The following tutorial reproduces the fit to Gaia BH3. This one can take longer to run since the orbit is ultra well constrained. For that reason, we don't run it automatically when building the docs. Please go ahead and run the code on your own computer; ETA=approx 20 minutes.
As a first step, we will load the astrometry data for Gaia-BH3 and plot it:
headers = [
:transit_id
:ccd_id
:obs_time_tcb
:centroid_pos_al
:centroid_pos_error_al
:parallax_factor_al
:scan_pos_angle
:outlier_flag
]
df = CSV.read(joinpath(@__DIR__, "astrom.dat"), DataFrame, skipto=7, header=headers, delim=' ', ignorerepeated=true)
df.epoch = jd2mjd.(df.obs_time_tcb)
scatter(
df.obs_time_tcb,
df.centroid_pos_al,
)We can now construct a likelihood object for this data. We must also supply the Gaia ID, which will be used to query the full Gaia solution for this object (for now, using DR3):
ref_epoch_mjd = 57936.375
gaia_bh3_astrom_obs = GaiaDR4Astrom(
df,
gaia_id=4373465352415301632,
variables=@variables begin
astrometric_jitter ~ LogUniform(0.00001, 10) # mas
ra_offset_mas ~ Normal(0, 10000)
dec_offset_mas ~ Normal(0, 10000)
pmra ~ Uniform(-1000, 1000) # mas/yr
pmdec ~ Uniform(-1000, 1000) # mas/yr
# ra_offset_mas = -0.00154
# dec_offset_mas = 0.0019354
# pmra = 5.421
# pmdec = -24.121
plx = system.plx
ref_epoch = $ref_epoch_mjd
end
)This object also has published RV data from Gaia, which we can load and use as normal:
using CSV, DataFrames
using OctofitterRadialVelocity
headers_rv = [
:transit_id,
:obs_time_tcb,
:radial_velocity_kms,
:radial_velocity_err_kms,
]
dfrv = CSV.read(joinpath(@__DIR__, "epochrv.dat"), DataFrame, skipto=7, header=headers_rv, delim=' ', ignorerepeated=true)
dfrv.epoch = jd2mjd.(dfrv.obs_time_tcb)
dfrv.rv = dfrv.radial_velocity_kms * 1e3
dfrv.σ_rv = dfrv.radial_velocity_err_kms * 1e3
# Calculate mean RV for the prior
mean_rv = mean(dfrv.rv)
rvlike = StarAbsoluteRVObs(
dfrv,
name="GaiaRV",
variables=@variables begin
offset ~ Normal(mean_rv, 10_000) # wide prior on RV offset centred on mean RV
jitter ~ LogUniform(0.01, 100_000) # RV jitter parameter
end
)
errorbars(
dfrv.obs_time_tcb,
dfrv.rv,
dfrv.σ_rv
)Now, we define a model that incorporates this data:
mjup2msol = Octofitter.mjup2msol
ref_epoch_mjd = Octofitter.meta_gaia_DR3.ref_epoch_mjd
orbit_ref_epoch = mean(gaia_dr4_obs.table.epoch)
b = Planet(
name="BH",
basis=Visual{KepOrbit},
observations=[],
variables=@variables begin
a ~ Uniform(0, 1000)
e ~ Uniform(0, 0.99)
ω ~ Uniform(0,2pi)
i ~ Sine()
Ω ~ Uniform(0,2pi)
θ ~ Uniform(0,2pi)
tp = θ_at_epoch_to_tperi(θ, $orbit_ref_epoch; M=system.M, e, a, i, ω, Ω)
mass = system.M_sec / mjup2msol
end
)
# DR3 catalog position and reference epoch
ref_epoch_mjd = Octofitter.meta_gaia_DR3.ref_epoch_mjd
sys = System(
name="gaiadr4test",
companions=[b],
observations=[gaia_bh3_astrom_obs, rvlike,],
variables=@variables begin
M_pri ~ truncated(Normal(0.76,0.05),lower=0.1) # M sun
M_sec ~ LogUniform(1, 1000) # M sun
M = M_pri + M_sec
# Note: keep these physically plausible to prevent numerical errors
plx ~ Uniform(0.01,100) # mas
end
)
model = Octofitter.LogDensityModel(sys, verbosity=4)We will initialize the model starting positions and visualize them:
# Note: you can see the required format for paramter initialization by running:
# nt = Octofitter.drawfrompriors(model.system);
# println(nt)
# then remove any derived parameters (parameters in your model that are on the right of an `=`)
init_chain = initialize!(model, (;
M_pri = 0.7792923132247755,
M_sec = 36.032664849109906,
plx = 1.6686144513164856,
observations = (GaiaDR4 = (astrometric_jitter = 0.027554101045898238,),
GaiaRV = (offset = -359481.6706770764,jitter = 2143.4793485877644)),
planets = (BH = (
a = 18.905647598089196,
e = 0.7583328001601555,
i = 1.9216027029499319,
),
)))
octoplot(model, init_chain, show_rv=true)If you don't pick the starting point, you cabn also just run Pigeons for 8-10 rounds, which is recommended anyways for convergence, and the sampler will find this result.
Now, we can perform the fit. It is a little slow since we have many hundreds of RV and astrometry data points.
using Pigeons
chain, pt = octofit_pigeons(model, n_rounds=6) # might need more rounds to convergeincrement_n_rounds!(pt,1)
chain,pt = octofit_pigeons(pt)Finally, we can visualize the results:
octoplot(model, chain, show_rv=true, mark_epochs_mjd=mjd.([
"2017"
"2022"
"2027"
]))# Picks MAP sample
Octofitter.gaiastarplot(model, chain)octocorner(model, chain, small=true)Octofitter.rvpostplot(model, chain)