Cross-Validation
Cross-validation asks how well a fitted model predicts data it has not seen. In Octofitter that question is answered by subsetting observations: every observation type that carries a data table knows how to produce a copy of itself restricted to a chosen set of rows, and the machinery below builds derived models out of those copies.
There are two granularities:
- Pointwise — one held-out data row at a time. This is what leave-one-out cross-validation and PSIS-LOO need, and
Octofitter.pointwise_likecomputes it from a single chain with no refitting. - Whole likelihoods, or groups of epochs — drop an entire instrument, keep only one, or feed the data in cumulatively. These produce new
Systems that you refit.
We will use a model with two relative-astrometry instruments so that both granularities have something to say:
using Octofitter
using Distributions
using CairoMakie
using Random
astrom_dat_gpi = Table(;
epoch = [50000., 50120, 50240, 50360],
ra = [-505.764, -502.57, -498.209, -492.678],
dec = [-66.9298, -37.4722, -7.92755, 21.6356],
σ_ra = fill(10.0, 4),
σ_dec = fill(10.0, 4),
cor = fill(0.0, 4),
)
astrom_dat_sphere = Table(;
epoch = [50480., 50600, 50720, 50840],
ra = [-485.977, -478.11, -469.08, -458.896],
dec = [51.1472, 80.5359, 109.729, 138.651],
σ_ra = fill(10.0, 4),
σ_dec = fill(10.0, 4),
cor = fill(0.0, 4),
)
A = Body(
name="A",
variables=@variables begin
mass ~ truncated(Normal(1.2, 0.1), lower=0.1) # M⊙
end
)
b = Body(
name="b",
about=A,
variables=@variables begin
a ~ truncated(Normal(10, 4), lower=0.1, upper=100)
e ~ Uniform(0.0, 0.5)
i ~ Sine()
ω ~ UniformCircular()
Ω ~ UniformCircular()
θ ~ UniformCircular()
epoch = 50420.0
end
)
gpi = RelAstromObs(astrom_dat_gpi; target=b, ref=A, name="GPI")
sphere = RelAstromObs(astrom_dat_sphere; target=b, ref=A, name="SPHERE")
sys = System(
name="Tutoria",
bodies=[A, b],
observations=[gpi, sphere],
variables=@variables begin
plx ~ truncated(Normal(50.0, 0.02), lower=0.1)
end
)
model = Octofitter.LogDensityModel(sys)
Random.seed!(0)
chain = octofit(model)
display(chain)[ Info: [Tutoria] observing_geometry = false (auto): worst accumulated bias 0.00535σ, 18.7× inside the 0.1σ limit — over 300 prior draws (seed 0xc70f177e5000001)
[ Info: [Tutoria] barycentric_lighttime = false (auto): changes no prediction at all — over 300 prior draws (seed 0xc70f177e5000001)
[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 11
┌ Info: Determined number type
└ T = Float64
ℓπcallback(θ): 0.000009 seconds
∇ℓπcallback(θ): 0.000025 seconds (1 allocation: 32 bytes)
┌ Info: Starting values not provided for all parameters! Guessing starting point using global optimization:
│ num_params = 11
└ num_fixed = 0
┌ Warning: Verbosity toggle: unrecognized_stop_reason
│ Unrecognized stop reason: Too many steps (101) without any function evaluations (probably search has converged). Defaulting to ReturnCode.Default.
└ @ OptimizationBase ~/.julia/packages/OptimizationBase/Jfw5O/src/utils.jl:170
┌ Info: Found sample of initial positions
│ logpost_range = (-58.67906911190986, -49.421754602529056)
└ mean_logpost = -52.54181566495545
[ Info: Sampling, beginning with adaptation phase...
Sampling 6%|██ | ETA: 0:00:02
iterations: 126
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.11
n_steps: 511
is_accept: true
acceptance_rate: 0.99668054377031
log_density: -52.976883771368065
hamiltonian_energy: 56.858482299755835
hamiltonian_energy_error: -0.001500454390011896
max_hamiltonian_energy_error: -0.22541906682441493
tree_depth: 9
numerical_error: false
step_size: 0.029251717125393477
nom_step_size: 0.029251717125393477
is_adapt: true
mass_matrix: DenseEuclideanMetric([0.000166793, 0.00552537, 0 ...])
Sampling 18%|█████▋ | ETA: 0:00:03
iterations: 362
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.07
n_steps: 63
is_accept: true
acceptance_rate: 0.6181754526760976
log_density: -53.25274932168541
hamiltonian_energy: 56.80291402917668
hamiltonian_energy_error: 0.676591459790572
max_hamiltonian_energy_error: 1.4693248999375328
tree_depth: 6
numerical_error: false
step_size: 0.06240966231504271
nom_step_size: 0.06240966231504271
is_adapt: true
mass_matrix: DenseEuclideanMetric([4.7762e-5, 0.01027, 0.0406 ...])
Sampling 29%|█████████ | ETA: 0:00:02
iterations: 579
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.07
n_steps: 127
is_accept: true
acceptance_rate: 0.9890373999938042
log_density: -61.810274158721825
hamiltonian_energy: 64.93868813163071
hamiltonian_energy_error: 0.023400524870965
max_hamiltonian_energy_error: -0.05686342239101805
tree_depth: 7
numerical_error: false
step_size: 0.024869148711011094
nom_step_size: 0.024869148711011094
is_adapt: true
mass_matrix: DenseEuclideanMetric([2.45567e-5, 0.00638275, 0. ...])
Sampling 36%|███████████▍ | ETA: 0:00:02
iterations: 730
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 127
is_accept: true
acceptance_rate: 0.9303838741815484
log_density: -61.96944267139411
hamiltonian_energy: 66.23164191840269
hamiltonian_energy_error: 0.18486408864788473
max_hamiltonian_energy_error: -0.5338990064359308
tree_depth: 7
numerical_error: false
step_size: 0.029579646653874792
nom_step_size: 0.029579646653874792
is_adapt: true
mass_matrix: DenseEuclideanMetric([2.45567e-5, 0.00638275, 0. ...])
Sampling 45%|█████████████▉ | ETA: 0:00:01
iterations: 899
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 63
is_accept: true
acceptance_rate: 0.847987918479339
log_density: -54.2356385302678
hamiltonian_energy: 59.403395917444264
hamiltonian_energy_error: -0.941901656997679
max_hamiltonian_energy_error: 2.3558195560499655
tree_depth: 6
numerical_error: false
step_size: 0.052925282550564494
nom_step_size: 0.052925282550564494
is_adapt: true
mass_matrix: DenseEuclideanMetric([2.45567e-5, 0.00638275, 0. ...])
Sampling 50%|███████████████▋ | ETA: 0:00:01
iterations: 1010
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 255
is_accept: true
acceptance_rate: 1.0
log_density: -51.50311776516283
hamiltonian_energy: 55.66964400149361
hamiltonian_energy_error: -0.04460341076187291
max_hamiltonian_energy_error: -0.12295636284893874
tree_depth: 7
numerical_error: false
step_size: 0.0232698008995888
nom_step_size: 0.0232698008995888
is_adapt: false
mass_matrix: DenseEuclideanMetric([1.00673e-5, 0.00851224, 0. ...])
Sampling 57%|█████████████████▋ | ETA: 0:00:01
iterations: 1136
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 127
is_accept: true
acceptance_rate: 0.9990754717031346
log_density: -58.586942325022534
hamiltonian_energy: 61.70239517311105
hamiltonian_energy_error: -0.03155832692620919
max_hamiltonian_energy_error: -0.05393025137446017
tree_depth: 7
numerical_error: false
step_size: 0.0232698008995888
nom_step_size: 0.0232698008995888
is_adapt: false
mass_matrix: DenseEuclideanMetric([1.00673e-5, 0.00851224, 0. ...])
Sampling 63%|███████████████████▌ | ETA: 0:00:01
iterations: 1259
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 127
is_accept: true
acceptance_rate: 0.994333682303415
log_density: -53.02230379460754
hamiltonian_energy: 60.88958969656021
hamiltonian_energy_error: -0.011748946609671407
max_hamiltonian_energy_error: -0.08380183091199456
tree_depth: 7
numerical_error: false
step_size: 0.0232698008995888
nom_step_size: 0.0232698008995888
is_adapt: false
mass_matrix: DenseEuclideanMetric([1.00673e-5, 0.00851224, 0. ...])
Sampling 69%|█████████████████████▎ | ETA: 0:00:01
iterations: 1371
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 63
is_accept: true
acceptance_rate: 0.9749188875992072
log_density: -53.88086581029095
hamiltonian_energy: 58.50002167324439
hamiltonian_energy_error: 0.03175878698564816
max_hamiltonian_energy_error: 0.05975585812588946
tree_depth: 5
numerical_error: false
step_size: 0.0232698008995888
nom_step_size: 0.0232698008995888
is_adapt: false
mass_matrix: DenseEuclideanMetric([1.00673e-5, 0.00851224, 0. ...])
Sampling 74%|███████████████████████ | ETA: 0:00:01
iterations: 1480
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 255
is_accept: true
acceptance_rate: 0.9187999048044264
log_density: -53.4313527555426
hamiltonian_energy: 58.03899198289097
hamiltonian_energy_error: 0.09957346627491148
max_hamiltonian_energy_error: 0.4248712849198242
tree_depth: 7
numerical_error: false
step_size: 0.0232698008995888
nom_step_size: 0.0232698008995888
is_adapt: false
mass_matrix: DenseEuclideanMetric([1.00673e-5, 0.00851224, 0. ...])
Sampling 80%|████████████████████████▋ | ETA: 0:00:00
iterations: 1590
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 63
is_accept: true
acceptance_rate: 0.9934245416436592
log_density: -55.46365575312548
hamiltonian_energy: 60.600751422428274
hamiltonian_energy_error: -0.030889044523135567
max_hamiltonian_energy_error: -0.06858871187085214
tree_depth: 5
numerical_error: false
step_size: 0.0232698008995888
nom_step_size: 0.0232698008995888
is_adapt: false
mass_matrix: DenseEuclideanMetric([1.00673e-5, 0.00851224, 0. ...])
Sampling 85%|██████████████████████████▎ | ETA: 0:00:00
iterations: 1696
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 127
is_accept: true
acceptance_rate: 0.9884022140385362
log_density: -56.886870150675534
hamiltonian_energy: 69.27361157346955
hamiltonian_energy_error: -0.14295618348931782
max_hamiltonian_energy_error: -0.4828002672769145
tree_depth: 7
numerical_error: false
step_size: 0.0232698008995888
nom_step_size: 0.0232698008995888
is_adapt: false
mass_matrix: DenseEuclideanMetric([1.00673e-5, 0.00851224, 0. ...])
Sampling 90%|████████████████████████████ | ETA: 0:00:00
iterations: 1808
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 127
is_accept: true
acceptance_rate: 0.9859230848943523
log_density: -53.072673976481575
hamiltonian_energy: 60.01933725123446
hamiltonian_energy_error: 0.02013727212997196
max_hamiltonian_energy_error: 0.038799098817335675
tree_depth: 7
numerical_error: false
step_size: 0.0232698008995888
nom_step_size: 0.0232698008995888
is_adapt: false
mass_matrix: DenseEuclideanMetric([1.00673e-5, 0.00851224, 0. ...])
Sampling 96%|█████████████████████████████▊ | ETA: 0:00:00
iterations: 1918
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 127
is_accept: true
acceptance_rate: 0.9993375272108742
log_density: -55.94716876856044
hamiltonian_energy: 58.93058233313833
hamiltonian_energy_error: 0.011385276452365645
max_hamiltonian_energy_error: -0.3355956054061835
tree_depth: 7
numerical_error: false
step_size: 0.0232698008995888
nom_step_size: 0.0232698008995888
is_adapt: false
mass_matrix: DenseEuclideanMetric([1.00673e-5, 0.00851224, 0. ...])
Sampling 100%|███████████████████████████████| Time: 0:00:02
iterations: 2000
ratio_divergent_transitions: 0.0
ratio_divergent_transitions_during_adaption: 0.06
n_steps: 127
is_accept: true
acceptance_rate: 0.9906855246995395
log_density: -54.01185953339626
hamiltonian_energy: 58.19613972519258
hamiltonian_energy_error: -0.0909103762306529
max_hamiltonian_energy_error: -0.6662004438883642
tree_depth: 7
numerical_error: false
step_size: 0.0232698008995888
nom_step_size: 0.0232698008995888
is_adapt: false
mass_matrix: DenseEuclideanMetric([1.00673e-5, 0.00851224, 0. ...])
[ Info: Sampling compete. Building chains.
Sampling report for chain:
mean_accept = 0.9598760932119917
ratio_divergent_transitions = 0.0
mean_tree_depth = 6.813
max_tree_depth_frac = 0.0
total_steps = 315840
μs/step (approx.) = 7.57Calculating Pointwise Likelihoods
After you have defined a model and sampled from its posterior (e.g. via octofit), you can see how each datapoint is influencing the posterior:
@time likelihood_mat, epochs = Octofitter.pointwise_like(model, chain)
size(likelihood_mat)(1000, 8)likelihood_mat is an N_sample × N_data matrix, and epochs labels each column with the epoch of the data row it came from.
The columns are ordered exactly as the data are defined in the model: observation by observation in the order you listed them in System(...; observations=[...]), and within each observation, row by row in table order.
epochs8-element Vector{Float64}:
50000.0
50120.0
50240.0
50360.0
50480.0
50600.0
50720.0
50840.0A ~ line written inside a @variables block, an LL += ... line, and the UnitLengthPrior that sits behind every UniformCircular are all prior terms: they reshape the prior rather than adding data, so they do not participate in this machinery and get no column.
ObsPriorONeil2019 is not a prior-shaped term in this sense — it wraps a real likelihood — so it still contributes one column per row of the observation it wraps.
The consequence is that the columns sum to the model's log-likelihood minus those prior terms — which is the quantity used by PSIS-LOO:
θ = Octofitter.mcmcchain2result(model, chain, 1)
lnlike = Octofitter.make_ln_like(model.system)
(sum(likelihood_mat[1, :]), lnlike(model.system, θ))(-53.699912453982805, -50.924538084313475)The difference is exactly the three UnitLengthPrior terms this model's three UniformCircular variables contribute.
Observations with no epochs get one column, labelled NaN. PhotometryObs carries data but has no epoch column; it contributes one column per photometry row, and those columns are labelled NaN in epochs.
Pareto-Smoothed Importance Sampling
Once you have likelihood_mat you can use the Julia package ParetoSmooth.jl to efficiently calculate a leave-one-out cross-validation score. This technique takes a single posterior chain and, using the pointwise likelihoods, estimates what the posterior would have been had each datapoint in turn been held out.
In broad terms, one might say that this test verifies that no individual datapoints are overly skewing the results.
psis_loo comes from ParetoSmooth.jl, not from Octofitter; Octofitter supplies the matrix it consumes. Note the transpose — ParetoSmooth wants N_data × N_sample:
using ParetoSmooth
result = psis_loo(
collect(likelihood_mat'),
chain_index=ones(Int, size(chain, 1))
)
display(result)[ Info: No source provided for samples; variables are assumed to be from a Markov Chain. If the samples are independent, specify this with keyword argument `source=:other`.
Results of PSIS-LOO-CV with 1000 Monte Carlo samples and 8 data points. Total Monte Carlo SE of 0.049.
┌───────────┬────────┬──────────┬───────┬─────────┐
│ │ total │ se_total │ mean │ se_mean │
├───────────┼────────┼──────────┼───────┼─────────┤
│ cv_elpd │ -53.93 │ 0.47 │ -6.74 │ 0.06 │
│ naive_lpd │ -53.28 │ 0.25 │ -6.66 │ 0.03 │
│ p_eff │ 0.65 │ 0.22 │ 0.08 │ 0.03 │
└───────────┴────────┴──────────┴───────┴─────────┘chain_index tells ParetoSmooth which chain each row of the matrix came from. pointwise_like flattens all chains into the sample dimension in order, so for a single chain ones(Int, size(chain,1)) is right; for several chains use repeat(1:size(chain,3), inner=size(chain,1)).
result.pointwise is a KeyedArray of per-row diagnostics, indexed by name:
result.pointwise(:pareto_k)1-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓ data ∈ 8-element UnitRange{Int64}
And data, 8-element view(::Matrix{Float64}, :, 5) with eltype Float64:
(1) 0.46241676654220937
(2) 0.3504903942096959
(3) 0.2833591528034171
(4) 0.18647078303428194
(5) 0.22012215445828545
(6) -0.028079382504153812
(7) 0.037712711460642166
(8) 0.3476159271561565The available statistics are :cv_elpd (the leave-one-out expected log predictive density of that row), :naive_lpd (the same quantity without holding the row out), :p_eff (the difference between them — how much of that row the fit has effectively absorbed), :mcse, and :pareto_k.
The diagnostic to look at is the Pareto shape parameter $\hat{k}$ for each point. Values above about 0.7 mean the importance-sampling approximation is unreliable for that point — usually because it is highly influential — and that point deserves a genuine refit with the row held out (see Holding out whole observations below).
using CairoMakie
nrows = size(likelihood_mat, 2)
fig = Figure(size=(650, 620))
ax = Axis(
fig[1,1],
xlabel="data row",
ylabel="Pareto k̂",
xticks=1:nrows
)
scatter!(ax, collect(result.pointwise(:pareto_k)))
hlines!(ax, [0.7], color=:red, linestyle=:dash)
ylims!(ax, -0.1, 0.9) # keep the threshold visible even when nothing is near it
ax = Axis(
fig[2,1],
xlabel="data row",
ylabel="MCSE",
xticks=1:nrows
)
scatter!(ax, collect(result.pointwise(:mcse)))
ax = Axis(
fig[3,1],
xlabel="data row",
ylabel="p_eff",
xticks=1:nrows
)
scatter!(ax, collect(result.pointwise(:p_eff)))
fig
Every $\hat{k}$ here is comfortably below the line, which is what a healthy fit looks like: no single astrometric epoch is holding the posterior up on its own, so the LOO score above can be trusted as it stands.
Finding an influential point
The diagnostic is only interesting when something is wrong, so here is the same model with one SPHERE measurement displaced by 60 mas — six times its quoted uncertainty — standing in for a mis-registered frame or a background object mistaken for the companion:
astrom_dat_sphere_bad = Table(;
epoch = astrom_dat_sphere.epoch,
ra = astrom_dat_sphere.ra,
dec = [51.1472, 80.5359 + 60.0, 109.729, 138.651], # row 2 nudged
σ_ra = fill(10.0, 4),
σ_dec = fill(10.0, 4),
cor = fill(0.0, 4),
)
sys_bad = System(
name="Tutoria_outlier",
bodies=[A, b],
observations=[
gpi,
RelAstromObs(astrom_dat_sphere_bad; target=b, ref=A, name="SPHERE"),
],
variables=@variables begin
plx ~ truncated(Normal(50.0, 0.02), lower=0.1)
end
)
model_bad = Octofitter.LogDensityModel(sys_bad)
Random.seed!(0)
chain_bad = octofit(model_bad, verbosity=0)
mat_bad, epochs_bad = Octofitter.pointwise_like(model_bad, chain_bad)
result_bad = psis_loo(collect(mat_bad'), chain_index=ones(Int, size(chain_bad, 1)))
display(result_bad)[ Info: [Tutoria_outlier] observing_geometry = false (auto): worst accumulated bias 0.00535σ, 18.7× inside the 0.1σ limit — over 300 prior draws (seed 0xc70f177e5000001)
[ Info: [Tutoria_outlier] barycentric_lighttime = false (auto): changes no prediction at all — over 300 prior draws (seed 0xc70f177e5000001)
[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 11
┌ Info: Determined number type
└ T = Float64
ℓπcallback(θ): 0.000005 seconds
∇ℓπcallback(θ): 0.000036 seconds (1 allocation: 32 bytes)
┌ Info: Starting values not provided for all parameters! Guessing starting point using global optimization:
│ num_params = 11
└ num_fixed = 0
┌ Info: Found sample of initial positions
│ logpost_range = (-74.72239484285932, -65.24196408862298)
└ mean_logpost = -68.20376665661955
[ Info: Resolving chain
[ Info: Planning pointwise terms
┌ Info: Calculating pointwise likelihoods
│ samples = 1000
└ points = 8
[ Info: No source provided for samples; variables are assumed to be from a Markov Chain. If the samples are independent, specify this with keyword argument `source=:other`.
┌ Warning: Some Pareto k values are high (>.7), indicating PSIS has failed to approximate the true distribution.
└ @ ParetoSmooth ~/.julia/packages/ParetoSmooth/j8cDa/src/InternalHelpers.jl:50
Results of PSIS-LOO-CV with 1000 Monte Carlo samples and 8 data points. Total Monte Carlo SE of 0.27.
┌───────────┬────────┬──────────┬───────┬─────────┐
│ │ total │ se_total │ mean │ se_mean │
├───────────┼────────┼──────────┼───────┼─────────┤
│ cv_elpd │ -72.05 │ 13.90 │ -9.01 │ 1.74 │
│ naive_lpd │ -65.45 │ 9.78 │ -8.18 │ 1.22 │
│ p_eff │ 6.60 │ 4.15 │ 0.83 │ 0.52 │
└───────────┴────────┴──────────┴───────┴─────────┘Two things changed. p_eff — the number of parameters the data are effectively buying — has jumped, because one row is now doing work no smooth orbit can absorb; and ParetoSmooth has warned that some $\hat{k}$ crossed 0.7, which it does on its own whenever the importance-sampling approximation has failed somewhere.
To find which rows those are, threshold the vector — and, because the columns are in model order, map them back to the observation and epoch they came from with the same walk over the observations used above:
k̂ = collect(result_bad.pointwise(:pareto_k))
flagged = findall(>(0.7), k̂)
data_obs_bad = filter(!Octofitter._isprior, sys_bad.observations)
rowinst = reduce(vcat, [fill(Octofitter.likelihoodname(o), max(Octofitter._nrows(o), 1))
for o in data_obs_bad])
[(row=i, instrument=rowinst[i], epoch=epochs_bad[i], k̂=round(k̂[i], digits=2))
for i in flagged]2-element Vector{@NamedTuple{row::Int64, instrument::String, epoch::Float64, k̂::Float64}}:
(row = 6, instrument = "SPHERE", epoch = 50600.0, k̂ = 0.87)
(row = 8, instrument = "SPHERE", epoch = 50840.0, k̂ = 0.78)Now plot it. The top panel is the diagnostic — one $\hat{k}$ per data row, with the 0.7 threshold drawn and the rows above it picked out and labelled by epoch. The bottom panel is a skypanel! of the same fit, zoomed to the data, so you can see where the flagged epochs sit relative to the posterior orbits:
fig = Figure(size=(700, 900))
ax = Axis(fig[1,1], xlabel="data row", ylabel="Pareto k̂",
xticks=1:length(k̂), title="PSIS-LOO influence diagnostic")
hlines!(ax, [0.7], color=:red, linestyle=:dash)
hlines!(ax, [0.5], color=:gray, linestyle=:dot)
fine = setdiff(eachindex(k̂), flagged)
scatter!(ax, fine, k̂[fine], color=:black, markersize=12)
scatter!(ax, flagged, k̂[flagged], color=:red, marker=:diamond, markersize=17)
for i in flagged
right = i > length(k̂) / 2 # keep the label inside the axis
text!(ax, i, k̂[i];
text = "epoch $(round(Int, epochs_bad[i]))\n($(rowinst[i]))",
align = (right ? :right : :left, :center),
offset = (right ? -12 : 12, 0),
color = :red, fontsize = 12)
end
ylims!(ax, -0.1, 1.15)
rowsize!(fig.layout, 1, Relative(0.3))
# The same draws Octofitter's own figures use, restricted to the data region.
series = PosteriorSeries(model_bad, chain_bad; N=50)
skyax = skypanel!(fig[2,1], series; colorbar=false).sky
ra_all = vcat(astrom_dat_gpi.ra, astrom_dat_sphere_bad.ra)
dec_all = vcat(astrom_dat_gpi.dec, astrom_dat_sphere_bad.dec)
scatter!(skyax, ra_all[flagged], dec_all[flagged],
marker=:circle, markersize=26, strokewidth=3, strokecolor=:red,
color=(:white, 0.0))
# Equal spans, so the sky panel's DataAspect gives a square view.
ra_lo, ra_hi = extrema(ra_all)
dec_lo, dec_hi = extrema(dec_all)
half = max(ra_hi - ra_lo, dec_hi - dec_lo) / 2 + 60
cx, cy = (ra_lo + ra_hi) / 2, (dec_lo + dec_hi) / 2
limits!(skyax, cx - half, cx + half, cy - half, cy + half)
fig
The corrupted epoch is the one the top panel picks out, and the ring around it in the bottom panel shows why: it sits off the track that every posterior draw wants to follow, so the posterior is being pulled by that single row and holding it out would move the answer.
It is not usually alone. The last epoch of the campaign crosses the line here too, and that is worth understanding, because nothing was done to it: the ends of a baseline are the rows an orbit fit leans on hardest, which is already visible in the p_eff panel of the healthy fit above — rows 1 and 8 are several times any interior row. Tilting the orbit to accommodate a displaced point therefore lands on whichever rows were carrying the fit anyway. Read $\hat{k}$ as pointing at a region of the dataset that has become load-bearing, not as a list of bad measurements.
What to do about a flagged row is a judgement call, not a rule. PSIS-LOO is telling you that its approximation failed there, not that the measurement is wrong. Refit with the row genuinely held out (below) and compare; if the posterior moves materially, the measurement deserves a look before the orbit does.
Holding out whole observations
The functions below return new System objects, which you then wrap in a LogDensityModel and refit. Every one of them keeps the model's prior-shaped terms intact — dropping the UnitLengthPrior behind a UniformCircular would change the prior, not the data set.
Leave one instrument out at a time:
kfold_systems = Octofitter.generate_kfold_systems(sys)
map(s -> Octofitter.likelihoodname.(s.observations), kfold_systems)2-element Vector{Tuple{String}}:
("SPHERE",)
("GPI",)Keep only one instrument at a time — useful for checking that two instruments agree before combining them:
per_like_systems = Octofitter.generate_systems_per_like(sys)
map(s -> Octofitter.likelihoodname.(s.observations), per_like_systems)2-element Vector{Tuple{String}}:
("GPI",)
("SPHERE",)Or select by an arbitrary predicate on the observation:
filtered = Octofitter.generate_system_filtered_like(
o -> Octofitter.likelihoodname(o) == "GPI", sys)
Octofitter.likelihoodname.(filtered.observations)("GPI",)Refit any of them the usual way:
gpi_only_model = Octofitter.LogDensityModel(filtered)
gpi_only_chain = octofit(gpi_only_model, verbosity=0)[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 11
┌ Info: Determined number type
└ T = Float64
ℓπcallback(θ): 0.000009 seconds
∇ℓπcallback(θ): 0.000022 seconds (1 allocation: 32 bytes)
┌ Info: Starting values not provided for all parameters! Guessing starting point using global optimization:
│ num_params = 11
└ num_fixed = 0
┌ Warning: Verbosity toggle: unrecognized_stop_reason
│ Unrecognized stop reason: Too many steps (101) without any function evaluations (probably search has converged). Defaulting to ReturnCode.Default.
└ @ OptimizationBase ~/.julia/packages/OptimizationBase/Jfw5O/src/utils.jl:170
┌ Info: Found sample of initial positions
│ logpost_range = (-34.68634330576207, -24.256933432869086)
└ mean_logpost = -27.114133818435384Scoring the held-out dataset
Refitting without a dataset is only half the exercise; the other half is asking how well the reduced posterior predicts the data you removed. Evaluate the full model's per-row likelihoods at the reduced model's draws, and keep the columns belonging to the held-out observation:
# Which columns belong to each observation, in the order described above:
# observation by observation, prior-shaped terms skipped, and one column per
# table row (or a single column for an observation that carries no table).
data_obs = filter(!Octofitter._isprior, sys.observations)
widths = [max(Octofitter._nrows(o), 1) for o in data_obs]
bounds = cumsum(widths)
starts = [1; bounds[1:end-1] .+ 1]
heldout = [i for (o, s, e) in zip(data_obs, starts, bounds)
if Octofitter.likelihoodname(o) != "GPI" for i in s:e]
# Per-row likelihoods of *all* the data, under the GPI-only posterior
mat_all, _ = Octofitter.pointwise_like(model, gpi_only_chain)
mat_heldout = mat_all[:, heldout]
# Log pointwise predictive density of the held-out rows: log mean exp over draws,
# summed over rows.
using StatsBase: mean
lppd = sum(log(mean(exp, col)) for col in eachcol(mat_heldout))-30.856260340233575lppd is on the same scale as a log Bayes factor per held-out dataset: larger is better, and comparing it across folds says which instruments the model is and is not able to predict from the others. It is not as directly interpretable as PSIS-LOO — there is no importance-sampling diagnostic to warn you when it is unreliable — but it needs no approximation either, because each fold is a genuine refit. Note that an instrument carrying its own offset/jitter variables is being predicted with those nuisances free, so a poor score there can mean a bad calibration rather than a bad orbit.
Epoch-level folds
Data rows are numbered globally — observation by observation, then table row by table row — so row 5 means "the fifth data row in the model" regardless of which observation it lives in.
One system per data row, each containing only that row:
per_epoch_systems, per_epoch_epochs = Octofitter.generate_system_per_epoch(sys)
(length(per_epoch_systems), per_epoch_epochs)(8, [50000.0, 50120.0, 50240.0, 50360.0, 50480.0, 50600.0, 50720.0, 50840.0])One system per data row, each containing rows 1 through i — for watching a posterior tighten as data accumulate:
cumulative_systems, cumulative_epochs = Octofitter.generate_cumulative_system_per_epoch(sys)
length.(cumulative_epochs)8-element Vector{Int64}:
1
2
3
4
5
6
7
8And the general form, which takes the row groups you want:
grouped, grouped_epochs = Octofitter.generate_systems_with_epoch_groups(
sys,
[[1, 2], [3, 4, 5], [6, 7, 8]],
g -> "_group_$g",
)
grouped_epochs3-element Vector{Vector{Float64}}:
[50000.0, 50120.0]
[50240.0, 50360.0, 50480.0]
[50600.0, 50720.0, 50840.0]Observations that cannot be subset
Not every likelihood decomposes into independent per-row terms. The canonical case is MarginalizedRVObs from OctofitterRadialVelocity: it integrates the instrument's zero point out analytically, which couples every point in that instrument, so "the likelihood of row 7" does not exist as a quantity. Asking for it produces an error naming the offending observation and its type, rather than a bare failure from inside the machinery:
Octofitter.pointwise_like(model_with_marginalized_rv, chain)
# ERROR: ArgumentError: Cannot subset the data of observation "HIRES" (a MarginalizedRVObs), …Use a plain RadialVelocityObs with an explicit offset variable if you need to cross-validate radial velocities.
Holding out rows from a RadialVelocityObs that carries a Gaussian process works with both GP backends — Celerite and AbstractGPs — and does all the behind the scenes work to do cross-validation correctly in the presence of a Gaussian process: the process is conditioned on the retained rows only, and each held-out point is then scored against the predictive mean and variance there, with that point's own measurement error and jitter added.
See also
- Posterior Predictive Checks — the qualitative counterpart: does the fitted model reproduce the data at all?
- Bayesian evidence — comparing whole models rather than individual points.
prior_only_model— a copy of the model with its data likelihoods replaced by no-ops.