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

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     = 8.3 seconds
Compute duration  = 8.3 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      r ⋯
      Symbol      Float64     Float64   Float64     Float64     Float64   Floa ⋯

           M       1.2033      0.0893    0.0020   1939.3017   2773.0591    1.0 ⋯
         plx      49.9998      0.0199    0.0003   4120.8259   2610.2398    1.0 ⋯
         b_e       0.0989      0.0694    0.0023    850.5128    787.4438    1.0 ⋯
         b_i       0.4507      0.1500    0.0061    605.8869    731.1533    1.0 ⋯
        b_ωy      -0.0761      0.6587    0.0373    342.9196   1388.2918    1.0 ⋯
        b_ωx      -0.0929      0.7541    0.0467    325.8367   2292.8517    1.0 ⋯
        b_Ωy      -0.1251      0.7502    0.0541    227.6519   1472.3253    1.0 ⋯
        b_Ωx      -0.0791      0.6580    0.0391    332.8647   1804.9848    1.0 ⋯
         b_P      29.0598      4.8976    0.2135    593.4195    333.2965    1.0 ⋯
        b_θy       0.0725      0.0088    0.0002   2988.4600   2887.4696    1.0 ⋯
        b_θx      -1.0001      0.0978    0.0018   3143.0895   2953.8144    1.0 ⋯
         b_ω      -0.4421      1.8026    0.1024    406.6162    706.6098    1.0 ⋯
         b_Ω      -0.6659      1.8637    0.0996    488.1653    637.6421    1.0 ⋯
         b_a      10.0121      1.1162    0.0490    574.3237    546.2951    1.0 ⋯
         b_θ      -1.4984      0.0052    0.0001   1952.1540   1272.8857    1.0 ⋯
        b_tp   44600.1600   1884.3469   93.7578    637.1098    319.7142    1.0 ⋯
                                                               2 columns omitted

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

           M       1.0361       1.1432       1.2000       1.2613       1.3852
         plx      49.9602      49.9866      50.0000      50.0133      50.0395
         b_e       0.0036       0.0403       0.0884       0.1468       0.2498
         b_i       0.1236       0.3555       0.4585       0.5548       0.7268
        b_ωy      -1.0632      -0.6763      -0.1228       0.5075       1.0362
        b_ωx      -1.0823      -0.8248      -0.2390       0.7110       1.0667
        b_Ωy      -1.0801      -0.8433      -0.2827       0.6650       1.0619
        b_Ωx      -1.0546      -0.6636      -0.1623       0.4851       1.0396
         b_P      20.8772      25.7885      28.4909      31.7857      40.4012
        b_θy       0.0567       0.0663       0.0721       0.0782       0.0913
        b_θx      -1.2065      -1.0635      -0.9976      -0.9319      -0.8230
         b_ω      -2.9623      -2.0674      -0.7547       1.1327       2.9160
         b_Ω      -3.0323      -2.4566      -1.0171       0.7854       2.9961
         b_a       8.1993       9.2610       9.8381      10.6161      12.6354
         b_θ      -1.5087      -1.5019      -1.4984      -1.4948      -1.4883
        b_tp   39669.4635   43700.1849   44955.3443   45899.7784   47374.1013

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