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),
)
obs_pri = ObsPriorAstromONeil2019(astrom_like)
@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 obs_pri

@system TutoriaPrime begin
    M ~ truncated(Normal(1.2, 0.1), lower=0.1)
    plx ~ truncated(Normal(50.0, 0.02), lower=0.1)
end b

model = Octofitter.LogDensityModel(TutoriaPrime)
LogDensityModel for System TutoriaPrime of dimension 11 and 12 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)
Example block output

Now run the fit:

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     = 6.87 seconds
Compute duration  = 6.87 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.2165      0.0962     0.0127     64.9982     37.7014    1. ⋯
         plx      49.9998      0.0190     0.0004   2753.4143   1719.2433    1. ⋯
         b_e       0.0984      0.0687     0.0054    188.0304   1343.7429    1. ⋯
         b_i       0.4805      0.1587     0.0203     57.7191     41.7867    1. ⋯
        b_ωy      -0.0761      0.6668     0.0638    121.5000   1711.1788    1. ⋯
        b_ωx      -0.0444      0.7529     0.0690    157.0360   1642.1563    1. ⋯
        b_Ωy      -0.0983      0.7883     0.0890     88.1294    453.4278    1. ⋯
        b_Ωx       0.0044      0.6200     0.0464    205.8828   1048.1021    1. ⋯
         b_P      29.9620      5.1809     0.6138     88.7165    129.7613    1. ⋯
        b_θy       0.0734      0.0090     0.0004    455.9070   3052.9545    1. ⋯
        b_θx      -1.0115      0.1008     0.0070    224.2407   3203.2068    1. ⋯
         b_ω      -0.4435      1.8054     0.1927     81.1189    115.0526    1. ⋯
         b_Ω      -0.6215      1.9029     0.1311    234.6825     67.3985    1. ⋯
         b_a      10.2603      1.2519     0.1799     66.6416    107.1422    1. ⋯
         b_θ      -1.4984      0.0050     0.0001   3589.2378   3073.9981    1. ⋯
        b_tp   44297.8994   2043.3335   307.3211     62.4712     48.5024    1. ⋯
                                                               2 columns omitted

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

           M       1.0428       1.1474       1.2082       1.2794       1.4226
         plx      49.9623      49.9864      49.9999      50.0125      50.0374
         b_e       0.0037       0.0410       0.0883       0.1470       0.2446
         b_i       0.1614       0.3755       0.4816       0.5892       0.7642
        b_ωy      -1.0333      -0.7100      -0.1263       0.5312       1.0323
        b_ωx      -1.0824      -0.7659      -0.1822       0.7550       1.0689
        b_Ωy      -1.1004      -0.8805      -0.2369       0.7657       1.0697
        b_Ωx      -1.0318      -0.4995      -0.0236       0.5162       1.0353
         b_P      21.0802      26.2573      29.1970      33.3555      41.0000
        b_θy       0.0567       0.0670       0.0732       0.0795       0.0913
        b_θx      -1.2067      -1.0859      -1.0079      -0.9392      -0.8269
         b_ω      -2.9610      -2.1604      -0.5457       1.1293       2.8207
         b_Ω      -3.0853      -2.5609      -0.1559       0.7714       3.0682
         b_a       8.2471       9.3838      10.0032      10.9854      13.0020
         b_θ      -1.5083      -1.5017      -1.4984      -1.4952      -1.4884
        b_tp   39208.1609   43184.5734   44703.4934   45747.6023   47332.4163

Plot the MCMC results:

octoplot(model, results_obspri)
Example block output

Compare this with the previous fit using uniform priors:

octoplot(model_with_uniform_priors, 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