Prior Predictive Checks

The prior predictive distributin of a Bayesian model what you get by sampling parameters directly from the priors and calculating where the model would place the data. For example, if sampling from relative astrometry, the prior predictive model is the distribution of (simulated) astrometry points corresponding to orbits drawn from the prior. For radial velocity data, these would be simulated RV points based on an RV curve drawn from the priors.

To generate a prior predictive distribution, one first needs to create a model. We will use the model and sample data from the Fit Astrometry tutorial:

using Octofitter
using CairoMakie
using PairPlots
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(50.0, 8),
    σ_dec = fill(50.0, 8),
    cor = fill(0.0, 8)
)
astrom_obs = PlanetRelAstromObs(astrom_dat, name="relastrom")
planet_b = Planet(
    name="b",
    basis=Visual{KepOrbit},
    observations=[astrom_obs],
    variables=@variables begin
        M = system.M
        a ~ truncated(Normal(10, 4), lower=0, upper=100)
        e ~ Uniform(0.0, 0.5)
        i ~ Sine()
        ω ~ UniformCircular()
        Ω ~ UniformCircular()
        θ_x ~ Normal()
        θ_y ~ Normal()
        θ = atan(θ_y, θ_x)
        tp = θ_at_epoch_to_tperi(θ, 50420; M, e, a, 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
)

We can now draw one sample from the prior:

prior_draw_system = generate_from_params(sys)
prior_draw_astrometry = prior_draw_system.planets.b.observations[1]

PlanetRelAstromObs Table with 5 columns and 8 rows:
     epoch  ra        dec       σ_ra  σ_dec
   ┌───────────────────────────────────────
 1 │ 50000  -72.4139  154.007   50.0  50.0
 2 │ 50120  -47.1304  133.54    50.0  50.0
 3 │ 50240  -20.0197  107.933   50.0  50.0
 4 │ 50360  7.7989    78.5787   50.0  50.0
 5 │ 50480  35.3856   46.7408   50.0  50.0
 6 │ 50600  61.9601   13.5402   50.0  50.0
 7 │ 50720  86.8766   -20.0344  50.0  50.0
 8 │ 50840  109.598   -53.109   50.0  50.0

And plot the generated astrometry:

Makie.scatter(prior_draw_astrometry.table.ra, prior_draw_astrometry.table.dec,color=:black, axis=(;autolimitaspect=1,xreversed=true))

We can repeat this many times to get a feel for our chosen priors in the domain of our data:

using Random
Random.seed!(1)


fig = Figure()
ax = Axis(
    fig[1,1], xlabel="ra offset [mas]", ylabel="dec offset [mas]",
    xreversed=true,
    aspect=1
)
for i in 1:50
    prior_draw_system = generate_from_params(sys)
    prior_draw_astrometry = prior_draw_system.planets.b.observations[1]
    Makie.scatter!(
        ax,
        prior_draw_astrometry.table.ra,
        prior_draw_astrometry.table.dec,
        color=Makie.cgrad(:turbo)[i/50],
    )
end


Makie.errorbars!(ax,astrom_dat.ra,astrom_dat.dec,astrom_dat.σ_dec,color=:black,linewidth=3)
Makie.errorbars!(ax,astrom_dat.ra,astrom_dat.dec,astrom_dat.σ_ra,direction=:x,color=:black,linewidth=3)

fig

The heavy black crosses are our actual data, while the colored ones are simulations drawn from our priors. Notice that our real data lies at a greater separation than most draws from the prior? That might mean the priors could be tweaked.