Observable-Based Priors
This tutorial shows how to fit an orbit to relative astrometry using the observable-based priors ofO'Neil et al. 2019. Please cite that paper if you use this functionality.
We will fit the same astrometry as in the previous tutorial, and just change our priors.
using Octofitter
using CairoMakie
using PairPlots
using Distributions
astrom_dat = Table(;
epoch = [50000, 50120, 50240, 50360, 50480, 50600, 50720, 50840],
ra = [-494.4, -495.0, -493.7, -490.4, -485.2, -478.1, -469.1, -458.3],
dec = [-76.7, -44.9, -12.9, 19.1, 51.0, 82.8, 114.3, 145.3],
σ_ra = [12.6, 10.4, 9.9, 8.7, 8.0, 6.9, 5.8, 4.2],
σ_dec = [12.6, 10.4, 9.9, 8.7, 8.0, 6.9, 5.8, 4.2],
cor = [0.2, 0.5, 0.1, -0.8, 0.3, -0.0, 0.1, -0.2]
)
astrom_like = PlanetRelAstromLikelihood(
astrom_dat,
name = "obs_prior_example",
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
)
# We wrap the likelihood in this prior
obs_pri_astrom_like = ObsPriorAstromONeil2019(astrom_like)
planet_b = Planet(
name="b",
basis=Visual{KepOrbit},
# NOTE! We only provide the wrapped obs_pri_astrom_like
likelihoods=[obs_pri_astrom_like],
variables=@variables begin
M = system.M
e ~ Uniform(0.0, 0.5)
i ~ Sine()
ω ~ UniformCircular()
Ω ~ UniformCircular()
# Results will be sensitive to the prior on period
P ~ LogUniform(0.1, 150) # Period, years
a = ∛(M * P^2)
θ_x ~ Normal()
θ_y ~ Normal()
θ = atan(θ_y, θ_x)
tp = θ_at_epoch_to_tperi(θ, 50420; M, e, a, i, ω, Ω)
end
)
sys = System(
name="TutoriaPrime",
companions=[planet_b],
likelihoods=[],
variables=@variables begin
M ~ truncated(Normal(1.2, 0.1), lower=0.1)
plx ~ truncated(Normal(50.0, 0.02), lower=0.1)
end
)
model = Octofitter.LogDensityModel(sys)LogDensityModel for System TutoriaPrime of dimension 11 and 10 epochs with fields .ℓπcallback and .∇ℓπcallback
Initialize the model starting points and confirm the data are entered correctly:
init_chain = initialize!(model,)
octoplot(model, init_chain)
Now run the fit:
results_obspri = octofit(model,iterations=5000,)Chains MCMC chain (5000×33×1 Array{Float64, 3}):
Iterations = 1:1:5000
Number of chains = 1
Samples per chain = 5000
Wall duration = 11.9 seconds
Compute duration = 11.9 seconds
parameters = M, plx, b_e, b_i, b_ωx, b_ωy, b_Ωx, b_Ωy, b_P, b_θ_x, b_θ_y, b_ω, b_Ω, b_M, b_a, b_θ, b_tp, b_obspri_obs_prior_example_jitter, b_obspri_obs_prior_example_northangle, b_obspri_obs_prior_example_platescale
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 ⋯
Symbol Float64 Float64 Float64 ⋯
M 1.2064 0.0891 0.0020 1 ⋯
plx 50.0004 0.0199 0.0003 4 ⋯
b_e 0.0957 0.0692 0.0019 1 ⋯
b_i 0.4570 0.1535 0.0063 ⋯
b_ωx 0.0613 0.6603 0.0494 ⋯
b_ωy 0.1014 0.7554 0.0561 ⋯
b_Ωx 0.1013 0.7525 0.0654 ⋯
b_Ωy 0.0521 0.6624 0.0497 ⋯
b_P 29.1587 4.9520 0.1899 ⋯
b_θ_x 0.0911 0.0476 0.0014 1 ⋯
b_θ_y -1.2540 0.6465 0.0183 1 ⋯
b_ω -0.0678 1.7089 0.1125 ⋯
b_Ω -0.3908 1.6866 0.1223 ⋯
b_M 1.2064 0.0891 0.0020 1 ⋯
b_a 10.0458 1.1504 0.0455 ⋯
b_θ -1.4983 0.0055 0.0001 1 ⋯
b_tp 44591.0406 1814.9164 78.4797 ⋯
⋮ ⋮ ⋮ ⋮ ⋱
4 columns and 3 rows omitted
Quantiles
parameters 2.5% 25.0% 50.0% ⋯
Symbol Float64 Float64 Float64 ⋯
M 1.0421 1.1429 1.2035 ⋯
plx 49.9616 49.9872 50.0007 ⋯
b_e 0.0030 0.0383 0.0842 ⋯
b_i 0.1323 0.3582 0.4629 ⋯
b_ωx -1.0350 -0.5307 0.1377 ⋯
b_ωy -1.0698 -0.6988 0.2731 ⋯
b_Ωx -1.0712 -0.6957 0.2258 ⋯
b_Ωy -1.0297 -0.5413 0.1538 ⋯
b_P 20.8266 25.8350 28.5190 ⋯
b_θ_x 0.0162 0.0567 0.0864 ⋯
b_θ_y -2.7158 -1.6508 -1.1747 ⋯
b_ω -2.9342 -1.7930 0.4218 ⋯
b_Ω -2.9684 -2.1285 0.2000 ⋯
b_M 1.0421 1.1429 1.2035 ⋯
b_a 8.2016 9.2865 9.8542 ⋯
b_θ -1.5089 -1.5019 -1.4985 ⋯
b_tp 40021.2534 43659.2873 44892.4481 ⋯
⋮ ⋮ ⋮ ⋮ ⋱
2 columns and 3 rows omitted
Plot the MCMC results:
octoplot(model, results_obspri)
Compare this with the previous fit using uniform priors:
octoplot(model_with_uniform_priors, results_unif_pri)
We can compare the results in a corner plot:
octocorner(model,results_unif_pri,results_obspri,small=true)