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.44 seconds
Compute duration  = 10.44 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.2051      0.0874    0.0019   2151.3258   2748.0350    0.9 ⋯
         plx      49.9999      0.0200    0.0003   4790.1753   2847.9926    1.0 ⋯
         b_e       0.0948      0.0684    0.0021    985.2859   1787.4257    1.0 ⋯
         b_i       0.4574      0.1556    0.0064    581.9840    722.3258    1.0 ⋯
        b_ωy       0.0189      0.6611    0.0355    410.7945   1997.5759    1.0 ⋯
        b_ωx       0.0293      0.7640    0.0511    274.9460   2376.8815    1.0 ⋯
        b_Ωy       0.0260      0.7703    0.0623    189.5369   1520.3264    1.0 ⋯
        b_Ωx       0.0094      0.6541    0.0546    175.2032   1045.3184    1.0 ⋯
         b_P      29.3175      5.1027    0.2111    604.3258   1050.0339    1.0 ⋯
        b_θy       0.0729      0.0091    0.0002   2771.6579   2726.9665    1.0 ⋯
        b_θx      -1.0053      0.1035    0.0020   2644.1614   2358.9609    0.9 ⋯
         b_ω      -0.1105      1.7465    0.0984    438.1638   1492.3421    1.0 ⋯
         b_Ω      -0.5067      1.7542    0.1271    263.1045   1206.4141    1.0 ⋯
         b_a      10.0758      1.1749    0.0500    567.1201   1112.4369    1.0 ⋯
         b_θ      -1.4983      0.0053    0.0001   3339.4804   2968.5909    1.0 ⋯
        b_tp   44752.0952   1738.0276   55.4514    974.2424    964.4621    0.9 ⋯
                                                               2 columns omitted

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

           M       1.0414       1.1450       1.2034       1.2637       1.3781
         plx      49.9597      49.9867      50.0003      50.0131      50.0391
         b_e       0.0042       0.0388       0.0830       0.1385       0.2492
         b_i       0.1348       0.3585       0.4665       0.5674       0.7386
        b_ωy      -1.0359      -0.5807       0.0272       0.6212       1.0507
        b_ωx      -1.0786      -0.7731       0.0923       0.7988       1.0790
        b_Ωy      -1.0744      -0.7970       0.0671       0.8168       1.0824
        b_Ωx      -1.0448      -0.5636       0.0511       0.5661       1.0411
         b_P      20.5943      25.7205      28.7099      32.3986      40.7267
        b_θy       0.0567       0.0666       0.0722       0.0788       0.0921
        b_θx      -1.2293      -1.0712      -0.9978      -0.9329      -0.8212
         b_ω      -2.9156      -1.8377       0.1651       1.2893       2.8894
         b_Ω      -3.0035      -2.3530       0.1043       0.7671       2.8975
         b_a       8.1272       9.2564       9.9031      10.7619      12.7081
         b_θ      -1.5087      -1.5019      -1.4983      -1.4947      -1.4881
        b_tp   40833.6948   43769.3959   44965.7566   45958.5231   47674.7839

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)