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 = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 50120, ra = -502.570356287689, dec = -37.47217527025044, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 50240, ra = -498.2089148883798, dec = -7.927548139010479, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 50360, ra = -492.67768482682357, dec = 21.63557115669823, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 50480, ra = -485.9770335870402, dec = 51.147204404903704, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 50600, ra = -478.1095526888573, dec = 80.53589069730698, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 50720, ra = -469.0801731788123, dec = 109.72870493064629, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 50840, ra = -458.89628893460525, dec = 138.65128697876773, σ_ra = 10, σ_dec = 10, cor=0),
)
@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, 50)
    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     = 21.92 seconds
Compute duration  = 21.92 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   Flo ⋯

           M       1.1600      0.0954     0.0020   2238.7514   2672.7462    1. ⋯
         plx      50.0001      0.0193     0.0003   5527.9851   3108.6697    0. ⋯
         b_e       0.1338      0.0915     0.0034    707.3535   1245.2915    0. ⋯
         b_i       0.6629      0.0941     0.0021   2129.7356   2162.1189    1. ⋯
        b_ωy      -0.1303      0.7867     0.2195     16.5512    153.1249    1. ⋯
        b_ωx      -0.0820      0.6137     0.1593     17.1335    173.3665    1. ⋯
        b_Ωy       0.0472      0.5354     0.1298     17.4381    200.2400    1. ⋯
        b_Ωx       0.1368      0.8420     0.2416     16.2992    147.2926    1. ⋯
         b_P      40.9328      5.6871     0.1420   1474.7744   1690.7853    1. ⋯
        b_θy       0.0744      0.0103     0.0001   5226.3276   3055.0762    1. ⋯
        b_θx      -1.0009      0.1000     0.0014   5207.6429   3161.7118    0. ⋯
         b_ω      -1.0684      1.6993     0.4452     20.5543    389.8454    1. ⋯
         b_Ω      -0.2269      1.6288     0.4653     17.0463    169.5432    1. ⋯
         b_a      12.4369      1.1576     0.0283   1579.9510   1661.8649    1. ⋯
         b_θ      -1.4965      0.0071     0.0001   5309.3387   3794.2545    1. ⋯
        b_tp   42907.5306   7271.5183   111.8130   4077.1741   2882.0251    1. ⋯
                                                               2 columns omitted

Quantiles
  parameters         2.5%        25.0%        50.0%        75.0%        97.5%
      Symbol      Float64      Float64      Float64      Float64      Float64

           M       0.9800       1.0940       1.1592       1.2244       1.3496
         plx      49.9621      49.9875      50.0002      50.0133      50.0382
         b_e       0.0040       0.0523       0.1218       0.2076       0.3056
         b_i       0.4700       0.6034       0.6678       0.7296       0.8273
        b_ωy      -1.0957      -0.8770      -0.3582       0.7517       1.0742
        b_ωx      -1.0299      -0.5693      -0.1992       0.4203       1.0274
        b_Ωy      -0.9739      -0.2866       0.0772       0.4011       0.9629
        b_Ωx      -1.0978      -0.8252       0.4972       0.9226       1.1231
         b_P      30.2860      36.3524      41.2991      45.9966      49.5858
        b_θy       0.0557       0.0673       0.0739       0.0810       0.0968
        b_θx      -1.2114      -1.0658      -0.9959      -0.9307      -0.8223
         b_ω      -3.0253      -2.6686      -1.8640       0.4768       1.6742
         b_Ω      -2.7594      -1.8434       0.5352       1.3076       1.6099
         b_a      10.3289      11.5001      12.4923      13.3734      14.3853
         b_θ      -1.5104      -1.5014      -1.4965      -1.4916      -1.4830
        b_tp   32729.3018   35795.2007   40519.1799   50199.5876   50399.7334
octoplot(model, results_obspri)
Example block output

Compare this with the previous fit using uniform priors:

octoplot(model, results_unif_pri)
Example block output

We can compare the results in a corner plot:

octocorner(model,results_unif_pri,results_obspri,small=true)
Example block output