Hipparcos Modelling
This tutorial explains how to model Hipparcos IAD data. The first example reproduces the catalog values of position, parallax, and proper motion. The second uses Hipparcos to constrain the mass of a directly imaged planet.
Reproduce Catalog Values
This is the so-called "Nielsen" test from Nielsen et al (2020) and available in Orbitize!.
We start by using a system with a planet with zero mass to fit the straight line motion.
using Octofitter
using Distributions
using CairoMakie
hip_like = Octofitter.HipparcosIADLikelihood(
hip_id=21547,
renormalize=true, # default: true
variables=@variables begin
# Optional: flux ratio for luminous companions, one entry per companion
# fluxratio ~ Product([Uniform(0, 1), Uniform(0, 1)]) # uncomment if needed for unresolved companions
end
)
planet_b = Planet(
name="b",
basis=AbsoluteVisual{KepOrbit},
variables=@variables begin
mass = 0.
e = 0.
ω = 0.
a = 1.
i = 0
Ω = 0.
tp = 0.
end
)
sys = System(
name="c_Eri_straight_line",
companions=[planet_b],
likelihoods=[hip_like],
variables=@variables begin
M = 1.0 # Host mass not important for this example
rv = 0.0 # system RV not significant for this example
plx ~ Uniform(10,100)
pmra ~ Uniform(-100, 100)
pmdec ~ Uniform(-100, 100)
# It is convenient to put a prior of the catalog value +- 10,000 mas on position
ra_hip_offset_mas ~ Normal(0, 10000)
dec_hip_offset_mas ~ Normal(0, 10000)
dec = $hip_like.hip_sol.dedeg + ra_hip_offset_mas/60/60/1000
ra = $hip_like.hip_sol.radeg + dec_hip_offset_mas/60/60/1000/cosd(dec)
ref_epoch = Octofitter.hipparcos_catalog_epoch_mjd
end
)
model = Octofitter.LogDensityModel(sys)LogDensityModel for System c_Eri_straight_line of dimension 5 and 104 epochs with fields .ℓπcallback and .∇ℓπcallback
Let's initialize the starting point for the chains to reasonable values
initialize!(model, (;
plx=34.,
pmra=44.25,
pmdec=-64.5,
ra_hip_offset_mas=0.,
dec_hip_offset_mas=0.,
))Chains MCMC chain (1000×18×1 Array{Float64, 3}):
Iterations = 1:1:1000
Number of chains = 1
Samples per chain = 1000
Wall duration = 3.28 seconds
Compute duration = 3.28 seconds
parameters = plx, pmra, pmdec, ra_hip_offset_mas, dec_hip_offset_mas, M, rv, dec, ra, ref_epoch, b_mass, b_e, b_ω, b_a, b_i, b_Ω, b_tp
internals = logpost
Summary Statistics
parameters mean std mcse ess_bulk ess_tail ⋯
Symbol Float64 Float64 Float64 Float64 Float64 ⋯
plx 34.0000 0.0000 NaN NaN NaN ⋯
pmra 44.2500 0.0000 NaN NaN NaN ⋯
pmdec -64.5000 0.0000 NaN NaN NaN ⋯
ra_hip_offset_mas 0.0000 0.0000 NaN NaN NaN ⋯
dec_hip_offset_mas 0.0000 0.0000 NaN NaN NaN ⋯
M 1.0000 0.0000 NaN NaN NaN ⋯
rv 0.0000 0.0000 NaN NaN NaN ⋯
dec -2.4734 0.0000 0.0000 NaN NaN ⋯
ra 69.4004 0.0000 NaN NaN NaN ⋯
ref_epoch 48348.5625 0.0000 NaN NaN NaN ⋯
b_mass 0.0000 0.0000 NaN NaN NaN ⋯
b_e 0.0000 0.0000 NaN NaN NaN ⋯
b_ω 0.0000 0.0000 NaN NaN NaN ⋯
b_a 1.0000 0.0000 NaN NaN NaN ⋯
b_i 0.0000 0.0000 NaN NaN NaN ⋯
b_Ω 0.0000 0.0000 NaN NaN NaN ⋯
b_tp 0.0000 0.0000 NaN NaN NaN ⋯
2 columns omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% ⋯
Symbol Float64 Float64 Float64 Float64 ⋯
plx 34.0000 34.0000 34.0000 34.0000 ⋯
pmra 44.2500 44.2500 44.2500 44.2500 ⋯
pmdec -64.5000 -64.5000 -64.5000 -64.5000 - ⋯
ra_hip_offset_mas 0.0000 0.0000 0.0000 0.0000 ⋯
dec_hip_offset_mas 0.0000 0.0000 0.0000 0.0000 ⋯
M 1.0000 1.0000 1.0000 1.0000 ⋯
rv 0.0000 0.0000 0.0000 0.0000 ⋯
dec -2.4734 -2.4734 -2.4734 -2.4734 ⋯
ra 69.4004 69.4004 69.4004 69.4004 ⋯
ref_epoch 48348.5625 48348.5625 48348.5625 48348.5625 483 ⋯
b_mass 0.0000 0.0000 0.0000 0.0000 ⋯
b_e 0.0000 0.0000 0.0000 0.0000 ⋯
b_ω 0.0000 0.0000 0.0000 0.0000 ⋯
b_a 1.0000 1.0000 1.0000 1.0000 ⋯
b_i 0.0000 0.0000 0.0000 0.0000 ⋯
b_Ω 0.0000 0.0000 0.0000 0.0000 ⋯
b_tp 0.0000 0.0000 0.0000 0.0000 ⋯
1 column omitted
We can now sample from the model using Hamiltonian Monte Carlo. This should only take about 15 seconds.
using Pigeons
chain,pt = octofit_pigeons(model, n_rounds=6)(chain = MCMC chain (64×21×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))Plot the posterior values:
octoplot(model,chain,show_astrom=false,show_astrom_time=false)
We now visualize the model fit compared to the Hipparcos catalog values:
using LinearAlgebra, StatsBase
fig = Figure(size=(1080,720))
j = i = 1
for prop in (
(;chain=:ra, hip=:radeg, hip_err=:e_ra),
(;chain=:dec, hip=:dedeg, hip_err=:e_de),
(;chain=:plx, hip=:plx, hip_err=:e_plx),
(;chain=:pmra, hip=:pm_ra, hip_err=:e_pmra),
(;chain=:pmdec, hip=:pm_de, hip_err=:e_pmde)
)
global i, j, ax
ax = Axis(
fig[j,i],
xlabel=string(prop.chain),
)
i+=1
if i > 3
j+=1
i = 1
end
unc = hip_like.hip_sol[prop.hip_err]
if prop.chain == :ra
unc /= 60*60*1000 * cosd(hip_like.hip_sol.dedeg)
end
if prop.chain == :dec
unc /= 60*60*1000
end
if prop.hip == :zero
n = Normal(0, unc)
else
mu = hip_like.hip_sol[prop.hip]
n = Normal(mu, unc)
end
n0,n1=quantile.(n,(1e-4, 1-1e-4))
nxs = range(n0,n1,length=200)
h = fit(Histogram, chain[prop.chain][:], nbins=55)
h = normalize(h, mode=:pdf)
barplot!(ax, (h.edges[1][1:end-1] .+ h.edges[1][2:end])./2, h.weights, gap=0, color=:red, label="posterior")
lines!(ax, nxs, pdf.(n,nxs), label="Hipparcos Catalog", color=:black, linewidth=2)
end
Legend(fig[i-1,j+1],ax,tellwidth=false)
fig
Constrain Planet Mass
We now allow the planet to have a non zero mass and have free orbit. We start by specifying relative astrometry data on the planet, collated by Jason Wang and co. on whereistheplanet.com.
astrom_dat = Table(;
epoch = [57009.1, 57052.1, 57053.1, 57054.3, 57266.4, 57332.2, 57374.2, 57376.2, 57415.0, 57649.4, 57652.4, 57739.1, 58068.3, 58442.2],
sep = [454.24, 451.81, 456.8, 461.5, 455.1, 452.88, 455.91, 455.01, 454.46, 454.81, 451.43, 449.39, 447.54, 434.22],
σ_sep = [1.88, 2.06, 2.57, 23.9, 2.23, 5.41, 6.23, 3.03, 6.03, 2.02, 2.67, 2.15, 3.02, 2.01],
pa = [2.98835, 2.96723, 2.97038, 2.97404, 2.91994, 2.89934, 2.89131, 2.89184, 2.8962, 2.82394, 2.82272, 2.79357, 2.70927, 2.61171],
σ_pa = [0.00401426, 0.00453786, 0.00523599, 0.0523599, 0.00453786, 0.00994838, 0.00994838, 0.00750492, 0.00890118, 0.00453786, 0.00541052, 0.00471239, 0.00680678, 0.00401426]
)
astrom_like1 = PlanetRelAstromLikelihood(
astrom_dat,
name="VLT/SPHERE",
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
)PlanetRelAstromLikelihood Table with 5 columns and 14 rows:
epoch sep σ_sep pa σ_pa
┌────────────────────────────────────────────
1 │ 57009.1 454.24 1.88 2.98835 0.00401426
2 │ 57052.1 451.81 2.06 2.96723 0.00453786
3 │ 57053.1 456.8 2.57 2.97038 0.00523599
4 │ 57054.3 461.5 23.9 2.97404 0.0523599
5 │ 57266.4 455.1 2.23 2.91994 0.00453786
6 │ 57332.2 452.88 5.41 2.89934 0.00994838
7 │ 57374.2 455.91 6.23 2.89131 0.00994838
8 │ 57376.2 455.01 3.03 2.89184 0.00750492
9 │ 57415.0 454.46 6.03 2.8962 0.00890118
10 │ 57649.4 454.81 2.02 2.82394 0.00453786
11 │ 57652.4 451.43 2.67 2.82272 0.00541052
12 │ 57739.1 449.39 2.15 2.79357 0.00471239
13 │ 58068.3 447.54 3.02 2.70927 0.00680678
14 │ 58442.2 434.22 2.01 2.61171 0.00401426We specify our full model:
planet_b_mass = Planet(
name="b",
basis=AbsoluteVisual{KepOrbit},
likelihoods=[astrom_like1],
variables=@variables begin
a ~ truncated(Normal(10,1),lower=0.1)
e ~ Uniform(0,0.99)
ω ~ Uniform(0, 2pi)
i ~ Sine()
Ω ~ Uniform(0, 2pi)
θ ~ Uniform(0, 2pi)
M = system.M
tp = θ_at_epoch_to_tperi(θ, 58442.2; M, e, a, i, ω, Ω)
mass = system.M_sec
end
)
sys_mass = System(
name="cEri",
companions=[planet_b_mass],
likelihoods=[hip_like],
variables=@variables begin
M_pri ~ truncated(Normal(1.75,0.05), lower=0.03) # Msol
M_sec ~ LogUniform(0.1, 100) # MJup
M = M_pri + M_sec*Octofitter.mjup2msol # Msol
rv = 12.60e3 # m/s
plx ~ Uniform(20,40)
pmra ~ Uniform(-100, 100)
pmdec ~ Uniform(-100, 100)
# It is convenient to put a prior of the catalog value +- 1000 mas on position
ra_hip_offset_mas ~ Normal(0, 1000)
dec_hip_offset_mas ~ Normal(0, 1000)
dec = $hip_like.hip_sol.dedeg + ra_hip_offset_mas/60/60/1000
ra = $hip_like.hip_sol.radeg + dec_hip_offset_mas/60/60/1000/cos(dec)
ref_epoch = Octofitter.hipparcos_catalog_epoch_mjd
end
)
model = Octofitter.LogDensityModel(sys_mass)LogDensityModel for System cEri of dimension 13 and 118 epochs with fields .ℓπcallback and .∇ℓπcallback
Initialize the starting points, and confirm the data are entered correcly:
init_chain = initialize!(model, (;
plx=34.,
pmra=44.25,
pmdec=-64.5,
ra_hip_offset_mas=0.,
dec_hip_offset_mas=0.,
))
octoplot(model, init_chain)
Now we sample:
using Pigeons
chain,pt = octofit_pigeons(model, n_rounds=8, explorer=SliceSampler())
chainChains MCMC chain (256×28×1 Array{Float64, 3}):
Iterations = 1:1:256
Number of chains = 1
Samples per chain = 256
Wall duration = 481.04 seconds
Compute duration = 481.04 seconds
parameters = M_pri, M_sec, plx, pmra, pmdec, ra_hip_offset_mas, dec_hip_offset_mas, M, rv, dec, ra, ref_epoch, b_a, b_e, b_ω, b_i, b_Ω, b_θ, b_M, b_tp, b_mass, b_VLT_SPHERE_jitter, b_VLT_SPHERE_northangle, b_VLT_SPHERE_platescale
internals = loglike, logpost, logprior, pigeons_logpotential
Summary Statistics
parameters mean std mcse ess_bulk ess_tail ⋯
Symbol Float64 Float64 Float64 Float64 Float64 ⋯
M_pri 1.7475 0.0492 0.0057 78.1539 178.0162 ⋯
M_sec 26.9069 35.8649 6.6931 30.0818 122.0133 ⋯
plx 33.9560 0.3400 0.0237 205.3933 197.3019 ⋯
pmra 44.7649 0.8434 0.1366 46.0978 246.3629 ⋯
pmdec -64.2144 0.3935 0.0290 167.6616 195.1355 ⋯
ra_hip_offset_mas -5.9463 8.1096 1.5016 42.4042 151.2808 ⋯
dec_hip_offset_mas -1.4275 1.8722 0.2884 40.1104 139.0824 ⋯
M 1.7732 0.0664 0.0101 46.7210 200.2393 ⋯
rv 12600.0000 0.0000 NaN NaN NaN ⋯
dec -2.4734 0.0000 0.0000 42.4042 151.2808 ⋯
ra 69.4004 0.0000 0.0000 40.1104 139.0824 ⋯
ref_epoch 48348.5625 0.0000 NaN NaN NaN ⋯
b_a 10.6531 0.7028 0.2057 163.1165 19.2279 ⋯
b_e 0.5513 0.0671 0.0175 41.9941 12.4342 ⋯
b_ω 2.2648 1.6037 0.1071 211.5730 191.8005 ⋯
b_i 2.4046 0.0797 0.0164 93.1482 16.1362 ⋯
b_Ω 1.8002 1.8319 0.1336 14.7048 25.7725 ⋯
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱
2 columns and 7 rows omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% ⋯
Symbol Float64 Float64 Float64 Float64 ⋯
M_pri 1.6578 1.7112 1.7506 1.7810 ⋯
M_sec 0.1130 0.5094 2.7644 66.2746 ⋯
plx 33.3412 33.7111 33.9676 34.2013 ⋯
pmra 43.7333 44.1439 44.4435 45.6177 ⋯
pmdec -64.9039 -64.4992 -64.2299 -63.9306 ⋯
ra_hip_offset_mas -20.5706 -14.4452 -0.6219 -0.1432 ⋯
dec_hip_offset_mas -5.4899 -2.6184 -0.3490 0.0059 ⋯
M 1.6596 1.7281 1.7643 1.8248 ⋯
rv 12600.0000 12600.0000 12600.0000 12600.0000 12 ⋯
dec -2.4734 -2.4734 -2.4734 -2.4734 ⋯
ra 69.4004 69.4004 69.4004 69.4004 ⋯
ref_epoch 48348.5625 48348.5625 48348.5625 48348.5625 48 ⋯
b_a 9.4950 10.2052 10.5708 10.8262 ⋯
b_e 0.3595 0.5103 0.5650 0.5976 ⋯
b_ω 0.9023 1.1213 1.2919 4.2788 ⋯
b_i 2.2425 2.3693 2.4195 2.4516 ⋯
b_Ω 0.1078 0.3884 0.6469 3.6923 ⋯
⋮ ⋮ ⋮ ⋮ ⋮ ⋱
1 column and 7 rows omitted
octoplot(model, chain, show_mass=true)
We see that we constrained both the orbit and the parallax. The mass is not strongly constrained by Hipparcos.