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     = 21.38 seconds
Compute duration  = 21.38 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.1606       0.0970     0.0020   2294.6372   2848.8208    1 ⋯
         plx      50.0001       0.0201     0.0003   4823.7848   2174.5579    0 ⋯
         b_e       0.2719       0.1336     0.0049    606.0626   1141.7799    1 ⋯
         b_i       0.7308       0.1471     0.0098    238.3376   1031.4444    1 ⋯
        b_ωy       0.0878       0.6558     0.1764     14.2274     94.9791    1 ⋯
        b_ωx      -0.0434       0.7604     0.2188     13.4459    114.9606    1 ⋯
        b_Ωy       0.0462       0.7586     0.2128     14.1056    112.3098    1 ⋯
        b_Ωx      -0.0700       0.6579     0.1916     12.7974     96.7960    2 ⋯
         b_P      70.6177      29.2478     1.7249    302.6526    477.0053    1 ⋯
        b_θy       0.0712       0.0087     0.0001   3606.8014   2989.1630    1 ⋯
        b_θx      -1.0017       0.0970     0.0016   3616.2053   3648.9332    1 ⋯
         b_ω      -0.5316       1.6113     0.4619     14.8380     88.4611    1 ⋯
         b_Ω      -0.8339       1.5986     0.4748     13.7254    127.0137    1 ⋯
         b_a      17.6132       4.8717     0.2920    295.2035    466.3562    1 ⋯
         b_θ      -1.4999       0.0053     0.0001   2657.1223   3273.5457    1 ⋯
        b_tp   38695.1478   14444.3998   376.6884   2086.3718    875.8126    1 ⋯
                                                               2 columns omitted

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

           M      0.9748       1.0954       1.1581       1.2253       1.3551
         plx     49.9611      49.9861      49.9999      50.0137      50.0399
         b_e      0.0215       0.1651       0.2862       0.3814       0.4845
         b_i      0.4361       0.6280       0.7333       0.8493       0.9781
        b_ωy     -1.0624      -0.2930       0.0582       0.7201       1.1068
        b_ωx     -1.1234      -0.8761       0.0220       0.7385       1.0914
        b_Ωy     -1.0821      -0.7979       0.0193       0.8727       1.0989
        b_Ωx     -1.0991      -0.6155      -0.1442       0.3423       1.0819
         b_P     32.1151      47.0704      64.5677      88.8796     137.8894
        b_θy      0.0557       0.0652       0.0704       0.0768       0.0898
        b_θx     -1.1976      -1.0654      -0.9979      -0.9342      -0.8225
         b_ω     -3.0357      -1.8275       0.0317       1.0067       1.7829
         b_Ω     -2.9514      -2.5503      -1.2518       0.3548       1.6373
         b_a     10.5737      13.6390      16.8631      20.8775      27.9471
         b_θ     -1.5100      -1.5035      -1.4999      -1.4962      -1.4893
        b_tp   3832.3395   29267.2421   49429.1785   50207.1312   50403.2402

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