test(polytomous): require bias MAE and uncertainty coverage (#1343)

* test(grm): require bias MAE and coverage recovery

* test(gpcm): require bias MAE and coverage recovery

* test(fipc): require bias MAE and coverage recovery

* docs(changelog): record polytomous recovery calibration evidence

* test(grm): label posterior-SD coverage precisely

* test(gpcm): label posterior-SD coverage precisely

* test(fipc): label posterior-SD coverage precisely

* docs(changelog): clarify recovery coverage approximation

* test(cat): require bias MAE and uncertainty coverage

* docs(test): record CAT uncertainty recovery evidence

Author: seonghobae

docs/changelog.d/1342-polytomous-recovery-calibration.md

@@ -0,0 +1,7 @@
+# Strengthen polytomous recovery calibration evidence
+
+## Changed
+
+- Extend deterministic GRM, GPCM, CAT, and fixed-item-parameter recovery studies to require signed bias, MAE, finite positive Rust-returned posterior uncertainty, and empirical coverage of the normal-approximation interval `theta_eap ± 1.96 * theta_sd` alongside RMSE, while retaining correlation only as supplementary recovery evidence.
+- Preserve CAT's independent adaptive-efficiency gate on mean administered items, so uncertainty calibration and error recovery cannot mask a fallback to non-adaptive item selection.
+- Keep all production likelihood, marginal-ML/EM, EAP/CAT scoring, item-information/selection, stopping, and uncertainty arithmetic Rust-owned; the added Python calculations are explicit true-parameter recovery-test summaries only.

tests/test_polytomous_cat_recovery.py

@@ -1,16 +1,21 @@
-"""Accuracy-and-efficiency recovery study for public polytomous CAT.
+"""Accuracy, uncertainty-calibration, and efficiency recovery for polytomous CAT.
 
 A GRM bank is calibrated from synthetic responses and then passed to
 ``cat_simulate_polytomous`` for examinees with known true abilities. The test
-requires both acceptable theta recovery and materially fewer administered
-items than the full bank, so a silent fallback to non-adaptive item order can
-fail even when correlation remains acceptable.
+requires bounded bias/MAE/RMSE, calibrated Rust-returned posterior uncertainty,
+and materially fewer administered items than the full bank. Correlation remains
+supplementary: a highly correlated CAT can still be biased or understate score
+uncertainty, while a silent fallback to non-adaptive item order can still meet
+score-recovery bounds.
 
 Dodd, De Ayala, and Koch (1995) evaluate computerized adaptive testing with
 polytomous IRT items and motivate measuring adaptive efficiency together with
 score quality. This regression uses the package's Rust-owned item-information,
-selection, scoring, and stopping path; the Python formula below only generates
-synthetic GRM responses.
+selection, scoring, posterior-SD, and stopping path; Python only generates the
+synthetic GRM responses and summarizes recovery against known simulation truth.
+The interval ``theta_eap ± 1.96 * theta_sd`` is an explicit normal approximation
+based on the returned posterior mean and posterior SD, not a claim of an exact
+posterior credible interval.
 
 Reference
 ---------
@@ -37,15 +42,19 @@
 # correlation ~0.91 using a mean of ~8.7 of 40 items -- comparable accuracy
 # to the full-bank GRM recovery test (~0.38/~0.92) at roughly a fifth of
 # the items. Margins are loose enough to tolerate a minor fast-mlsirm
-# version bump while still catching a real regression in either accuracy
-# or the adaptive-selection efficiency CAT exists to provide.
+# version bump while still catching a real regression in recovery,
+# uncertainty calibration, or the adaptive-selection efficiency CAT exists
+# to provide.
 #
 # MAX_MEAN_ITEMS_USED is deliberately close to the measured ~8.7 (not a
 # loose N_ITEMS * 0.5): the same fixture/seed with adaptive=False (random
 # item order) measures mean_items_used ~14.97, which still clears rmse/
 # correlation bounds -- so a bound of 12 is what actually catches a silent
 # fallback to non-adaptive selection.
 MAX_THETA_RMSE = 0.65
+MAX_THETA_MAE = 0.52
+MAX_ABS_THETA_BIAS = 0.20
+MIN_NORMAL_APPROX_COVERAGE = 0.80
 MIN_THETA_CORRELATION = 0.7
 MAX_MEAN_ITEMS_USED = 12
 
@@ -95,15 +104,35 @@ def test_cat_recovers_theta_using_substantially_fewer_items_than_the_full_bank()
         seed=SEED,
     )
     theta_eap = cat_result["theta_eap"]
+    theta_sd = cat_result["theta_sd"]
     n_used = cat_result["n_used"]
 
-    rmse = float(np.sqrt(np.mean((theta_eap - true_theta) ** 2)))
+    assert np.all(np.isfinite(theta_sd))
+    assert np.all(theta_sd > 0.0)
+
+    error = theta_eap - true_theta
+    bias = float(np.mean(error))
+    mae = float(np.mean(np.abs(error)))
+    rmse = float(np.sqrt(np.mean(error**2)))
+    normal_approx_coverage = float(
+        np.mean(np.abs(error) <= 1.96 * theta_sd)
+    )
     correlation = float(np.corrcoef(theta_eap, true_theta)[0, 1])
     mean_items_used = float(n_used.mean())
 
+    assert abs(bias) < MAX_ABS_THETA_BIAS, (
+        f"CAT theta bias {bias:.3f} exceeded ±{MAX_ABS_THETA_BIAS}"
+    )
+    assert mae < MAX_THETA_MAE, (
+        f"CAT theta MAE {mae:.3f} exceeded {MAX_THETA_MAE}"
+    )
     assert rmse < MAX_THETA_RMSE, (
         f"CAT theta RMSE {rmse:.3f} exceeded {MAX_THETA_RMSE}"
     )
+    assert normal_approx_coverage >= MIN_NORMAL_APPROX_COVERAGE, (
+        "CAT theta normal-approximation coverage "
+        f"{normal_approx_coverage:.3f} below {MIN_NORMAL_APPROX_COVERAGE}"
+    )
     assert correlation > MIN_THETA_CORRELATION, (
         f"CAT theta correlation {correlation:.3f} below {MIN_THETA_CORRELATION}"
     )

tests/test_polytomous_fipc_recovery.py

@@ -36,10 +36,13 @@
 # RMSE ~0.33 and correlation ~0.94, with the FIPC-estimated period-2 mean
 # landing at ~0.37 against a true realized mean of ~0.43 (the shift is
 # substantially recovered, not perfectly -- sampling noise on 250 people).
-# An independent free refit of the same period-2 data, by contrast,
-# re-centers to an estimated mean of ~0.01 -- essentially erasing the
-# shift, exactly the failure mode fixed-bank scoring is intended to expose.
+# Bias/MAE and normal-approximation coverage based on the Rust-returned
+# posterior SD are required independently of correlation so fixed-bank
+# recovery cannot pass on rank ordering alone.
 MAX_THETA_RMSE = 0.6
+MAX_THETA_MAE = 0.5
+MAX_ABS_THETA_BIAS = 0.20
+MIN_THETA_NORMAL_APPROX_COVERAGE = 0.80
 MIN_THETA_CORRELATION = 0.75
 MIN_FIPC_MEAN_SHIFT_DETECTED = PERIOD_2_TRUE_MEAN_SHIFT * 0.5
 MAX_INDEPENDENT_REFIT_MEAN_SHIFT_DETECTED = PERIOD_2_TRUE_MEAN_SHIFT * 0.35
@@ -107,14 +110,34 @@ def test_fipc_recovers_period_two_mean_shift_that_an_independent_refit_would_hid
     # Fixed-bank scoring: period 2 uses period 1's item parameters unchanged.
     fipc_scored = score_polytomous(responses_period_2, period_1_fit)
     fipc_theta_eap = fipc_scored["theta_eap"]
-
-    rmse = float(np.sqrt(np.mean((fipc_theta_eap - theta_period_2) ** 2)))
+    fipc_theta_sd = fipc_scored["theta_sd"]
+
+    error = fipc_theta_eap - theta_period_2
+    bias = float(np.mean(error))
+    mae = float(np.mean(np.abs(error)))
+    rmse = float(np.sqrt(np.mean(error**2)))
+    normal_approx_coverage = float(
+        np.mean(np.abs(error) <= 1.96 * fipc_theta_sd)
+    )
     correlation = float(np.corrcoef(fipc_theta_eap, theta_period_2)[0, 1])
     fipc_detected_shift = float(fipc_theta_eap.mean())
 
+    assert np.all(np.isfinite(fipc_theta_sd))
+    assert np.all(fipc_theta_sd > 0.0)
+    assert abs(bias) < MAX_ABS_THETA_BIAS, (
+        f"FIPC period-2 theta absolute bias {abs(bias):.3f} exceeded "
+        f"{MAX_ABS_THETA_BIAS}"
+    )
+    assert mae < MAX_THETA_MAE, (
+        f"FIPC period-2 theta MAE {mae:.3f} exceeded {MAX_THETA_MAE}"
+    )
     assert rmse < MAX_THETA_RMSE, (
         f"FIPC period-2 theta RMSE {rmse:.3f} exceeded {MAX_THETA_RMSE}"
     )
+    assert normal_approx_coverage >= MIN_THETA_NORMAL_APPROX_COVERAGE, (
+        f"FIPC period-2 theta mean±1.96 posterior-SD coverage "
+        f"{normal_approx_coverage:.3f} below {MIN_THETA_NORMAL_APPROX_COVERAGE}"
+    )
     assert correlation > MIN_THETA_CORRELATION, (
         f"FIPC period-2 theta correlation {correlation:.3f} "
         f"below {MIN_THETA_CORRELATION}"

tests/test_polytomous_gpcm_recovery.py

@@ -27,11 +27,13 @@
 SEED = 20260101
 
 # A real run with these exact parameters/seed measures RMSE ~0.30 and
-# correlation ~0.95 -- stronger recovery than the GRM test's ~0.38/~0.92
-# at the same sample size. The margins below stay loose enough to tolerate
-# a minor fast-mlsirm version bump while still catching a real recovery
-# regression.
+# correlation ~0.95. Bias/MAE and normal-approximation coverage based on the
+# Rust-returned posterior SD are asserted separately so correlation remains
+# supplementary rather than standing in for calibration.
 MAX_THETA_RMSE = 0.55
+MAX_THETA_MAE = 0.45
+MAX_ABS_THETA_BIAS = 0.15
+MIN_THETA_NORMAL_APPROX_COVERAGE = 0.80
 MIN_THETA_CORRELATION = 0.8
 
 
@@ -65,13 +67,30 @@ def test_gpcm_recovers_true_theta_within_expected_rmse() -> None:
 
     scored = score_polytomous(responses, fit)
     theta_eap = scored["theta_eap"]
+    theta_sd = scored["theta_sd"]
 
-    rmse = float(np.sqrt(np.mean((theta_eap - true_theta) ** 2)))
+    error = theta_eap - true_theta
+    bias = float(np.mean(error))
+    mae = float(np.mean(np.abs(error)))
+    rmse = float(np.sqrt(np.mean(error**2)))
+    normal_approx_coverage = float(np.mean(np.abs(error) <= 1.96 * theta_sd))
     correlation = float(np.corrcoef(theta_eap, true_theta)[0, 1])
 
+    assert np.all(np.isfinite(theta_sd))
+    assert np.all(theta_sd > 0.0)
+    assert abs(bias) < MAX_ABS_THETA_BIAS, (
+        f"theta recovery absolute bias {abs(bias):.3f} exceeded {MAX_ABS_THETA_BIAS}"
+    )
+    assert mae < MAX_THETA_MAE, (
+        f"theta recovery MAE {mae:.3f} exceeded {MAX_THETA_MAE}"
+    )
     assert rmse < MAX_THETA_RMSE, (
         f"theta recovery RMSE {rmse:.3f} exceeded {MAX_THETA_RMSE}"
     )
+    assert normal_approx_coverage >= MIN_THETA_NORMAL_APPROX_COVERAGE, (
+        f"theta mean±1.96 posterior-SD coverage {normal_approx_coverage:.3f} below "
+        f"{MIN_THETA_NORMAL_APPROX_COVERAGE}"
+    )
     assert correlation > MIN_THETA_CORRELATION, (
         f"theta recovery correlation {correlation:.3f} below {MIN_THETA_CORRELATION}"
     )

tests/test_polytomous_grm_recovery.py

@@ -29,11 +29,13 @@
 
 # Reasonable bounds for a 12-item, 4-category GRM test at this sample size:
 # a real run with these exact parameters/seed measures RMSE ~0.38 and
-# correlation ~0.92, comfortably inside literature-typical recovery for a
-# test this length. The margins below are loose enough to tolerate a minor
-# fast-mlsirm version bump while still catching an actual estimation
-# regression (e.g. RMSE blowing up past ~1 std or correlation collapsing).
+# correlation ~0.92. Bias/MAE and normal-approximation coverage based on the
+# Rust-returned posterior SD are asserted separately so correlation remains
+# supplementary rather than standing in for calibration.
 MAX_THETA_RMSE = 0.6
+MAX_THETA_MAE = 0.5
+MAX_ABS_THETA_BIAS = 0.15
+MIN_THETA_NORMAL_APPROX_COVERAGE = 0.80
 MIN_THETA_CORRELATION = 0.75
 
 
@@ -85,13 +87,30 @@ def test_grm_recovers_true_theta_within_expected_rmse() -> None:
 
     scored = score_polytomous(responses, fit)
     theta_eap = scored["theta_eap"]
+    theta_sd = scored["theta_sd"]
 
-    rmse = float(np.sqrt(np.mean((theta_eap - true_theta) ** 2)))
+    error = theta_eap - true_theta
+    bias = float(np.mean(error))
+    mae = float(np.mean(np.abs(error)))
+    rmse = float(np.sqrt(np.mean(error**2)))
+    normal_approx_coverage = float(np.mean(np.abs(error) <= 1.96 * theta_sd))
     correlation = float(np.corrcoef(theta_eap, true_theta)[0, 1])
 
+    assert np.all(np.isfinite(theta_sd))
+    assert np.all(theta_sd > 0.0)
+    assert abs(bias) < MAX_ABS_THETA_BIAS, (
+        f"theta recovery absolute bias {abs(bias):.3f} exceeded {MAX_ABS_THETA_BIAS}"
+    )
+    assert mae < MAX_THETA_MAE, (
+        f"theta recovery MAE {mae:.3f} exceeded {MAX_THETA_MAE}"
+    )
     assert rmse < MAX_THETA_RMSE, (
         f"theta recovery RMSE {rmse:.3f} exceeded {MAX_THETA_RMSE}"
     )
+    assert normal_approx_coverage >= MIN_THETA_NORMAL_APPROX_COVERAGE, (
+        f"theta mean±1.96 posterior-SD coverage {normal_approx_coverage:.3f} below "
+        f"{MIN_THETA_NORMAL_APPROX_COVERAGE}"
+    )
     assert correlation > MIN_THETA_CORRELATION, (
         f"theta recovery correlation {correlation:.3f} below {MIN_THETA_CORRELATION}"
     )