Recompute the ringdown posterior calibration

The derivations behind each test are in the course chapters on credible regions, rank tests and model mismatch. The calibration recipe follows Talts and colleagues, and the importance-sampling check follows Dax and colleagues. Talts and colleagues, 2018 Dax and colleagues, 2023

Get the files and the environment

Start with the README. It lists the model, run order, timings and fingerprints. The scripts sit in one folder, with settings in inputs/ and results in outputs/. This command recreates that layout, using curl's own list syntax:

mkdir ac87 && cd ac87
curl -fsS --create-dirs -o "#1" "https://sandboxphysics.com/stories/simulation-based-inference/materials/\
{README.md,requirements.txt,ringdown_model.py,neural_posterior.py,calibration_statistics.py,\
calibrate_linear_gaussian.py,validate_grid_posterior.py,train_neural_posterior.py,validate_neural_posterior.py,\
inputs/ringdown-config.json,inputs/training-settings.json,inputs/linear-gaussian-settings.json,\
inputs/grid-validation-settings.json,inputs/neural-validation-settings.json,\
outputs/frozen-neural-posterior.json,outputs/training-log.json,outputs/failure-cases.json}"
python3 -m pip install -r requirements.txt

The only dependency is NumPy, pinned to 2.0.2 in requirements.txt. The published runs used Python 3.9.6 on an arm64 macOS CPU. Keep the saved outputs you want to compare against, such as the conjugate Gaussian results, in a separate folder, because each script overwrites its own output.

Units and conventions

The settings files in the model configuration also fix the three mismatch families. Their ratios are teaching choices. The second-mode ratio is a teaching value, and the glitch is a sine-Gaussian wave packet standing in for a detector transient.

Step 1: calibrate a model whose answer is known

python3 calibrate_linear_gaussian.py
shasum -a 256 outputs/linear-gaussian-sbc.json

This runs in about 2 s. The prior is $\theta\sim\mathcal N(0,1)$ and four measurements have unit Gaussian noise, so the exact posterior is normal with variance $0.2$. The script tests that posterior and four wrong ones over $4000$ simulations with $99$ draws each, binned into 20 bars. The fingerprint should be 8bcd4246…cebd4, the value in the README.

Open the output and read three fields for each algorithm. counts is the histogram; compare it with uniform99Band, which is $[165, 236]$. analyticExpectedCounts is the histogram predicted without simulation from the known form of each wrong posterior. centralIntervalCoverage gives the simulated fraction of truths inside the 50% and 90% intervals, a Wilson 95% interval for that fraction, and the analytic value.

The exact posterior in this run fails the uniformity test, with a tail probability of $0.0093$. This is the first entry in the failure log. repeatStudy explains it: 200 independent repeats flag the exact posterior twice and every wrong posterior 200 times. A single run at a 1% threshold will raise false alarms about once in a hundred.

Step 2: check the exact grid posterior

python3 validate_grid_posterior.py

This takes about 15 s and writes grid-validation.json. The grid posterior is the reference for everything after it, so it gets four checks that do not reuse its code path.

Each ratio is the grid posterior standard deviation divided by the Fisher prediction. The Fisher matrix is a high signal-to-noise limit, so agreement is expected only in the bottom rows. At $\sigma=0.6$ the posterior is not Gaussian and is cut by the prior box.

A Metropolis sampler run on the same likelihood gives a frequency mean of $249.993\,\mathrm{Hz}$ and standard deviation $3.32\,\mathrm{Hz}$, against $249.994$ and $3.28$ from the grid. SBC of the grid over 2000 simulations gives rank tail probabilities of $0.336$ and $0.071$ and 90% coverage of $0.897$. Finally the first-order bias formula for a small unmodelled mode agrees with the grid shift to better than 1% at strength $0.02$.

Step 3: the frozen network and its training cost

The frozen network is a mixture of five Gaussians in logit-transformed coordinates. Its inputs are 20 projections of the record onto a template basis plus $\ln(\sigma/0.3)$. Two hidden layers of 128 tanh units give $23{,}198$ parameters. Its fingerprint is bb213976…ed575.

Training is optional:

python3 train_neural_posterior.py

It took $413.8\,\mathrm{s}$ on one CPU and saw $30.72$ million simulations, each drawn fresh from the well-specified model. The final validation loss in the training log is $1.512$. No mismatch family is ever shown to the network. Retraining from seed 87 should give a network with the same behaviour; identical bytes are expected only on the same NumPy build.

Step 4: test the network

python3 validate_neural_posterior.py

This takes about 70 s and writes neural-validation.json and browser-parity-cases.json. The second file fixes inputs and outputs that the Lab's browser code must reproduce. Its fingerprint should be 8f53fbf8…47ba0.

Read the validation file in this order.

1. sbc: rank histograms over 2000 records drawn from the prior. The tail probabilities are $0.462$ for frequency and $0.800$ for damping time.
2. wellSpecifiedSummary: coverage at levels from $0.05$ to $0.95$ for network and grid, median importance-sampling efficiency $0.966$, median Jensen-Shannon divergence to the grid $0.0045$ nats, and residual tests flagged in $0.9\%$ of records.
3. conditionalAndHeldOut: nine groups of 400 records. Three noise bands and short damping times stay inside training. Two noise levels lie outside it, and three mismatch families are never seen in training.
4. labPresetsAtFixedSource: the six Lab presets, 400 records each, which the Lab draws as dotted reference marks.

Each coverage carries a Wilson 95% interval. With 400 records that interval is about $\pm0.03$ near $0.9$, so only differences larger than that count as evidence.

Failures kept in the record

The failure log keeps six entries. Two matter most for a reader repeating this work.

The first is a sampler bug. The network's sampler first used the wrong correlation term, yet one-dimensional rank tests passed with tail probabilities $0.462$ and $0.781$. A comparison between sampler moments and the density's own moments caught it, with an error of $0.223$ that fell to $0.0038$ after the fix. Marginal SBC cannot see a wrong correlation.

The second is a limitation that no check here detects. An unmodelled second mode of strength $0.6$ lowers both posteriors' 90% coverage to about $0.74$, while importance-sampling efficiency stays near $0.98$ and the residual test flags about 6% of records. The other entries record a coarse quadrature that mispredicted rank histograms, the low signal-to-noise failure of the Fisher limit, and spurious matrix warnings from one NumPy build.

Checklist for your own simulation-based analysis

• Calibrate first on a model whose posterior you can write down, and predict the wrong-posterior histograms before simulating them.
• Record seeds, versions, timings and fingerprints, and keep a failed run instead of reseeding it.
• Check the reference posterior itself against a limit and an independent sampler.
• Test joint coverage and the sampler against the density, not only one-dimensional ranks.
• Split coverage by region of parameter space and by noise level.
• Run held-out and mismatched records, and report which check flags each one.