GitHub

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

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