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 = 3.88 seconds
Compute duration = 3.88 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 r ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Floa ⋯
M 1.1863 0.0978 0.0063 237.3985 446.4125 1.0 ⋯
plx 49.9998 0.0183 0.0006 1073.2701 614.9544 1.0 ⋯
b_a 11.9373 2.2558 0.1411 244.7130 434.4190 1.0 ⋯
b_e 0.1507 0.1088 0.0059 346.6565 464.3563 1.0 ⋯
b_i 0.6484 0.1595 0.0136 182.2556 48.3513 1.0 ⋯
b_ωy -0.1467 0.7151 0.0606 152.3857 553.8900 1.0 ⋯
b_ωx 0.0439 0.6995 0.0517 191.5746 627.1462 1.0 ⋯
b_Ωy 0.1647 0.7489 0.1014 72.7548 351.9880 1.0 ⋯
b_Ωx 0.2064 0.6242 0.0790 75.2757 277.2631 1.0 ⋯
b_θy 0.0745 0.0104 0.0004 594.4027 664.9113 1.0 ⋯
b_θx -1.0018 0.1022 0.0041 599.0980 646.8931 1.0 ⋯
b_ω -0.1394 1.9782 0.1197 391.5292 656.9181 1.0 ⋯
b_Ω -0.0124 1.6764 0.1724 212.5306 391.0570 1.0 ⋯
b_θ -1.4966 0.0073 0.0002 946.5030 445.8015 1.0 ⋯
b_tp 43169.6849 4963.2786 297.2918 282.7560 478.1124 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.0017 1.1205 1.1865 1.2460 1.3890
plx 49.9657 49.9873 49.9995 50.0126 50.0370
b_a 8.5265 10.3040 11.5296 13.1584 17.4011
b_e 0.0083 0.0596 0.1272 0.2237 0.3900
b_i 0.2524 0.5636 0.6605 0.7610 0.9090
b_ωy -1.1099 -0.8182 -0.3119 0.5708 1.0193
b_ωx -1.0261 -0.6289 0.0438 0.7191 1.0612
b_Ωy -1.0562 -0.6371 0.3930 0.8619 1.0715
b_Ωx -0.9877 -0.3066 0.3619 0.7359 1.0482
b_θy 0.0561 0.0671 0.0736 0.0814 0.0958
b_θx -1.2140 -1.0694 -0.9994 -0.9297 -0.8196
b_ω -3.0535 -2.1561 0.1623 1.5202 3.0482
b_Ω -2.9961 -1.6903 0.4344 1.0967 2.9229
b_θ -1.5118 -1.5014 -1.4967 -1.4917 -1.4821
b_tp 31296.7229 40292.9830 44175.3197 46940.4300 49827.1601
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(13.7, 0.119, 0.75, -3.11, 0.49, 47800.0, 1.28), 49.99961318468417, 4.1253319148313506e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.2, 0.00517, 0.613, 1.44, 1.21, 46700.0, 1.25), 50.006450964884905, 4.124767825352005e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(13.5, 0.0218, 0.789, 2.14, 0.356, 44100.0, 1.15), 50.020366220760415, 4.1236203487529033e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.2, 0.278, 0.571, 1.33, 5.23, 41800.0, 1.2), 49.98660625505593, 4.1264053604186657e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(14.8, 0.401, 0.824, -2.65, 2.96, 35300.0, 1.31), 49.99422346501079, 4.125776653863974e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(9.95, 0.0674, 0.384, 1.96, 4.4, 41000.0, 0.896), 50.0165692328313, 4.1239333917490267e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(8.65, 0.269, 0.444, 1.39, 6.2, 45000.0, 1.01), 49.978123956661754, 4.127105694860846e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(12.6, 0.00953, 0.705, -1.53, 0.479, 36100.0, 1.14), 50.00977092362149, 4.1244939976834273e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(17.5, 0.13, 0.925, -2.38, 0.387, 48700.0, 1.35), 49.99321183838555, 4.1258601401085937e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(9.28, 0.0847, 0.445, 2.29, 4.86, 44100.0, 1.12), 49.96741956577682, 4.127989834025228e6)
⋮
Visual{KepOrbit{Float64}, Float64}(KepOrbit(12.1, 0.103, 0.805, 2.04, 0.746, 46200.0, 1.32), 49.97860449106273, 4.127066013555554e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.2, 0.172, 0.402, -2.57, 3.13, 41700.0, 0.979), 49.99858141682606, 4.1254170449441085e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(12.4, 0.051, 0.745, 3.03, 3.7, 39800.0, 1.19), 50.01856340882588, 4.1237689758119704e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.0, 0.171, 0.734, 1.05, 0.44, 43800.0, 1.24), 50.00071388512431, 4.1252411010348755e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.6, 0.14, 0.627, 2.87, 1.2, 49000.0, 1.06), 50.004343280850705, 4.124941684395438e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(11.8, 0.0729, 0.66, 0.72, 0.292, 41900.0, 1.19), 49.999263683036055, 4.125360751462074e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(10.3, 0.24, 0.691, -1.66, 3.25, 44600.0, 1.16), 50.020186470918155, 4.1236351671724156e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(16.5, 0.0644, 0.907, 1.58, 3.48, 28200.0, 1.15), 50.00467916186937, 4.1249139771960694e6)
Visual{KepOrbit{Float64}, Float64}(KepOrbit(9.37, 0.198, 0.489, -1.31, 1.88, 43600.0, 1.23), 49.998395498995784, 4.1254323852080978e6)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
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.