feat(eval): integrated autocorrelation time and effective sample size (T2-1a)

[Telschow]

Part of #112 (not closing it). No wiring: no verdict, table, golden file or page changes, and the estimator is not used anywhere yet. navkit run --out on origin/main and on this branch is identical after stripping runtime_s, wall_s and realtime_factor (15 cases; the cleaned JSON is the same size and compares equal).

What it adds

navkit.eval.statistics (exported from navkit.eval):

• integrated_autocorrelation_time(x, c=5.0) -> float
• effective_sample_size(x, c=5.0) -> float
• autocorrelation_estimate(x, c=5.0) -> AutocorrelationEstimate, a frozen result with tau, window, truncated, n, c, for a caller that wants to know whether the window rule was met. The two public functions return plain floats, as specified.

Method, and why Sokal windowing

Remove the mean. Autocovariance by FFT, zero-padded to a power of two of at least 2n (NumPy only). rho_k is the autocovariance at lag k over the one at lag zero (the biased estimator, sums over n, which keeps the estimated sequence positive semidefinite). tau(M) = 1 + 2 * sum_{k=1..M} rho_k. The window is the smallest M with M >= c * tau(M), searched up to n // 2. If none meets the rule, the value at n // 2 is returned and truncated is set. n_eff = n / tau, clipped to [1, n], and tau is clipped at 1 (antithetic series are not rewarded; the docstring says so).

Why windowing: a fixed window cannot work. Too short cuts off a slow tail; too long adds the noise of the sample autocorrelation at large lags, which does not shrink with the lag. The rule grows the window with the answer. Method as decided; not changed.

Refusals, each with a test: fewer than 10 samples, NaN or inf, a constant series, a non-1-D input, a non-positive c.

AR(1) spread I measured, and the tolerances I chose

tau = (1 + a) / (1 - a) exactly for x_t = a x_{t-1} + e_t. Relative deviation of the estimate from the closed form over eight seeds (100 to 107), scratch script, no series truncated:

a	n	mean deviation	sd	max abs deviation	tolerance in the test
0.5	200 000	+0.0013	0.0220	0.0378	0.08
0.9	200 000	-0.0012	0.0317	0.0523	0.12
0.99	200 000	+0.0151	0.0938	0.1558	not used (too noisy)
0.99	2 000 000	-0.0061	0.0369	0.0786	0.15
0.99	5 000 000	-0.0024	0.0334	0.0478	not used

Each tolerance is at least about 3 sd, and about twice the largest deviation seen; the bias is under 2% in all cases, so the spread, not a bias, sets the tolerance. a = 0.99 gets its own n (2 million) as asked; the test seeds (11, 12, 13) are different from the measured ones. The effective sample size is checked at 1.5 times the tolerance.

Other known-answer tests: white noise (tau within 0.05 of 1), chi-square(3) i.i.d. (a consistent NEES; within 0.06 of 1), slow AR(1) plus white noise (tau between the two pure cases and equal to 1 + (tau_AR - 1) * share within 12%, a closed form and not just an ordering), scale and offset invariance (3x + 7 and -x), tau never below 1 for an antithetic series, determinism, and that c moves the window.

Non-stationary series. A trend never meets the rule: it is cut at n // 2, marked truncated, and has a small n_eff. A random walk is the surprise: for the draws I used the rule is met late (a window of thousands of lags, a tenth of the length or more), so truncated is not set, yet tau is large and n_eff is a few samples. The test asserts the result over six draws (tau over 100, n_eff under n/100) and does not ask for a flag. Callers must not read truncated=False as proof of stationarity.

Property test: n_eff in [1, n], tau at least 1, both finite, for any AR(1) series of 10 to 400 samples with a in (-0.99, 0.99) (Hypothesis, already a dev dependency), plus a seeded loop of heavy-tailed, step, trend, periodic and random-walk series.

MEASURE (MEASUREMENT, synthetic; not field performance)

Every case in configs/benchmark.yaml through navkit.benchmark.run_case (as tests/test_observability.py does, keeping the estimator result). The series is the per-epoch NEES the calibration report scores: normalized_error_squared(est.positions, ref_on_est.positions, position_cov), 3 degrees of freedom, all 6001 epochs. The outage rows use the epochs inside the configured outage window, the half-open window_mask the report uses. The window rule cannot see lags past half the series, so in a 3000-epoch outage window it stops at 1500. trunc marks a window cut by that limit.

Case	Mean NEES	n	tau	n_eff	window	trunc
gnss_only	4.29	6001	232.0	25.9	1163	no
dead_reckoning	no covariance, not scored
vision_anchor_in_measurement_noise	164	6001	1542.3	3.9	3000	yes
vision_only	264	6001	1571.8	3.8	3000	yes
outage_control	4.13	6001	513.0	11.7	2567	no
outage_control, outage window	5.93	3000	867.0	3.5	1500	yes
outage_visual	286	6001	425.4	14.1	2128	no
outage_visual, outage window	479	3000	869.7	3.4	1500	yes
outage_visual_degraded_camera	204	6001	440.5	13.6	2204	no
outage_visual_degraded_camera, outage window	321	3000	884.5	3.4	1500	yes
outage_visual_clone	2.49	6001	413.3	14.5	2068	no
outage_visual_clone, outage window	3.13	3000	712.7	4.2	1500	yes
outage_visual_degraded_camera_rereferenced	192	6001	439.5	13.7	2198	no
outage_visual_degraded_camera_rereferenced, outage window	308	3000	896.2	3.3	1500	yes
outage_visual_degraded_camera_clone	2.18	6001	346.0	17.3	1730	no
outage_visual_degraded_camera_clone, outage window	2.66	3000	227.1	13.2	1139	no
vision_only_clone	2.29	6001	519.4	11.6	2600	no
outage_visual_clone_outliers	2.49	6001	378.9	15.8	1896	no
outage_visual_clone_outliers, outage window	2.98	3000	685.6	4.4	1500	yes
outage_visual_clone_scale_drift	2.88	6001	409.5	14.7	2049	no
outage_visual_clone_scale_drift, outage window	3.57	3000	769.9	3.9	1500	yes
outage_visual_clone_correlated	13.2	6001	379.5	15.8	1899	no
outage_visual_clone_correlated, outage window	18.5	3000	240.2	12.5	1202	no
outage_visual_clone_correlated_inflated	4.88	6001	499.0	12.0	2497	no
outage_visual_clone_correlated_inflated, outage window	7.41	3000	242.0	12.4	1213	no

Does the audit's "roughly 5 to 13 effective samples per 6001-epoch run" hold? Roughly, as an order of magnitude: single digits to low tens. Over the 14 scored whole-run series n_eff is 3.8 to 25.9; 11 of 14 fall between 11 and 18, and the two vision-only anchor rows are 3.8 and 3.9. So "5 to 13" understates the top (gnss_only is 25.9) and the bottom. In every case the run carries a few to a few tens of independent samples out of 6001, which is what matters: the verdict that treats epochs as independent is overconfident by a factor of hundreds. Caveats: the NEES series of the divergent cases (the anchor rows, the vision-only rows, the outage-window series of every outage case) are not stationary, and where trunc is yes the tau is a lower bound set by the window limit, so n_eff there is an upper bound. The window-only n_eff is also affected by there being fewer epochs than the whole run: 1500 is half of 3000.

Mutation checks (not committed)

Window always the largest allowed (window = taus.size - 1):

AssertionError: (0.5, 11)
assert 1.0 == 3.0 ± 0.24
AssertionError: (0.9, 11)
assert 5.621421133339564 == 19.000000000000004 ± 2.28
assert 40.1367054837378 == 198.99999999999983 ± 29.85

(the AR(1) known-answer tests fail for all three coefficients). Mean subtraction removed (centred = arr):

assert 44946.316451073384 == 1.0 ± 0.06           (chi-square(3) series: tau near 1)
assert 19141.077409352463 == 21.051320836838244 ± 2.1e-08   (scale and offset invariance)
Failed: DID NOT RAISE ValueError                  (the constant series is no longer refused; 3 tests)
FAILED tests/test_effective_sample_size.py::test_an_iid_chi_square_series_has_tau_near_one

Both restored; all tests pass on the unmutated code.

Checks

ruff check, ruff format --check, mypy src, metrics.py --check (test count regenerated, 1629), check_doc_tables.py and check_links.py pass. The new file runs in about 1 s. Full suite: 1598 passed, 29 skipped, 2 xfailed in 636.10s (0:10:36) (the skips are the dataset-dependent tests; the datasets are absent in this clone). The coverage ratchet is left to CI.

🤖 Generated with Claude Code