Posterior Predictive Checks

A posterior predictive check compares our true data with simulated data drawn from the posterior. This allows us to evaluate if the model is able to reproduce our observations appropriately. Samples drawn from the posterior predictive distribution should match the locations of the original data.

To demonstrate, we will fit a model to relative astrometry data:

using Octofitter
using Distributions

astrom_like = PlanetRelAstromLikelihood(Table(;
    epoch= [50000,50120,50240,50360,50480,50600,50720,50840,],
    ra = [-505.764,-502.57,-498.209,-492.678,-485.977,-478.11,-469.08,-458.896,],
    dec = [-66.9298,-37.4722,-7.92755,21.6356, 51.1472, 80.5359, 109.729, 138.651,],
    σ_ra = fill(10.0, 8),
    σ_dec = fill(10.0, 8),
    cor = fill(0.0, 8)
))
@planet b Visual{KepOrbit} begin
    a ~ truncated(Normal(10, 4), lower=0, upper=100)
    e ~ Uniform(0.0, 0.5)
    i ~ Sine()
    ω ~ UniformCircular()
    Ω ~ UniformCircular()
    θ ~ UniformCircular()
    tp = θ_at_epoch_to_tperi(system,b,50420)  # reference epoch for θ. Choose an MJD date near your data.
end astrom_like
@system Tutoria begin
    M ~ truncated(Normal(1.2, 0.1), lower=0)
    plx ~ truncated(Normal(50.0, 0.02), lower=0)
end b
model = Octofitter.LogDensityModel(Tutoria)

using Random
Random.seed!(0)
chain = octofit(model)
Chains MCMC chain (1000×28×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 5.85 seconds
Compute duration  = 5.85 seconds
parameters        = M, plx, b_a, b_e, b_i, b_ωy, b_ωx, b_Ωy, b_Ωx, b_θy, b_θx, b_ω, b_Ω, 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      rh ⋯
      Symbol      Float64     Float64    Float64    Float64    Float64   Float ⋯

           M       1.1786      0.0968     0.0039   615.4503   493.5200    0.99 ⋯
         plx      49.9998      0.0185     0.0006   903.2070   782.7594    1.00 ⋯
         b_a      11.9229      2.4987     0.3018    97.4455   233.3716    1.02 ⋯
         b_e       0.1554      0.1112     0.0073   248.1913   213.0305    1.01 ⋯
         b_i       0.6356      0.1671     0.0133   139.9958   148.0941    1.02 ⋯
        b_ωy      -0.0221      0.7485     0.0605   174.1941   677.9194    1.00 ⋯
        b_ωx       0.0188      0.6788     0.0483   222.0595   639.7208    1.01 ⋯
        b_Ωy      -0.1289      0.7340     0.1275    42.4438   381.2656    1.04 ⋯
        b_Ωx      -0.1149      0.6695     0.0955    65.3341   499.6503    1.01 ⋯
        b_θy       0.0741      0.0105     0.0004   856.4565   609.3305    1.00 ⋯
        b_θx      -0.9997      0.1000     0.0035   819.3240   631.6835    1.00 ⋯
         b_ω      -0.0624      1.8601     0.1161   326.9396   574.1878    1.00 ⋯
         b_Ω      -0.8725      1.7599     0.2498    73.2361   232.7969    1.01 ⋯
         b_θ      -1.4968      0.0072     0.0002   895.0744   662.9977    1.00 ⋯
        b_tp   42893.2184   5510.3072   321.3524   314.8469   268.5010    1.02 ⋯
                                                               2 columns omitted

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

           M       0.9932       1.1156       1.1809       1.2451       1.3611
         plx      49.9649      49.9867      49.9997      50.0131      50.0353
         b_a       8.2422      10.1200      11.5010      13.2805      18.2017
         b_e       0.0057       0.0633       0.1296       0.2265       0.4246
         b_i       0.2428       0.5446       0.6576       0.7532       0.9092
        b_ωy      -1.0813      -0.7826      -0.0811       0.7373       1.0816
        b_ωx      -1.0478      -0.6258       0.0393       0.6422       1.0394
        b_Ωy      -1.0947      -0.8418      -0.2751       0.6458       1.0373
        b_Ωx      -1.0652      -0.7051      -0.2830       0.5196       1.0071
        b_θy       0.0557       0.0668       0.0738       0.0806       0.0958
        b_θx      -1.2049      -1.0629      -0.9962      -0.9286      -0.8197
         b_ω      -2.9956      -1.8527       0.0945       1.3880       2.9622
         b_Ω      -3.0292      -2.5166      -1.5728       0.6648       2.7011
         b_θ      -1.5100      -1.5019      -1.4969      -1.4917      -1.4830
        b_tp   28946.9835   40296.1702   44360.3166   46794.4787   49979.6066

We now have our posterior as approximated by the MCMC chain. Convert these posterior samples into orbit objects:

# Instead of creating orbit objects for all rows in the chain, just pick
# every twentieth row.
ii = 1:20:1000
orbits = Octofitter.construct_elements(chain, :b, ii)
50-element Vector{Visual{KepOrbit{Float64}, Float64}}:
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(8.12, 0.348, 0.51, -1.14, 2.16, 45200.0, 1.17), 50.006499281962306, 4.1247638399355267e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(9.52, 0.0493, 0.458, 0.309, 1.21, 45000.0, 1.11), 50.007498447661675, 4.124681425844144e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(15.6, 0.265, 0.761, 0.267, 3.84, 48900.0, 1.2), 50.012922328322304, 4.1242341058561215e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(17.6, 0.48, 0.641, 0.45, 5.04, 25600.0, 1.09), 49.98170668335251, 4.1268098607905484e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(15.8, 0.345, 0.865, 2.88, 0.985, 48500.0, 1.23), 50.01613173657431, 4.1239694642192544e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(16.3, 0.0828, 0.885, 0.633, 3.39, 47900.0, 1.17), 49.99498147889525, 4.1257140996656273e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.3, 0.015, 0.305, 1.15, 3.44, 50000.0, 1.09), 50.02284415716656, 4.1234160806997875e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(13.1, 0.187, 0.793, -0.528, 3.86, 47500.0, 1.31), 49.96920810596063, 4.127842081519708e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.6, 0.0333, 0.664, -2.9, 3.8, 41700.0, 1.16), 50.01366547160223, 4.124172824667636e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.6, 0.118, 0.759, -0.954, 3.92, 47000.0, 1.33), 50.00286303821135, 4.1250637956945733e6)
 ⋮
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.5, 0.00462, 0.552, 2.49, 1.03, 47900.0, 1.1), 49.98702139606487, 4.126371090721517e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(8.77, 0.203, 0.363, -0.819, 1.74, 44600.0, 1.14), 50.00690239142661, 4.124730589898784e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.0, 0.00456, 0.58, 2.47, 3.96, 40900.0, 1.14), 49.99067379809432, 4.1260696111654127e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.0, 0.0715, 0.545, -1.5, 4.16, 47000.0, 1.27), 50.009575813364336, 4.1245100892233257e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(9.92, 0.0211, 0.551, -2.88, 4.42, 44900.0, 1.17), 50.00380560677614, 4.1249860385036077e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(8.45, 0.199, 0.391, 2.86, 4.56, 45200.0, 1.13), 50.019108875334304, 4.1237240054413388e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.5, 0.157, 0.58, 0.606, 4.71, 37500.0, 1.09), 50.024900432863, 4.1232465875034058e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(15.1, 0.121, 0.797, -0.508, 0.0563, 32300.0, 1.04), 49.989282992231566, 4.126184407006878e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.0, 0.028, 0.545, -0.575, 1.24, 43000.0, 1.09), 50.00657689676722, 4.1247574379228195e6)

Calculate and plot the location the planet would be at each observation epoch:

using CairoMakie

epochs = astrom_like.table.epoch' # transpose

x = raoff.(orbits, epochs)[:]
y = decoff.(orbits, epochs)[:]

fig = Figure()
ax = Axis(
    fig[1,1], xlabel="ra offset [mas]", ylabel="dec offset [mas]",
    xreversed=true,
    aspect=1
)
for orbit in orbits
    Makie.lines!(ax, orbit, color=:lightgrey)
end

Makie.scatter!(
    ax,
    x, y,
    markersize=3,
)
fig

Makie.scatter!(ax, astrom_like.table.ra, astrom_like.table.dec,color=:black, label="observed")
Makie.scatter!(ax, [0],[0], marker='⋆', color=:black, markersize=20,label="")
Makie.xlims!(400,-700)
Makie.ylims!(-200,200)
fig
Example block output

Looks like a great match to the data! Notice how the uncertainty around the middle point is lower than the ends. That's because the orbit's posterior location at that epoch is also constrained by the surrounding data points. We can know the location of the planet in hindsight better than we could measure it!

You can follow this same procedure for any kind of data modelled with Octofitter.