From 098fbb692a004096afa3fc8682c6da30a1c92c28 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 26 Jul 2026 00:19:09 +0900 Subject: [PATCH 1/2] Add confidence-interval (ACI) classification for CAT Implement mlsirm_core::exposure::ci_classify: single-cut binary-response classification by interim EAP (fixed 41-point [-4,4] grid, standard-normal log prior) with SE = EAP posterior SD; interval theta_hat +/- z_crit * se against theta_cut with STRICT first-crossing decisions and full theta/se/lower/upper counterfactual traces. Verified against R catIrt termCI.R/eapEst.R/catIrt.Rd at commit c9e979e4812c27d95d367a7f097edfe8e93ac8eb (READ). Kingsbury & Weiss (1983), Thompson (2007), and Eggen & Straetmans (2000) NOT method-section verified; cited as historical/background only. Tests: pinned 17-digit independent-Python oracle (all four traces at 1e-12), below/continue paths, full error paths, MC-500 #[ignore] structural invariants (passed). Mutation kills executed: M1 swapped decisions (2 tests fail), M2 point-estimate-vs-cut (2 fail), M3 variance-instead-of-SD (1 fail), M4 n_used off-by-one (2 fail); original restored green. PyO3 py_ci_classify + Python wrapper with validate-before-cast (complex rejection, exact 0/1 responses); pytest TestCiClassify 5 passed. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- CHANGELOG.md | 13 ++ crates/fast-mlsirm-py/src/lib.rs | 63 ++++++++- crates/mlsirm-core/src/exposure.rs | 180 ++++++++++++++++++++++++ python/fast_mlsirm/__init__.py | 2 + python/fast_mlsirm/exposure.py | 98 +++++++++++++ tests/test_paper_features.py | 108 ++++++++++++++ tests/unit/exposure_tests.rs | 218 ++++++++++++++++++++++++++++- 7 files changed, 676 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a31b799e5..6ae203f5e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -119,6 +119,19 @@ ### Added +- **Confidence-interval (ACI) classification for CAT** + (`fast_mlsirm.ci_classify`; in Rust `mlsirm_core::exposure::ci_classify`, + PyO3 `py_ci_classify`): single-cut binary-response classification by + interim EAP ability estimate on a fixed 41-point `[-4, 4]` grid with + standard-normal prior, SE = EAP posterior SD, interval + `theta_hat +/- z_crit * se` vs `theta_cut` with STRICT first-crossing + decisions -> `"above"`/`"below"`/`"continue"` with 1-based `n_used`; full + theta/se/lower/upper traces are returned as offline diagnostics (entries + past `n_used` are counterfactual replay values). Verified against R catIrt + `termCI.R`/`eapEst.R`/`catIrt.Rd` at commit + `c9e979e4812c27d95d367a7f097edfe8e93ac8eb` (READ); Kingsbury & Weiss + (1983), Thompson (2007), and Eggen & Straetmans (2000) were NOT + method-section verified and are historical/background context only. - **Wald SPRT classification for CAT** (`fast_mlsirm.sprt_classify`; in `mlsirm_core::exposure`). Single-cut binary-response sequential probability ratio test: point hypotheses at `theta_cut -/+ delta`, cumulative binary diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 4e11a2f23..e88890802 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -40,10 +40,10 @@ use mlsirm_core::dif::{ }; use mlsirm_core::exposure::{ a_stratified as core_a_stratified, ccat_select as core_ccat_select, - epv_select as core_epv_select, kl_information as core_kl_information, - kl_select as core_kl_select, owen_cat as core_owen_cat, owen_update as core_owen_update, - sprt_classify as core_sprt_classify, sympson_hetter as core_sympson_hetter, AStratifiedConfig, - SympsonHetterConfig, + ci_classify as core_ci_classify, epv_select as core_epv_select, + kl_information as core_kl_information, kl_select as core_kl_select, owen_cat as core_owen_cat, + owen_update as core_owen_update, sprt_classify as core_sprt_classify, + sympson_hetter as core_sympson_hetter, AStratifiedConfig, SympsonHetterConfig, }; use mlsirm_core::facets::fit_facets as core_fit_facets; use mlsirm_core::factor::{ @@ -2482,6 +2482,60 @@ fn py_sprt_classify( Ok(out.into()) } +/// Single-cut binary-response confidence-interval (ACI) classification for +/// CAT (`mlsirm_core::exposure::ci_classify`). Interim EAP on a fixed +/// 41-point [-4, 4] grid with standard-normal prior; SE is the EAP posterior +/// SD; interval `theta_hat +/- z_crit * se` vs `theta_cut` with STRICT +/// first-crossing decisions ("above"/"below"/"continue"). Trace entries past +/// `n_used` are offline counterfactual replay values. +/// +/// References (APA 7th; see the core module comment for read/not-read +/// source status): +/// Nydick, S. W. (2014). catIrt (R package). (READ: termCI.R, eapEst.R, +/// catIrt.Rd at commit c9e979e4812c27d95d367a7f097edfe8e93ac8eb) +/// Kingsbury, G. G., & Weiss, D. J. (1983). In D. J. Weiss (Ed.), New +/// horizons in testing (pp. 257-283). Academic Press. (NOT read; +/// historical origin) +/// Thompson, N. A. (2007). Practical Assessment, Research & Evaluation, +/// 12(1). (NOT read for the CI method section; background only) +#[pyfunction] +fn py_ci_classify( + py: Python<'_>, + a: PyReadonlyArray1<'_, f64>, + b: PyReadonlyArray1<'_, f64>, + c: PyReadonlyArray1<'_, f64>, + responses: PyReadonlyArray1<'_, u8>, + theta_cut: f64, + z_crit: f64, +) -> PyResult> { + let res = core_ci_classify( + a.as_slice()?, + b.as_slice()?, + c.as_slice()?, + responses.as_slice()?, + theta_cut, + z_crit, + ) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("decision", res.decision)?; + out.set_item("n_used", res.n_used)?; + out.set_item( + "theta_trace", + numpy::PyArray1::from_slice(py, &res.theta_trace), + )?; + out.set_item("se_trace", numpy::PyArray1::from_slice(py, &res.se_trace))?; + out.set_item( + "lower_trace", + numpy::PyArray1::from_slice(py, &res.lower_trace), + )?; + out.set_item( + "upper_trace", + numpy::PyArray1::from_slice(py, &res.upper_trace), + )?; + Ok(out.into()) +} + /// Horn's parallel analysis for principal-component retention /// (`mlsirm_core::parallel`; oracle: CRAN paran 1.5.6, PCA path). `data` is /// a flattened row-major `n_persons * n_items` matrix; `centile` is 0 for @@ -6303,6 +6357,7 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(py_ccat_select, m)?)?; m.add_function(wrap_pyfunction!(py_epv_select, m)?)?; m.add_function(wrap_pyfunction!(py_sprt_classify, m)?)?; + m.add_function(wrap_pyfunction!(py_ci_classify, m)?)?; m.add_function(wrap_pyfunction!(guttman_lambdas, m)?)?; m.add_function(wrap_pyfunction!(tenberge_mu, m)?)?; m.add_function(wrap_pyfunction!(cronbach_alpha, m)?)?; diff --git a/crates/mlsirm-core/src/exposure.rs b/crates/mlsirm-core/src/exposure.rs index e675e37a3..bd9ca6c61 100644 --- a/crates/mlsirm-core/src/exposure.rs +++ b/crates/mlsirm-core/src/exposure.rs @@ -1602,3 +1602,183 @@ pub fn sprt_classify( llr_trace, }) } + +// ================ Confidence-interval (ACI) classification for CAT ========= +// +// `ci_classify` implements the confidence-interval classification stopping +// rule for a single cut score with binary responses: after each response, +// compute the interim EAP ability estimate and its posterior SD, form the +// interval theta_hat +/- z_crit * SE, and classify as soon as the whole +// interval lies STRICTLY on one side of the cut (first strict crossing +// decides). EAP uses the crate CAT convention shared with `eap_interim`: +// a fixed uniform grid of 41 points on [-4, 4], standard-normal log prior +// -0.5 * theta^2 (no quadrature-weight multiplier), and the D = 1 logistic +// 3PL P_i(theta) = c_i + (1 - c_i) / (1 + exp(-a_i (theta - b_i))). +// +// For each prefix k = 1..n: +// theta_hat_k = sum_q w_q theta_q / sum_q w_q +// se_k = sqrt(sum_q w_q (theta_q - theta_hat_k)^2 / sum_q w_q) +// lower_k = theta_hat_k - z_crit * se_k, upper_k = theta_hat_k + z_crit * se_k +// lower_k > theta_cut -> "above" (n_used = k); +// upper_k < theta_cut -> "below" (n_used = k); else continue. +// No crossing -> "continue" with n_used = len(responses). Traces are filled +// for ALL supplied responses; entries after n_used are offline counterfactual +// replay values (live CAT would stop at n_used). +// +// CITATION GOVERNANCE / SCOPE (adversarial spec review, +// ci_classify_spec_review.md): the implemented confidence-interval stopping +// rule was verified against catIrt R/termCI.R, R/eapEst.R, and man/catIrt.Rd +// at commit c9e979e4812c27d95d367a7f097edfe8e93ac8eb (READ): form the +// interval theta_hat +/- z * SEM, where the EAP SEM is the posterior SD +// (sqrt(E[theta^2] - theta_hat^2)), and classify only when the full interval +// lies strictly within a category. The fixed 41-point [-4, 4] EAP grid and +// the caller-supplied z_crit (catIrt computes qnorm((1 + conf.lev) / 2) from +// a confidence level; passing that value is equivalent) are repository +// implementation choices. Thompson (2007), Kingsbury & Weiss (1983), and +// Eggen & Straetmans (2000) were NOT method-section verified in this +// iteration and are historical/background context only. +// +// References (APA 7th): +// Nydick, S. W. (2014). catIrt: An R package for simulating IRT-based +// computerized adaptive tests. (READ: R/termCI.R interval rule and +// strict within-bounds comparisons; R/eapEst.R posterior-SD SEM; +// man/catIrt.Rd conf.lev parameterization and first-satisfied-criterion +// termination) +// Kingsbury, G. G., & Weiss, D. J. (1983). A comparison of IRT-based +// adaptive mastery testing and a sequential mastery testing procedure. +// In D. J. Weiss (Ed.), New horizons in testing (pp. 257-283). Academic +// Press. (NOT read; historical origin of ability-confidence-interval +// classification) +// Thompson, N. A. (2007). A practitioner's guide for variable-length +// computerized classification testing. Practical Assessment, Research & +// Evaluation, 12(1). (NOT read for the CI method section in this +// iteration; background only) +// Eggen, T. J. H. M., & Straetmans, G. J. J. M. (2000). Computerized +// adaptive testing for classifying examinees into three categories. +// Educational and Psychological Measurement, 60(5), 713-734. (NOT read; +// historical) + +/// Result of [`ci_classify`]. `decision` is `"above"`, `"below"`, or +/// `"continue"`; `n_used` is the 1-based count of responses consumed by the +/// first strict interval crossing (or all responses when no crossing +/// occurs); the four traces hold the interim EAP estimate, posterior SD, +/// and interval bounds after every supplied response (entries past `n_used` +/// are offline counterfactuals). +#[derive(Debug, Clone)] +pub struct CiResult { + pub decision: &'static str, + pub n_used: usize, + pub theta_trace: Vec, + pub se_trace: Vec, + pub lower_trace: Vec, + pub upper_trace: Vec, +} + +/// Single-cut binary-response confidence-interval (ACI) classification (see +/// module comment above for the exact verified contract and source status). +pub fn ci_classify( + a: &[f64], + b: &[f64], + c: &[f64], + responses: &[u8], + theta_cut: f64, + z_crit: f64, +) -> Result { + let n = a.len(); + if n == 0 { + return Err("ci_classify: item pool is empty".into()); + } + if b.len() != n || c.len() != n || responses.len() != n { + return Err(format!( + "ci_classify: length mismatch (a: {}, b: {}, c: {}, responses: {})", + n, + b.len(), + c.len(), + responses.len() + )); + } + for i in 0..n { + if !a[i].is_finite() || a[i] <= 0.0 { + return Err(format!("ci_classify: a[{i}] must be finite and > 0")); + } + if !b[i].is_finite() { + return Err(format!("ci_classify: b[{i}] must be finite")); + } + if !c[i].is_finite() || !(0.0..1.0).contains(&c[i]) { + return Err(format!("ci_classify: c[{i}] must be finite and in [0, 1)")); + } + if responses[i] > 1 { + return Err(format!("ci_classify: responses[{i}] must be 0 or 1")); + } + } + if !theta_cut.is_finite() { + return Err("ci_classify: theta_cut must be finite".into()); + } + if !z_crit.is_finite() || z_crit <= 0.0 { + return Err("ci_classify: z_crit must be finite and > 0".into()); + } + + const Q: usize = 41; + let grid: Vec = (0..Q) + .map(|q| -4.0 + 8.0 * q as f64 / (Q - 1) as f64) + .collect(); + // Cumulative log posterior weights, updated one response at a time. + let mut logw: Vec = grid.iter().map(|&t| -0.5 * t * t).collect(); + + let mut theta_trace = Vec::with_capacity(n); + let mut se_trace = Vec::with_capacity(n); + let mut lower_trace = Vec::with_capacity(n); + let mut upper_trace = Vec::with_capacity(n); + let mut decision = "continue"; + let mut n_used = n; + for i in 0..n { + for (q, &t) in grid.iter().enumerate() { + let p = p3pl(t, a[i], b[i], c[i]).clamp(1e-12, 1.0 - 1e-12); + logw[q] += if responses[i] == 1 { + p.ln() + } else { + (1.0 - p).ln() + }; + } + let m = logw.iter().cloned().fold(f64::NEG_INFINITY, f64::max); + let mut den = 0.0; + let mut num = 0.0; + for (q, &t) in grid.iter().enumerate() { + let w = (logw[q] - m).exp(); + den += w; + num += w * t; + } + let theta_hat = num / den; + let mut ss = 0.0; + for (q, &t) in grid.iter().enumerate() { + let w = (logw[q] - m).exp(); + ss += w * (t - theta_hat) * (t - theta_hat); + } + let se = (ss / den).sqrt(); + let lower = theta_hat - z_crit * se; + let upper = theta_hat + z_crit * se; + theta_trace.push(theta_hat); + se_trace.push(se); + lower_trace.push(lower); + upper_trace.push(upper); + // First STRICT crossing decides (catIrt termCI.R uses strict + // within-bounds comparisons; equality means continue). + if decision == "continue" { + if lower > theta_cut { + decision = "above"; + n_used = i + 1; + } else if upper < theta_cut { + decision = "below"; + n_used = i + 1; + } + } + } + Ok(CiResult { + decision, + n_used, + theta_trace, + se_trace, + lower_trace, + upper_trace, + }) +} diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index e42a2e96a..002cc3104 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -98,6 +98,7 @@ ccat_select as ccat_select, epv_select as epv_select, sprt_classify as sprt_classify, + ci_classify as ci_classify, ) from .rsm import fit_rsm as fit_rsm, RsmFit as RsmFit from .mixed import fit_mixed_items as fit_mixed_items, MixedFormatFit as MixedFormatFit, MixedItemParameters as MixedItemParameters @@ -256,6 +257,7 @@ "ccat_select", "epv_select", "sprt_classify", + "ci_classify", "AStratifiedResult", "omega_total_1f", "omega_total_1f_from_data", diff --git a/python/fast_mlsirm/exposure.py b/python/fast_mlsirm/exposure.py index 5bd71ef96..6eb1744fa 100644 --- a/python/fast_mlsirm/exposure.py +++ b/python/fast_mlsirm/exposure.py @@ -713,4 +713,102 @@ def sprt_classify( "n_used": int(r["n_used"]), "llr": float(r["llr"]), "llr_trace": np.asarray(r["llr_trace"]), + } +def ci_classify( + a: np.ndarray, + b: np.ndarray, + c: np.ndarray | None = None, + *, + responses: np.ndarray, + theta_cut: float, + z_crit: float, +) -> dict: + """Single-cut binary-response confidence-interval (ACI) classification. + + After each response, computes the interim EAP ability estimate on a + fixed uniform grid of 41 points on ``[-4, 4]`` with a standard-normal + log prior (``-0.5 * theta**2``, no quadrature-weight multiplier) under + the D = 1 logistic 3PL + ``P_i(theta) = c_i + (1 - c_i) / (1 + exp(-a_i (theta - b_i)))``, plus + the EAP posterior SD as the standard error, and forms the interval + ``theta_hat +/- z_crit * se``. The FIRST STRICT crossing decides: + ``lower > theta_cut`` -> ``"above"``, ``upper < theta_cut`` -> + ``"below"`` (``n_used = k``, 1-based); no crossing -> ``"continue"`` + with ``n_used = len(responses)``. Equality with the cut means continue. + All numerics run in the Rust core + (``mlsirm_core::exposure::ci_classify``). + + Traces are returned for ALL supplied responses as offline diagnostics; + entries past ``n_used`` are counterfactual replay values (a live CAT + would stop at ``n_used``). ``z_crit`` is the normal critical value; for + a confidence level ``L`` pass ``qnorm((1 + L) / 2)`` (catIrt's + ``conf.lev`` parameterization), e.g. 1.6448536269514722 for L = 0.90. + + Source status: the interval stopping rule was verified against R catIrt + ``termCI.R``/``eapEst.R``/``catIrt.Rd`` at commit + c9e979e4812c27d95d367a7f097edfe8e93ac8eb (READ): interval + ``theta_hat +/- z * SEM`` with the EAP SEM equal to the posterior SD, + classifying only when the whole interval lies strictly within a + category. The fixed 41-point grid and caller-supplied ``z_crit`` are + repository implementation choices. Kingsbury & Weiss (1983), Thompson + (2007), and Eggen & Straetmans (2000) were NOT method-section verified + in this iteration and are cited as historical/background context only. + + References (APA 7th ed.): + Kingsbury, G. G., & Weiss, D. J. (1983). A comparison of IRT-based + adaptive mastery testing and a sequential mastery testing + procedure. In D. J. Weiss (Ed.), *New horizons in testing* + (pp. 257-283). Academic Press. (NOT read; historical origin.) + Thompson, N. A. (2007). A practitioner's guide for variable-length + computerized classification testing. *Practical Assessment, + Research & Evaluation, 12*(1). + https://doi.org/10.7275/fq3r-zz60 (NOT read for the CI method + section in this iteration; background only.) + Eggen, T. J. H. M., & Straetmans, G. J. J. M. (2000). Computerized + adaptive testing for classifying examinees into three + categories. *Educational and Psychological Measurement, 60*(5), + 713-734. (NOT read; historical.) + """ + from . import _core + + # Reject complex input BEFORE the dtype casts: the casts would silently + # discard imaginary parts (complex laundering). + for name, arr in (("a", a), ("b", b), ("c", c), ("responses", responses)): + if arr is not None and np.iscomplexobj(np.asarray(arr)): + raise ValueError(f"{name} must be real-valued") + a = np.asarray(a, dtype=np.float64) + b = np.asarray(b, dtype=np.float64) + if a.ndim != 1 or b.ndim != 1: + raise ValueError("a and b must be 1-D arrays") + if c is None: + c = np.zeros_like(a) + c = np.asarray(c, dtype=np.float64) + if c.ndim != 1: + raise ValueError("c must be a 1-D array") + # Validate responses BEFORE the uint8 cast (casts truncate/wrap). + resp = np.asarray(responses) + if resp.ndim != 1: + raise ValueError("responses must be a 1-D array") + if resp.dtype == np.bool_: + resp = resp.astype(np.uint8) + else: + resp_f = np.asarray(resp, dtype=np.float64) + if not np.all(np.isin(resp_f, (0.0, 1.0))): + raise ValueError("responses must contain only 0 and 1") + resp = resp_f.astype(np.uint8) + r = _core.py_ci_classify( + np.ascontiguousarray(a), + np.ascontiguousarray(b), + np.ascontiguousarray(c), + np.ascontiguousarray(resp), + float(theta_cut), + float(z_crit), + ) + return { + "decision": str(r["decision"]), + "n_used": int(r["n_used"]), + "theta_trace": np.asarray(r["theta_trace"]), + "se_trace": np.asarray(r["se_trace"]), + "lower_trace": np.asarray(r["lower_trace"]), + "upper_trace": np.asarray(r["upper_trace"]), } \ No newline at end of file diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index 904c70aa5..a414baa2a 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -6716,4 +6716,112 @@ def test_core_validation_propagates(self): delta=0.5, alpha=0.6, beta=0.5, + ) + +class TestCiClassify: + """ci_classify wrapper: every assert reads crate-returned dict values.""" + + def test_pinned_oracle(self): + from fast_mlsirm import ci_classify + + r = ci_classify( + np.full(6, 1.5), + np.array([-1.5, -0.9, -0.3, 0.3, 0.9, 1.5]), + np.zeros(6), + responses=np.array([1, 1, 1, 1, 1, 0]), + theta_cut=0.0, + z_crit=1.6448536269514722, + ) + assert r["decision"] == "above" + assert r["n_used"] == 5 + theta = np.array( + [ + 0.18783548849905624, + 0.40433637208107137, + 0.65453031321107147, + 0.93795666218057705, + 1.251068565832161, + 1.004851105902542, + ] + ) + se = np.array( + [ + 0.91459937771477151, + 0.84249780260178286, + 0.78082991898905685, + 0.72897935456728113, + 0.68628322205747161, + 0.60091214158918393, + ] + ) + np.testing.assert_allclose(r["theta_trace"], theta, rtol=0, atol=1e-12) + np.testing.assert_allclose(r["se_trace"], se, rtol=0, atol=1e-12) + np.testing.assert_allclose( + r["lower_trace"], theta - 1.6448536269514722 * se, rtol=0, atol=1e-12 + ) + np.testing.assert_allclose( + r["upper_trace"], theta + 1.6448536269514722 * se, rtol=0, atol=1e-12 + ) + # First-strict-crossing anchor: k=5 is the first lower bound > cut. + assert r["lower_trace"][3] <= 0.0 < r["lower_trace"][4] + + def test_default_c_and_bool_responses(self): + from fast_mlsirm import ci_classify + + r = ci_classify( + np.array([1.0, 1.1]), + np.array([0.0, 0.2]), + responses=np.array([True, False]), + theta_cut=0.0, + z_crit=1.96, + ) + assert r["decision"] == "continue" + assert r["n_used"] == 2 + assert r["theta_trace"].shape == (2,) + assert np.all(r["se_trace"] > 0.0) + assert np.all(r["lower_trace"] <= 0.0) + assert np.all(r["upper_trace"] >= 0.0) + + def test_rejects_bad_responses(self): + from fast_mlsirm import ci_classify + + with pytest.raises(ValueError, match="0 and 1"): + ci_classify( + np.array([1.0, 1.0]), + np.zeros(2), + responses=np.array([1, 2]), + theta_cut=0.0, + z_crit=1.96, + ) + with pytest.raises(ValueError, match="0 and 1"): + ci_classify( + np.array([1.0, 1.0]), + np.zeros(2), + responses=np.array([1.0, 0.5]), + theta_cut=0.0, + z_crit=1.96, + ) + + def test_rejects_complex_input(self): + from fast_mlsirm import ci_classify + + with pytest.raises(ValueError, match="real-valued"): + ci_classify( + np.array([1.0 + 1j, 1.0]), + np.zeros(2), + responses=np.array([1, 0]), + theta_cut=0.0, + z_crit=1.96, + ) + + def test_core_validation_propagates(self): + from fast_mlsirm import ci_classify + + with pytest.raises(ValueError, match="z_crit"): + ci_classify( + np.array([1.0, 1.0]), + np.zeros(2), + responses=np.array([1, 0]), + theta_cut=0.0, + z_crit=0.0, ) \ No newline at end of file diff --git a/tests/unit/exposure_tests.rs b/tests/unit/exposure_tests.rs index a6d86600e..477bf23dc 100644 --- a/tests/unit/exposure_tests.rs +++ b/tests/unit/exposure_tests.rs @@ -28,8 +28,9 @@ //! `sh_controls_max_exposure`. use crate::exposure::{ - a_stratified, ccat_select, eap_interim, epv_select, kl_information, kl_select, owen_cat, - owen_update, p3pl, sprt_classify, sympson_hetter, AStratifiedConfig, Lcg, SympsonHetterConfig, + a_stratified, ccat_select, ci_classify, eap_interim, epv_select, kl_information, kl_select, + owen_cat, owen_update, p3pl, sprt_classify, sympson_hetter, AStratifiedConfig, Lcg, + SympsonHetterConfig, }; fn pool30() -> (Vec, Vec, Vec) { @@ -1861,3 +1862,216 @@ fn sprt_extreme_parameters_stay_finite() { assert!((rc.llr - (-0.2_f64.ln())).abs() < 1e-12, "llr = {}", rc.llr); assert!(rc.llr.is_finite()); } + +// ---------- ci_classify (confidence-interval / ACI classification) ---------- + +/// Pinned 17-digit oracle from the adversarial spec review +/// (ci_classify_spec_review.md): independent Python recomputation of the +/// approved 41-point [-4,4] EAP posterior-SD rule. Every assert reads crate +/// outputs (decision, n_used, the four traces). Kills mutants: M1 swapped +/// decisions (expects "above"), M2 point-estimate-vs-cut (theta_trace[0] > 0 +/// would decide at k=1, oracle n_used = 5), M3 variance-instead-of-SD (that +/// mutant crosses at k=4; the lower_trace[3] <= 0 < lower_trace[4] anchor +/// pins the first strict crossing to k=5), M4 n_used off-by-one, M5 +/// final-CI-only (the counterfactual tail is also above the cut, so n_used +/// and the crossing-index anchors are the discriminating asserts, not the +/// final decision alone). +#[test] +fn ci_classify_pinned_oracle() { + let a = [1.5; 6]; + let b = [-1.5, -0.9, -0.3, 0.3, 0.9, 1.5]; + let c = [0.0; 6]; + let responses = [1u8, 1, 1, 1, 1, 0]; + let r = ci_classify(&a, &b, &c, &responses, 0.0, 1.6448536269514722).unwrap(); + assert_eq!(r.decision, "above"); + assert_eq!(r.n_used, 5); + let theta_exp = [ + 0.18783548849905624, + 0.40433637208107137, + 0.65453031321107147, + 0.93795666218057705, + 1.251068565832161, + 1.004851105902542, + ]; + let se_exp = [ + 0.91459937771477151, + 0.84249780260178286, + 0.78082991898905685, + 0.72897935456728113, + 0.68628322205747161, + 0.60091214158918393, + ]; + let lower_exp = [ + -1.316546615142645, + -0.98144919422711663, + -0.62982061107030296, + -0.26110767315215866, + 0.12223311891498612, + 0.016438590330396186, + ]; + let upper_exp = [ + 1.6922175921407576, + 1.7901219383892593, + 1.938881237492446, + 2.1370209975133125, + 2.3799040127493356, + 1.9932636214746879, + ]; + assert_eq!(r.theta_trace.len(), 6); + assert_eq!(r.se_trace.len(), 6); + assert_eq!(r.lower_trace.len(), 6); + assert_eq!(r.upper_trace.len(), 6); + for k in 0..6 { + assert!( + (r.theta_trace[k] - theta_exp[k]).abs() < 1e-12, + "theta[{k}] = {}", + r.theta_trace[k] + ); + assert!( + (r.se_trace[k] - se_exp[k]).abs() < 1e-12, + "se[{k}] = {}", + r.se_trace[k] + ); + assert!( + (r.lower_trace[k] - lower_exp[k]).abs() < 1e-12, + "lower[{k}] = {}", + r.lower_trace[k] + ); + assert!( + (r.upper_trace[k] - upper_exp[k]).abs() < 1e-12, + "upper[{k}] = {}", + r.upper_trace[k] + ); + } + // First-strict-crossing anchor (kills M3/M5): no crossing before k=5. + assert!(r.lower_trace[3] <= 0.0 && r.lower_trace[4] > 0.0); + for k in 0..4 { + assert!(r.lower_trace[k] <= 0.0 && r.upper_trace[k] >= 0.0); + } +} + +/// "below" decision on all-wrong responses with a positive cut, and a +/// "continue" outcome when z_crit is too wide to ever cross. Asserts read +/// crate decision/n_used/bound traces. +#[test] +fn ci_classify_below_and_continue() { + let a = [1.5; 6]; + let b = [-1.5, -0.9, -0.3, 0.3, 0.9, 1.5]; + let c = [0.0; 6]; + let wrong = [0u8; 6]; + let r = ci_classify(&a, &b, &c, &wrong, 0.5, 1.6448536269514722).unwrap(); + assert_eq!(r.decision, "below"); + assert!(r.n_used <= 6); + let k = r.n_used - 1; + assert!(r.upper_trace[k] < 0.5, "upper = {}", r.upper_trace[k]); + for j in 0..k { + assert!(r.upper_trace[j] >= 0.5 || r.lower_trace[j] > 0.5); + } + // Huge z_crit: interval always straddles any interior cut -> continue. + let rc = ci_classify(&a, &b, &c, &wrong, 0.5, 100.0).unwrap(); + assert_eq!(rc.decision, "continue"); + assert_eq!(rc.n_used, 6); + for j in 0..6 { + assert!(rc.lower_trace[j] <= 0.5 && rc.upper_trace[j] >= 0.5); + } +} + +/// Full validation error paths; each assert reads the crate Err string. +#[test] +fn ci_classify_error_paths() { + let ok_a = [1.0]; + let ok_b = [0.0]; + let ok_c = [0.0]; + let ok_r = [1u8]; + assert!(ci_classify(&[], &[], &[], &[], 0.0, 1.96) + .unwrap_err() + .contains("empty")); + assert!(ci_classify(&ok_a, &[0.0, 1.0], &ok_c, &ok_r, 0.0, 1.96) + .unwrap_err() + .contains("length mismatch")); + assert!(ci_classify(&[-1.0], &ok_b, &ok_c, &ok_r, 0.0, 1.96) + .unwrap_err() + .contains("a[0]")); + assert!(ci_classify(&[f64::NAN], &ok_b, &ok_c, &ok_r, 0.0, 1.96) + .unwrap_err() + .contains("a[0]")); + assert!( + ci_classify(&ok_a, &[f64::INFINITY], &ok_c, &ok_r, 0.0, 1.96) + .unwrap_err() + .contains("b[0]") + ); + assert!(ci_classify(&ok_a, &ok_b, &[1.0], &ok_r, 0.0, 1.96) + .unwrap_err() + .contains("c[0]")); + assert!(ci_classify(&ok_a, &ok_b, &[-0.1], &ok_r, 0.0, 1.96) + .unwrap_err() + .contains("c[0]")); + assert!(ci_classify(&ok_a, &ok_b, &ok_c, &[2], 0.0, 1.96) + .unwrap_err() + .contains("responses[0]")); + assert!(ci_classify(&ok_a, &ok_b, &ok_c, &ok_r, f64::NAN, 1.96) + .unwrap_err() + .contains("theta_cut")); + assert!(ci_classify(&ok_a, &ok_b, &ok_c, &ok_r, 0.0, 0.0) + .unwrap_err() + .contains("z_crit")); + assert!(ci_classify(&ok_a, &ok_b, &ok_c, &ok_r, 0.0, f64::NAN) + .unwrap_err() + .contains("z_crit")); +} + +/// MC-500 structural invariants on random pools/responses. All asserts read +/// crate outputs: trace lengths, SE positivity/monotonic bounds, decision +/// consistency with the returned interval at n_used, and no crossing before +/// n_used. +#[test] +#[ignore = "500-rep Monte Carlo; run explicitly"] +fn ci_classify_mc500_invariants() { + let mut rng = Lcg(20260220); + for rep in 0..500 { + let n = 3 + (rng.next_f64() * 18.0) as usize; + let mut a = Vec::with_capacity(n); + let mut b = Vec::with_capacity(n); + let mut c = Vec::with_capacity(n); + let mut resp = Vec::with_capacity(n); + for _ in 0..n { + a.push(0.5 + 2.0 * rng.next_f64()); + b.push(-2.5 + 5.0 * rng.next_f64()); + c.push(0.25 * rng.next_f64()); + resp.push(if rng.next_f64() < 0.5 { 1u8 } else { 0u8 }); + } + let cut = -1.5 + 3.0 * rng.next_f64(); + let z = 0.5 + 2.0 * rng.next_f64(); + let r = ci_classify(&a, &b, &c, &resp, cut, z).unwrap(); + assert_eq!(r.theta_trace.len(), n, "rep {rep}"); + assert_eq!(r.se_trace.len(), n, "rep {rep}"); + assert!(r.n_used >= 1 && r.n_used <= n, "rep {rep}"); + for k in 0..n { + assert!( + r.se_trace[k].is_finite() && r.se_trace[k] > 0.0, + "rep {rep}" + ); + assert!( + r.theta_trace[k] > -4.0 && r.theta_trace[k] < 4.0, + "rep {rep}" + ); + let lo = r.theta_trace[k] - z * r.se_trace[k]; + let hi = r.theta_trace[k] + z * r.se_trace[k]; + assert!((r.lower_trace[k] - lo).abs() < 1e-12, "rep {rep}"); + assert!((r.upper_trace[k] - hi).abs() < 1e-12, "rep {rep}"); + } + let k = r.n_used - 1; + match r.decision { + "above" => assert!(r.lower_trace[k] > cut, "rep {rep}"), + "below" => assert!(r.upper_trace[k] < cut, "rep {rep}"), + "continue" => assert_eq!(r.n_used, n, "rep {rep}"), + other => panic!("rep {rep}: unexpected decision {other}"), + } + for j in 0..k { + assert!( + r.lower_trace[j] <= cut && r.upper_trace[j] >= cut, + "rep {rep}: crossing before n_used at {j}" + ); + } + } +} From d0a3573dc8b90524bb97ef7e0c440becc1b42ba0 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 26 Jul 2026 00:27:31 +0900 Subject: [PATCH 2/2] Complete NOT-read citation list in py_ci_classify docstring Impl-review doc-defect: the PyO3 docstring omitted Eggen & Straetmans (2000) from the NOT-read historical citations required by the approved citation-governance wording (Rust core, Python wrapper, and CHANGELOG already carried it). Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- crates/fast-mlsirm-py/src/lib.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index e88890802..d5cab3350 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -2498,6 +2498,8 @@ fn py_sprt_classify( /// historical origin) /// Thompson, N. A. (2007). Practical Assessment, Research & Evaluation, /// 12(1). (NOT read for the CI method section; background only) +/// Eggen, T. J. H. M., & Straetmans, G. J. J. M. (2000). Educational and +/// Psychological Measurement, 60(5), 713-734. (NOT read; historical) #[pyfunction] fn py_ci_classify( py: Python<'_>,