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_like = PlanetRelAstromLikelihood(
(epoch = 50000, ra = -494.4, dec = -76.7, σ_ra = 12.6, σ_dec = 12.6, cor= 0.2),
(epoch = 50120, ra = -495.0, dec = -44.9, σ_ra = 10.4, σ_dec = 10.4, cor= 0.5),
(epoch = 50240, ra = -493.7, dec = -12.9, σ_ra = 9.9, σ_dec = 9.9, cor= 0.1),
(epoch = 50360, ra = -490.4, dec = 19.1, σ_ra = 8.7, σ_dec = 8.7, cor= -0.8),
(epoch = 50480, ra = -485.2, dec = 51.0, σ_ra = 8.0, σ_dec = 8.0, cor= 0.3),
(epoch = 50600, ra = -478.1, dec = 82.8, σ_ra = 6.9, σ_dec = 6.9, cor= -0.0),
(epoch = 50720, ra = -469.1, dec = 114.3, σ_ra = 5.8, σ_dec = 5.8, cor= 0.1),
(epoch = 50840, ra = -458.3, dec = 145.3, σ_ra = 4.2, σ_dec = 4.2, cor= -0.2),
)
@planet b Visual{KepOrbit} begin
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 = ∛(system.M * b.P^2)
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,b,50420)
end astrom_like ObsPriorAstromONeil2019(astrom_like)
@system TutoriaPrime begin
M ~ truncated(Normal(1.2, 0.1), lower=0)
plx ~ truncated(Normal(50.0, 0.02), lower=0)
end b
model = Octofitter.LogDensityModel(TutoriaPrime)
LogDensityModel for System TutoriaPrime of dimension 11 with fields .ℓπcallback and .∇ℓπcallback
Now run the fit:
using Random
Random.seed!(0)
results_obspri = octofit(model,iterations=5000,)
Chains MCMC chain (5000×29×1 Array{Float64, 3}):
Iterations = 1:1:5000
Number of chains = 1
Samples per chain = 5000
Wall duration = 48.26 seconds
Compute duration = 48.26 seconds
parameters = M, plx, b_e, b_i, b_ωy, b_ωx, b_Ωy, b_Ωx, b_P, b_θy, b_θx, b_ω, b_Ω, b_a, b_θ, b_tp
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 ess_bulk ess_tail ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Fl ⋯
M 1.1568 0.1022 0.0022 2111.2420 2525.0555 1 ⋯
plx 50.0000 0.0199 0.0003 5202.2572 3049.3321 1 ⋯
b_e 0.2826 0.1391 0.0038 1245.4620 1399.8625 1 ⋯
b_i 0.7251 0.1586 0.0107 227.8076 989.6467 1 ⋯
b_ωy -0.0528 0.6734 0.1627 17.8597 166.8063 1 ⋯
b_ωx -0.0303 0.7438 0.2007 15.2318 152.0828 1 ⋯
b_Ωy 0.0242 0.7421 0.1937 16.2568 129.3530 1 ⋯
b_Ωx 0.0532 0.6819 0.1808 16.0138 178.0855 1 ⋯
b_P 72.2710 31.4203 1.8334 322.5249 1108.1747 1 ⋯
b_θy 0.0712 0.0092 0.0002 2843.0765 2646.4868 1 ⋯
b_θx -1.0021 0.1018 0.0019 2864.5395 2525.0034 1 ⋯
b_ω -0.6826 1.7474 0.4694 17.9601 214.6823 1 ⋯
b_Ω -0.6535 1.7001 0.4851 15.1541 126.8492 1 ⋯
b_a 17.8389 5.2188 0.3078 316.2102 1030.7748 1 ⋯
b_θ -1.4999 0.0055 0.0001 2414.2848 3272.9043 1 ⋯
b_tp 38199.1498 14860.2848 370.4062 2147.8643 1793.3817 1 ⋯
2 columns omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 0.9608 1.0871 1.1554 1.2234 1.3620
plx 49.9609 49.9865 50.0002 50.0137 50.0384
b_e 0.0181 0.1702 0.3050 0.3991 0.4880
b_i 0.4103 0.6221 0.7303 0.8490 0.9893
b_ωy -1.0943 -0.6679 -0.0565 0.3913 1.0760
b_ωx -1.1062 -0.8331 -0.0640 0.7813 1.0986
b_Ωy -1.0918 -0.8093 0.0274 0.8358 1.0881
b_Ωx -1.0862 -0.3887 0.1882 0.5754 1.1182
b_P 31.2928 46.8876 65.1439 91.3860 141.8722
b_θy 0.0549 0.0647 0.0705 0.0770 0.0911
b_θx -1.2199 -1.0672 -0.9953 -0.9314 -0.8208
b_ω -3.0718 -2.1139 -1.4212 1.1156 1.9890
b_Ω -2.9651 -2.6137 0.1957 0.6631 1.6960
b_a 10.3759 13.6626 16.9471 21.2361 28.8588
b_θ -1.5107 -1.5036 -1.4999 -1.4962 -1.4889
b_tp 1171.4499 28684.6853 49341.0311 50200.2325 50399.2597
octoplot(model, results_obspri)
![Example block output](62c3c051.png)
Compare this with the previous fit using uniform priors:
octoplot(model, results_unif_pri)
![Example block output](7e484d47.png)
We can compare the results in a corner plot:
octocorner(model,results_unif_pri,results_obspri,small=true)
![Example block output](44b33f66.png)