What I check before I believe a posterior

A fitted posterior always looks reasonable. That is the problem.

It comes out smooth, it has a sensible mode, its credible intervals are the right sort of width, and there is nothing in its appearance that distinguishes a posterior which has correctly characterised your uncertainty from one which is confidently wrong. Fluency is not correctness. If the model is misspecified, the sampler under-converged, or the prior transform quietly different from the prior you intended, you get output that looks exactly the same and means something entirely different.

So the question “should I believe this?” has to be answered with evidence, and the evidence has to be assembled deliberately, because none of it arrives on its own.

What follows is the set of checks I ended up with on a parameter-estimation pipeline in single-molecule microscopy, where the fit recovers a molecule’s position and orientation from a small patch of camera pixels. The domain is incidental. Every check below applies to any Bayesian fit where you can simulate from the model, and most of them apply where you can’t.

They are roughly in the order I learned to need them, which is also the order of increasing embarrassment.

1. The rank statistic, and why the obvious one is wrong

The standard test for calibration is simulation-based: draw parameters from the prior, simulate data, fit, and ask where the true value falls within the resulting posterior. Repeat a few hundred times. If the posteriors are calibrated, those ranks are uniform on [0,1]. A pile-up near the edges means the posteriors are too narrow; a pile-up in the middle means they are too wide.

The usual way to compute that rank is the probability integral transform — evaluate the posterior CDF at the truth. It is simple, it is fast, and it is only valid if your posterior is unimodal and your parameter lives on the real line.

Mine were neither. Azimuth is circular, so there is no CDF: the choice of where to cut the circle determines the answer, and any cut is arbitrary. Inclination is frequently multimodal at low photon counts, and a CDF rank on a bimodal posterior reports something close to 0.5 for a truth sitting in the valley between the two modes — which is exactly the failure you most want to catch, scored as a perfect result.

What works instead is a highest-density-region rank: the fraction of posterior mass sitting at a density greater than or equal to the density at the truth. That is uniform on [0,1] for a correctly calibrated posterior whatever shape it has, because it only ever asks about level sets. It handles multimodality natively — a truth in the valley scores near 1, correctly — and it works on a circle, because density has no branch cut even where the coordinate does. It piles up near 1 when the posterior is too narrow or mis-centred, which is the direction of failure that matters.

The cost is that you need a density estimate rather than an ordering, so it is noisier and it inherits every problem density estimation has. Which brings me to the part I got wrong.

2. The null you compare it against

Having computed a set of ranks, you test them against uniformity. I originally did this against a Gaussian toy problem — generate ranks from a known-good Gaussian posterior at the same sample count, and compare.

That captures binning and discretisation effects, and it is wrong in exactly the place it matters. It captures nothing about posterior shape, and the shape is most unusual near a prior bound, which is where a great deal of this posterior lives. The toy null is a good approximation precisely where I did not need one.

The fix is to generate the null from the posteriors themselves, using the same estimator: for each fit, draw a synthetic truth from its own posterior samples and compute the rank of that. Under correct calibration, observed ranks and self-null ranks are draws from the same distribution, so they can be compared with a two-sample test rather than against assumed uniformity. I use a two-sample chi-square in the Bhattacharyya form, which folds in the fact that the null is itself an estimate carrying Monte Carlo error rather than an exact expectation.

This replaced a hard-coded verdict — something like flag it if the top bin exceeds 1.3× the null — with a number and a p-value. Hard-coded thresholds in validation code are a smell. They encode a sample size and a noise level that were true once.

3. The mistake that made a correct sampler look broken

Here is the one I am most glad I found, and it took an embarrassingly long time.

If you draw a synthetic truth from a set of posterior samples and then test it against an interval or a density estimate built from that same set, the test point helped define the thing it is being tested against. The rank is biased toward the high-density end. Coverage is biased upward.

I assumed this was a second-order effect. It is not, and it scales inversely with effective sample size in a way that bites hardest exactly where samplers are already struggling.

At an effective sample size of 40, the in-sample calculation reports 71.7% achievable coverage at a nominal 68%, where the true value is 63.3%. That is an eight-percentage-point error, and it points the wrong way: a perfectly calibrated posterior appears to be under-covering by eight points, which reads as over-confidence and sends you off to debug a sampler that is working. At an effective sample size of 20 the error is fifteen points.

Excluding the test point — leave-one-out, which costs one line — tracks the true coverage to within about two points everywhere.

I had been reading a systematically wrong calibration target for weeks and treating the discrepancy as a real defect in the fit.

4. Coverage, plural, and against the right target

Two further things about coverage, both of which I initially had wrong.

One number is not enough. A posterior that is too narrow in the core and too heavy in the tails can land on 68% coverage at a nominal 68% and be badly wrong at 50% and 95%. It can hit one level by accident. It cannot sit on the curve at 50, 68, 90 and 95 by accident. So the diagnostic is a calibration curve, not a calibration number — and with Wilson intervals rather than the normal approximation, because at the 95% level the observed fraction is close to 1 and the normal approximation for a binomial proportion is at its worst precisely there.

The shape of the failure is diagnostic too. A curve sagging uniformly below expectation is straightforward over-confidence. A curve that matches at 50 and fails at 95 is a tail problem — too few live points, or wings the sampler never explored — and those have different fixes.

And the nominal level is the wrong target. The shortest-interval estimator is biased narrow at finite sample count, because it takes a minimum over many candidate windows and the minimum of a noisy quantity is smaller than the quantity. With a thousand stored samples, a perfectly calibrated posterior covers about 46% at a nominal 50%.

Read against the diagonal, that convicts a correct sampler of over-confidence. So the comparison has to be against the achievable coverage — what a known-good posterior attains with this estimator at this sample count — obtained the same way as the rank null, and with the same leave-one-out correction.

The general lesson: your estimator has finite-sample behaviour, and validating against the theoretical ideal rather than against your estimator’s own achievable performance will manufacture failures that do not exist.

5. Under-coverage is not a diagnosis

By this point I had a decent set of statistics and a real problem with them: they all live in parameter space, and they all report failure the same way.

A badly converged sampler shows up as under-coverage. A prior that does not match the truth distribution shows up as under-coverage. A misspecified forward model shows up as under-coverage. A poor point estimator shows up as under-coverage. Four different causes, four completely different fixes, one symptom.

Three checks separate them, and none is a statistic about the posterior.

Insertion order. For nested sampling there is a test — Fowlie, Handley and Su — on the order in which new live points are inserted into the sorted set. Under correct sampling from the constrained prior, that order is uniform, and it is a statement about the sampler and nothing else. A prior mismatch cannot produce a non-uniform insertion order; a misspecified model cannot either. Around 1% of runs should return p below 0.01, the median p should sit near 0.5, and the pooled p-values should pass a KS test against uniform. This is the only diagnostic I have that points at one cause unambiguously.

Replication. Run the same data with a different seed. Run-to-run scatter must be small compared with the reported posterior width. If it is not, the width is measuring the sampler rather than the posterior — and, crucially, coverage cannot detect this, because both runs are wrong in the same way and the average looks fine.

Prior recovery. Run the pipeline with the likelihood switched off. With a constant likelihood the posterior must equal the prior, exactly. This is the only check that validates the prior transform — the code that maps a unit hypercube to parameter space — as opposed to checking whether the prior was a good idea. If your unit-cube transform has a bug, every other diagnostic in this document reports something plausible and wrong.

One subtlety that cost me an afternoon: my positional prior is patch-local, while the stored coordinates are global. Pooling global positions across emitters gives a mixture over every patch centre in the field, which is uniform on nothing at all and fails the test on a perfectly correct prior. Subtract the centre first. Validation code has coordinate bugs too.

6. Leaving parameter space entirely

Everything above compares posteriors to truths. That is only possible on simulations, and it can only ever validate the pipeline against the simulator — never the simulator against reality.

For that you have to render the model back into data space and look at residuals. On a simulated set drawn from the same forward model, the Pearson residuals must be exactly noise. Structure in them means a patch-extraction, background-join, normalisation or coordinate-convention bug, and no parameter-space statistic can see any of those, however many of them you run.

Three details make this test actually work rather than merely exist.

Use the right discrepancy. Chi-square assumes Gaussian errors. In the wings of a patch the expected count is a few photons, so Gaussian is badly wrong — and the wings are exactly where a coordinate bug shows up first, because that is where the model gradient is steepest with respect to position. Poisson deviance is the natural measure for count data and does not fall over there.

Get the standard error right. For known mean, the Pearson residual has variance exactly 1 and kurtosis 3 + 1/μ, so the sample variance over N pixels has standard error √((2 + ⟨1/μ⟩)/N). The 1/μ term is not a refinement. At eight photons of background it inflates the standard error by about 6%, and dropping it turns a healthy run into a failure that is 6% too significant.

Do not use a fixed threshold. With 25,000 pixels the standard error on the residual variance is about 0.009, so a 1.5% inflation is a 1.7σ fluctuation and a 15% inflation is 17σ. The same cutoff is either far too loose or far too tight depending only on how many emitters you happened to render. Report a z-score, and report the effect size alongside it in physical units — root-mean-square photons of mis-prediction per pixel — so that the size of the problem stays legible independently of its statistical significance. A 20σ result that corresponds to a fifth of a photon is not the same thing as a 20σ result that corresponds to thirty.

And report per-emitter as well as pooled. A global z-score cannot distinguish “every patch is slightly wrong” from “two patches are catastrophically wrong”. The first is a convention or normalisation error affecting everything. The second is a handful of emitters near a frame edge or sitting on top of a neighbour. Same aggregate number, unrelated causes.

7. The bound that refuses to answer

My favourite piece of the whole suite, and the one I would port to any other project first.

The Cramér–Rao bound gives an information-theoretic floor on the variance of any unbiased estimator. Comparing your empirical scatter against it tells you what fraction of the available information you are actually extracting. It is a genuinely useful number and a very satisfying one to put in a paper.

It is also a number you can fool yourself with completely, because if your scatter comes out below the bound, that is not an excellent result. It is impossible. It means the bound is wrong.

Mine was, by a factor of 2.7, because I had computed it in global rather than patch-local coordinates and the forward model clamps at the patch edge. Had the code simply printed the ratio, I would have had a beautiful efficiency figure that was pure artefact.

So the comparison is gated. If empirical scatter beats the bound anywhere, the code prints that the bound is wrong and declines to report an efficiency at all.

I now think this generalises into a rule. Any diagnostic capable of producing a physically impossible result should detect that case and refuse to answer, rather than printing a number that happens to look good. Impossible outputs are the cheapest bug detector you will ever build, and the temptation to accept a flattering one is much stronger than you expect at the point where you are tired and the deadline is close.

There is a second signature in the same comparison worth naming, because it took me a while to interpret. If empirical scatter is far above the bound while the posterior width sits right at it, the posterior is as tight as the information allows and the errors are nonetheless large. That is not an information shortage. The information is there and is being discarded — it is the signature of mode collapse, where the fit is confidently landing on the wrong one of two modes. An information shortage and a confident wrong answer look identical if you only report one of those two ratios.

8. Small honesties in the tooling

Two things that are not statistics but belong in the same file.

Check whether your effective sample size means anything. If the stored weights are uniform — because the samples are equal-weight resampled draws — then the Kish ESS is n by construction and carries no information about the sampler whatsoever. Reporting it anyway gives you a reassuring number that is a property of your storage format. Test for weight degeneracy before quoting ESS.

Say what you did not check. My results file was resumed across a schema change, so rows fitted before the upgrade carry no sampler diagnostics. It would have been easy to report on the populated subset and let the header imply the whole set. That quietly turns a partial audit into a claimed complete one. The code now prints, loudly, which rows it could not check and states that a calibration claim over the whole set requires a refit.

A validation suite that overstates its own coverage is worse than no validation suite, because it converts an unknown into a false known.

The list

Condensed, for anyone who wants it without the narrative:

1. Rank statistic based on highest-density regions, not the CDF, if anything is multimodal or circular.
2. Compare ranks against a null generated from your own posteriors with the same estimator — not against assumed uniformity, and not against a Gaussian toy.
3. Leave the test point out. In-sample coverage and rank estimates are badly biased at low effective sample size.
4. Coverage as a curve at several levels, with Wilson intervals, against achievable rather than nominal coverage.
5. At least one diagnostic that points at the sampler alone. Insertion order for nested sampling; replicate runs otherwise.
6. A likelihood-off run, to check that the prior in the sampler is the prior you meant.
7. Residuals in data space, with the right discrepancy measure, the right standard error, no fixed threshold, and a per-unit breakdown as well as a pooled one.
8. Gate any bound that can return an impossible value.
9. Check that your ESS is not tautological.
10. State what you did not check.

None of this is exotic. Most of it is one or two lines once you know to want it. The reason it is worth writing down is that the failure mode being defended against is not a crash or an obviously silly number — it is a smooth, plausible, well-presented posterior that means something other than what you think it means, and which will pass every informal look you give it.

The validation code on that project ended up longer than the model. In hindsight that ratio was about right.