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     = 10.47 seconds
Compute duration  = 10.47 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.2020      0.0884    0.0018   2303.2019   2758.6381    0.9 ⋯
         plx      49.9997      0.0199    0.0003   5168.0853   2581.8784    1.0 ⋯
         b_e       0.1006      0.0721    0.0023   1104.7399    772.1626    0.9 ⋯
         b_i       0.4524      0.1570    0.0055    798.5893   1046.9439    1.0 ⋯
        b_ωy       0.0358      0.6647    0.0373    350.4932   1234.5297    1.0 ⋯
        b_ωx       0.0528      0.7631    0.0540    261.7209   2088.5049    1.0 ⋯
        b_Ωy       0.0679      0.7723    0.0611    204.7625   1253.8436    1.0 ⋯
        b_Ωx       0.0539      0.6432    0.0405    308.2016   1784.5715    1.0 ⋯
         b_P      29.4420      5.4539    0.2487    683.9480    524.0619    1.0 ⋯
        b_θy       0.0727      0.0092    0.0002   2877.8512   2628.9731    0.9 ⋯
        b_θx      -1.0026      0.1016    0.0020   2574.0240   2321.0410    0.9 ⋯
         b_ω      -0.1024      1.7347    0.1044    401.6630   1462.0168    1.0 ⋯
         b_Ω      -0.3993      1.7427    0.1040    412.3514   1500.4011    1.0 ⋯
         b_a      10.0930      1.2385    0.0545    674.8492    495.9866    1.0 ⋯
         b_θ      -1.4984      0.0053    0.0001   3897.8116   3311.7056    0.9 ⋯
        b_tp   44575.6427   1932.8438   85.5781    765.9297    580.7568    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.0337       1.1397       1.2003       1.2622       1.3793
         plx      49.9605      49.9859      49.9996      50.0130      50.0384
         b_e       0.0036       0.0419       0.0882       0.1456       0.2662
         b_i       0.1290       0.3485       0.4581       0.5644       0.7454
        b_ωy      -1.0554      -0.5625       0.0782       0.6395       1.0596
        b_ωx      -1.0798      -0.7433       0.1552       0.8101       1.0899
        b_Ωy      -1.0753      -0.7666       0.1658       0.8382       1.0844
        b_Ωx      -1.0231      -0.5029       0.1370       0.6012       1.0457
         b_P      21.0011      25.8447      28.5169      32.2881      42.1592
        b_θy       0.0564       0.0662       0.0722       0.0785       0.0920
        b_θx      -1.2179      -1.0677      -0.9974      -0.9319      -0.8184
         b_ω      -2.9450      -1.8106       0.2769       1.2491       2.9135
         b_Ω      -3.0126      -2.2243       0.1963       0.8313       2.8990
         b_a       8.2144       9.2721       9.8425      10.7180      13.0621
         b_θ      -1.5088      -1.5020      -1.4984      -1.4948      -1.4881
        b_tp   40086.8519   43695.4157   44862.2418   45810.8068   47612.5242

octoplot(model, results_obspri)

Compare this with the previous fit using uniform priors:

octoplot(model_with_uniform_priors, results_unif_pri)

We can compare the results in a corner plot:

octocorner(model,results_unif_pri,results_obspri,small=true)