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 = 10.39 seconds
Compute duration = 10.39 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.2100 0.0865 0.0063 ⋯
plx 49.9993 0.0200 0.0006 ⋯
b_e 0.0998 0.0653 0.0025 ⋯
b_i 0.4868 0.1873 0.0355 ⋯
b_ωx 0.1293 0.6975 0.1036 ⋯
b_ωy 0.0361 0.7238 0.0550 ⋯
b_Ωx 0.1553 0.7900 0.0847 ⋯
b_Ωy 0.0102 0.6026 0.0395 ⋯
b_P 31.0566 7.4633 1.7707 ⋯
b_θ_x 0.0828 0.0478 0.0049 ⋯
b_θ_y -1.1405 0.6575 0.0688 ⋯
b_ω -0.1426 1.6474 0.1134 ⋯
b_Ω -0.4538 1.6402 0.1143 ⋯
b_M 1.2100 0.0865 0.0063 ⋯
b_a 10.4639 1.7058 0.4063 ⋯
b_θ -1.4981 0.0052 0.0001 ⋯
b_tp 43815.2276 2917.3232 723.1626 ⋯
⋮ ⋮ ⋮ ⋮ ⋱
4 columns and 3 rows omitted
Quantiles
parameters 2.5% 25.0% 50.0% ⋯
Symbol Float64 Float64 Float64 ⋯
M 1.0439 1.1463 1.2127 ⋯
plx 49.9572 49.9865 49.9999 ⋯
b_e 0.0070 0.0509 0.0874 ⋯
b_i 0.1160 0.3617 0.4837 ⋯
b_ωx -1.0326 -0.5211 0.2245 ⋯
b_ωy -1.0668 -0.6981 0.0436 ⋯
b_Ωx -1.0663 -0.7308 0.4058 ⋯
b_Ωy -1.0327 -0.4683 0.1054 ⋯
b_P 20.9571 25.8700 29.1032 ⋯
b_θ_x 0.0143 0.0438 0.0759 ⋯
b_θ_y -2.5957 -1.5727 -1.0506 ⋯
b_ω -2.8774 -1.7619 0.0494 ⋯
b_Ω -2.9893 -2.1438 0.1256 ⋯
b_M 1.0439 1.1463 1.2127 ⋯
b_a 8.2159 9.2827 9.9738 ⋯
b_θ -1.5084 -1.5015 -1.4979 ⋯
b_tp 36877.1742 43045.3937 44730.2800 ⋯
⋮ ⋮ ⋮ ⋮ ⋱
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)