Uncertainties

How much can you trust a fitted parameter? OITOOLS offers two answers, and they do not measure the same thing:

This page reports a simulation study that measures how well each of them actually performs, and compares them with PMOIRED on identical data. The code is in demos/bootstrap_validation; RESULTS.md there carries the same tables.

How the calibration is measured

An uncertainty is a claim about a distribution you never see: the spread of best-fit values you would get if you repeated the observation. So the study repeats it. 120 independent noise realisations ("universes") of the same observation are generated, each is fitted, and what each method claims about a single universe is compared with what the ensemble actually does:

sigma_true    = std over universes of the best-fit value      <- the truth
sigma_report  = mean over universes of the quoted 1 sigma     <- the claim
ratio         = sigma_report / sigma_true                     <- 1.00 = calibrated
coverage      = fraction of universes with |best - truth| < sigma_report
                                                              <- 0.683 = calibrated

A ratio below 1 means overconfident error bars. sigma_true is itself estimated from 120 universes, so it carries about 6.5% noise — but that noise is common to every method within a regime, so the differences between methods in a column are far better determined than 6.5%.

The model is a binary at CHARA/MIRC H-band resolution (a partially resolved primary plus an unresolved companion at 1.44 mas with 5% of the flux), fitted for four parameters against V² and T3PHI. The uv coverage, MJDs, telescope configurations and error bars come from a real dataset; only the observable values are replaced.

Five regimes stress different assumptions:

Choosing a resampling scheme

bootstrap_fit(...; mode=...) selects how block multiplicities are drawn. Wall times are per replicate, single-threaded, on the benchmark dataset (1440 data points, 180 blocks, 4 free parameters); reproduce with demos/bootstrap_validation/benchmark.jl.

The block columns give the calibration ratio; bold marks the best in each.

• :halfsample is the best-calibrated scheme wherever blocks are numerous, and the cheapest, because each replicate fits half the data. Its failure mode with few blocks is over-conservatism, which is the safe direction.
• :replacement is the default: the textbook bootstrap, the most uniform across regimes, never off by more than 14% including the single-snapshot case.
• :weights is consistently the most optimistic of the three, and clearly so when blocks are scarce. It is included for completeness.
• :pmoired exists only to reproduce PMOIRED results, and is biased low by about √2 (see the warning above and the tables below).

A replicate costs about the same as the original fit: 3.4 ms against 3.8 ms for a single Levenberg-Marquardt fit, of which only 0.85 ms is the resampling — the rest is the fit itself, so the scheme you choose barely matters for runtime. bootstrap_fit spreads replicates over all available Julia threads; 500 replicates of this problem take about 0.6 s on 16 threads.

Granularity matters as much as the scheme. :config (the default, one block per MJD and telescope configuration) works from 18 blocks upward; :epoch becomes erratic when there are few epochs (0.92–1.83 on a single snapshot); :point destroys the correlation structure and underestimates by design.

A note on MJD precision

Blocking is keyed on the MJD, and Float32 cannot resolve it: near MJD 55000 its representable values are 3.9 × 10⁻³ d apart — 5.6 minutes — so exposures taken minutes apart would collapse onto one value and their blocks would silently merge.

readoifits therefore stores v2_mjd, t3_mjd, vis_mjd and flux_mjd as Float64 regardless of T, so the block structure is identical whether you read with the default T=Float32 or with T=Float64. (uv_mjd follows T — it is the array handed to the model evaluator as $MJD, where the extra precision buys nothing.) Nothing to do on your side; the numbers below are unaffected by the precision you choose.

Results

ratio = quoted sigma / true scatter, one column per parameter; the last column is the coverage averaged over the four parameters (0.683 when calibrated).

Regime: ideal

Regime: systematics

Regime: underestimated

Regime: fewblocks

Regime: mismatch

What the numbers say

1. The two packages' analytic errors are the same number. OITOOLS analytic (LM) and PMOIRED analytic (dpfit) agree to two decimals in every regime. Both rescale the covariance to χ²ᵣ = 1, so the two implementations are cross-validated against each other.

2. PMOIRED's bootstrap is overconfident by a factor 1/√2, in every regime. PMOIRED bootstrapFit lands at 0.64–0.74 everywhere, and OITOOLS' emulation of the same scheme (mode=:pmoired) reproduces it to ±0.02 in all 20 (regime, parameter) cells. The deficit is the scheme — two independent half-samples give block multiplicities with variance ½ instead of 1 — not an implementation detail. Coverage is ~0.50 instead of 0.683.

3. Only resampling sees correlated systematics. With a 5% per-(epoch, baseline) calibration error the analytic covariance is 3–5× too small even after its χ²ᵣ rescaling, which is itself a factor 3.8 here.

4. χ²ᵣ rescaling handles pure underestimation. With noise at 2× the quoted errors, both analytic estimates come back calibrated (0.95–1.01).

5. Adding noise to the data is not a bootstrap. Refitting perturb_data replicates ("parametric Monte Carlo") is dominated everywhere: identical to the analytic covariance in the well-behaved regimes, exactly half in underestimated — the factor it was lied to about — and 11–25× too small under systematics, because by construction it believes the error bars. This is why it is not offered as an uncertainty API.

6. No uncertainty is meaningful when the model is wrong. With an unmodelled third component the best fit is biased by 10–36 true σ; the bootstraps inflate their error bars by 5–11× (they do sense the unmodelled structure) and coverage still collapses to 0.22–0.25, the analytic errors to 0.05. Resampling measures precision, never accuracy — a small error bar on a bad model is still a bad answer.

Cross-package performance

Same file, same model, same number of replicates:

PMOIRED's replicate costs more than its own single fit because each replicate fits a dataset built from two copies; its multi-process scaling is limited by process start-up and by pickling the data to the workers.

Reproducing

cd demos/bootstrap_validation
OUT=/tmp/bootstrap_validation

julia --project=../.. make_universes.jl \$OUT 120   # simulate
python3 sanitize.py \$OUT                           # make files readable by PMOIRED
julia --project=../.. -t auto run_oitools.jl \$OUT 100
PMOIRED_MAX_THREADS=12 python3 run_pmoired.py \$OUT 100
python3 analyze.py \$OUT                            # the tables above
julia --project=../.. benchmark.jl \$OUT            # the wall times above

See demos/bootstrap_validation/README.md for the details, including the two interoperability fixes needed to make PMOIRED read OITOOLS-written OIFITS files.