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_dat = 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)
)
astrom_obs = PlanetRelAstromObs(astrom_dat, name="simulated astrom")

planet_b = Planet(
    name="b",
    basis=Visual{KepOrbit},
    observations=[astrom_obs],
    variables=@variables begin
        a ~ truncated(Normal(10, 4), lower=0.1, upper=100)
        e ~ Uniform(0.0, 0.5)
        i ~ Sine()
        ω ~ UniformCircular()
        Ω ~ UniformCircular()
        θ ~ UniformCircular()
        tp = θ_at_epoch_to_tperi(θ,50420; system.M, a, e, i, ω, Ω)  # reference epoch for θ. Choose an MJD date near your data.
    end
)

sys = System(
    name="Tutoria",
    companions=[planet_b],
    observations=[],
    variables=@variables begin
        M ~ truncated(Normal(1.2, 0.1), lower=0.1)
        plx ~ truncated(Normal(50.0, 0.02), lower=0.1)
    end
)

model = Octofitter.LogDensityModel(sys)

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     = 3.17 seconds
Compute duration  = 3.17 seconds
parameters        = M, plx, b_a, b_e, b_i, b_ωx, b_ωy, b_Ωx, b_Ωy, b_θx, b_θy, 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

Use `describe(chains)` for summary statistics and quantiles.

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(model, chain, :b, ii)

50-element Vector{Visual{KepOrbit{Float64}, Float64}}:
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(8.37, 0.255, 0.397, -1.1, 2.22, 45300.0, 1.19), 49.96811135412299, 4.127928806140017e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(13.4, 0.112, 0.792, 2.85, 0.635, 47500.0, 1.41), 49.98341132826016, 4.126665242843802e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(13.7, 0.0274, 0.81, 1.58, 0.375, 42200.0, 1.1), 50.0234554331363, 4.123361820192522e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.7, 0.276, 0.731, -2.31, 3.03, 40900.0, 1.2), 50.000928409413795, 4.1252195270890687e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(12.5, 0.0912, 0.567, 1.19, 3.69, 34600.0, 1.03), 49.98547487657874, 4.1264948818910606e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(13.1, 0.432, 0.734, 1.0, 5.18, 36300.0, 1.2), 50.057617039473854, 4.120547849579864e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.0, 0.242, 0.477, -2.66, 2.93, 40500.0, 1.13), 49.98355433317848, 4.1266534362919503e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(8.6, 0.205, 0.402, -0.291, 1.32, 44900.0, 1.16), 50.01872928411929, 4.123751426699768e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.7, 0.0513, 0.483, 0.88, 0.418, 43700.0, 1.16), 49.99594220106413, 4.1256309445590596e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(16.4, 0.0371, 0.938, -2.79, 0.373, 46700.0, 1.05), 49.99726170708446, 4.125522062698671e6)
 ⋮
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.0, 0.0621, 0.574, 1.19, 0.317, 43900.0, 1.16), 50.00757393246178, 4.1246713253010297e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(8.78, 0.315, 0.477, -1.79, 2.92, 44800.0, 1.14), 50.03796966986167, 4.1221657794666993e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.9, 0.233, 0.655, 0.76, 5.0, 38300.0, 1.2), 50.00734389561336, 4.1246902990420554e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(14.2, 0.303, 0.892, -0.126, 4.67, 50000.0, 1.21), 49.98799223115194, 4.1262870749698663e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(12.0, 0.118, 0.579, 0.878, 4.08, 36300.0, 1.12), 49.972214009264654, 4.1275899084410323e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.3, 0.00092, 0.659, -2.96, 4.44, 44300.0, 1.15), 50.02174910750462, 4.1235024749694536e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(9.78, 0.28, 0.52, 1.28, 5.59, 42900.0, 1.21), 50.01230349197991, 4.124281263712908e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.7, 0.0601, 0.666, 1.44, 0.318, 44300.0, 1.35), 50.008313552679645, 4.1246103216380887e6)
 Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.7, 0.0641, 0.644, 0.492, 0.757, 43300.0, 1.12), 50.00122949893774, 4.125194686492227e6)

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

using CairoMakie

epochs = astrom_obs.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_obs.table.ra, astrom_obs.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

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.