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 = -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     = 13.55 seconds
Compute duration  = 13.55 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.0019   2453.0346   2995.8774    1. ⋯
         plx      49.9997      0.0202     0.0003   4870.5825   2711.9755    1. ⋯
         b_e       0.1286      0.0900     0.0037    566.2207   1178.1371    1. ⋯
         b_i       0.6607      0.0927     0.0024   1519.9631   1770.4303    1. ⋯
        b_ωy       0.4232      0.6741     0.1829     20.6592     53.0489    1. ⋯
        b_ωx       0.3091      0.5382     0.1377     19.2080     62.6672    1. ⋯
        b_Ωy      -0.2441      0.4951     0.1164     19.8263     70.5232    1. ⋯
        b_Ωx      -0.4760      0.7021     0.1999     20.5009     55.0935    1. ⋯
         b_P      40.7825      5.6070     0.1564   1179.2672   1379.9865    1. ⋯
        b_θy       0.0743      0.0100     0.0002   4321.3781   3182.4946    1. ⋯
        b_θx      -1.0001      0.0978     0.0017   3450.2107   2677.3433    1. ⋯
         b_ω       0.0300      1.3574     0.3679     22.5583     55.2060    1. ⋯
         b_Ω      -1.3974      1.3439     0.3896     18.4444     58.2862    1. ⋯
         b_a      12.4075      1.1468     0.0313   1328.6982   2083.2826    1. ⋯
         b_θ      -1.4966      0.0072     0.0001   5037.5313   3852.1382    1. ⋯
        b_tp   43108.7541   7240.2071   118.4670   3290.8692   2474.0076    1. ⋯
                                                               2 columns omitted

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

           M       0.9724       1.0955       1.1601       1.2225       1.3513
         plx      49.9596      49.9862      49.9995      50.0136      50.0386
         b_e       0.0042       0.0504       0.1133       0.2014       0.3046
         b_i       0.4643       0.6027       0.6653       0.7272       0.8253
        b_ωy      -1.0291       0.1571       0.7334       0.9254       1.1147
        b_ωx      -0.9542       0.0700       0.3866       0.7234       1.0612
        b_Ωy      -1.0144      -0.6328      -0.2405      -0.0095       0.9110
        b_Ωx      -1.1500      -0.9660      -0.8076      -0.3443       1.0406
         b_P      30.4627      36.2813      41.0368      45.5837      49.5555
        b_θy       0.0561       0.0671       0.0739       0.0809       0.0946
        b_θx      -1.2042      -1.0611      -0.9952      -0.9337      -0.8209
         b_ω      -2.9368       0.0862       0.4020       0.8428       1.5691
         b_Ω      -2.8248      -2.2498      -1.8115      -1.5597       1.5063
         b_a      10.3557      11.4799      12.4586      13.3324      14.3907
         b_θ      -1.5109      -1.5015      -1.4966      -1.4918      -1.4825
        b_tp   32745.2750   35887.8836   48851.0975   50210.8197   50401.6567

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