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.37 seconds
Compute duration = 3.37 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
Summary Statistics
parameters mean std mcse ess_bulk ess_tail r ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Floa ⋯
M 1.1873 0.0883 0.0033 718.5044 722.0107 0.9 ⋯
plx 49.9999 0.0201 0.0006 1045.9552 476.3477 0.9 ⋯
b_a 11.7600 2.1190 0.1310 253.3586 322.4434 1.0 ⋯
b_e 0.1477 0.1052 0.0067 254.5021 555.2656 1.0 ⋯
b_i 0.6460 0.1470 0.0094 274.5324 223.3653 1.0 ⋯
b_ωx 0.0256 0.7466 0.0550 208.7238 554.0880 0.9 ⋯
b_ωy 0.0334 0.6737 0.0578 148.5245 640.6825 1.0 ⋯
b_Ωx 0.0917 0.7542 0.1048 68.7277 424.5998 1.0 ⋯
b_Ωy 0.0879 0.6620 0.0800 82.5368 412.8118 1.0 ⋯
b_θx 0.0750 0.0102 0.0003 893.4947 675.8711 0.9 ⋯
b_θy -1.0022 0.0954 0.0037 664.6501 722.0107 0.9 ⋯
b_ω -0.1160 1.8140 0.1397 230.1912 522.4317 1.0 ⋯
b_Ω -0.3547 1.7099 0.1997 109.7777 313.3970 1.0 ⋯
b_θ -1.4961 0.0077 0.0002 1022.3334 741.3550 0.9 ⋯
b_tp 43428.6886 4758.3669 277.3695 298.1882 344.6987 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.0151 1.1296 1.1885 1.2497 1.3523
plx 49.9609 49.9865 49.9997 50.0133 50.0418
b_a 8.4653 10.1590 11.4897 13.0444 16.5646
b_e 0.0067 0.0607 0.1295 0.2180 0.3874
b_i 0.3039 0.5605 0.6631 0.7544 0.8845
b_ωx -1.0707 -0.7556 0.1146 0.7774 1.0744
b_ωy -1.0509 -0.5724 0.0225 0.6778 1.0569
b_Ωx -1.0789 -0.7088 0.2300 0.8220 1.0921
b_Ωy -1.0152 -0.4970 0.2034 0.6637 1.0521
b_θx 0.0559 0.0682 0.0743 0.0818 0.0960
b_θy -1.1992 -1.0623 -0.9996 -0.9357 -0.8231
b_ω -3.0305 -1.8010 0.0651 1.2052 2.9832
b_Ω -2.9941 -2.1182 0.2637 0.8884 2.8076
b_θ -1.5107 -1.5014 -1.4962 -1.4908 -1.4811
b_tp 32621.6606 40864.2312 44433.0924 46896.8460 50032.3173
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(10.4, 0.0605, 0.592, 2.36, 1.21, 48400.0, 1.26), 50.00054097470339, 4.125251491807883e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(9.1, 0.177, 0.398, 1.07, 0.285, 45000.0, 1.11), 50.0255597050733, 4.1231883753652074e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.1, 0.0895, 0.484, -2.69, 3.82, 43600.0, 1.07), 50.00966538480354, 4.1244988275761213e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(16.3, 0.0826, 0.857, 1.21, 3.36, 49700.0, 1.24), 49.9883604255012, 4.126256682383043e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.9, 0.00997, 0.655, 0.679, 0.483, 42400.0, 1.22), 49.96955409492288, 4.127809622941057e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(9.85, 0.0302, 0.615, 0.665, 1.28, 45700.0, 1.2), 49.9809509997178, 4.1268683792803576e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(14.9, 0.385, 0.602, 0.0599, 6.02, 33200.0, 1.19), 49.99977267759148, 4.125314880472218e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(12.8, 0.0417, 0.765, 1.09, 0.327, 42000.0, 1.21), 50.01107294155081, 4.124382743960786e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(9.79, 0.0281, 0.603, -0.305, 1.42, 44100.0, 1.13), 50.03455280937631, 4.1224472822397854e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(9.98, 0.169, 0.594, 1.1, 0.426, 44800.0, 1.25), 49.95435475420174, 4.129065569198391e6)
⋮
Visual{KepOrbit{Float64}, Float64}(KepOrbit(14.3, 0.0127, 0.874, 2.01, 0.528, 43900.0, 1.29), 49.95205644509205, 4.1292555487445313e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(15.3, 0.335, 0.831, 3.07, 1.35, 49700.0, 1.28), 50.00437897635054, 4.124934865097492e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.1, 0.112, 0.647, 1.41, 0.881, 46100.0, 1.15), 50.01671242780641, 4.1239177114012996e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.9, 0.0371, 0.637, 0.966, 3.86, 50300.0, 1.15), 49.9832450326463, 4.1266789723711526e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(12.3, 0.148, 0.71, -2.55, 3.17, 40500.0, 1.24), 50.00797999880363, 4.1246378328424976e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(13.2, 0.204, 0.643, -2.65, 1.27, 34100.0, 1.13), 49.98654503045293, 4.126406538428195e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.9, 0.189, 0.546, 0.538, 6.25, 41200.0, 1.12), 49.999535001967836, 4.1253344903903278e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(13.7, 0.0737, 0.695, 1.78, 3.35, 33400.0, 1.08), 50.016484810725395, 4.12393647869603e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.4, 0.0496, 0.634, -0.734, 3.84, 46900.0, 1.13), 49.992989136501514, 4.1258746438207207e6)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.