feat: Wald SPRT classification for CAT (single-cut binary, reduced scope)

CHANGELOG.md

@@ -119,6 +119,23 @@
 
 ### Added
 
+- **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
+  log-likelihood ratio under the D=1 logistic 3PL, and inclusive
+  first-crossing decisions against the log Wald boundaries
+  `A = ln((1-beta)/alpha)`, `B = ln(beta/(1-alpha))` -> `"above"`/`"below"`/
+  `"continue"` with 1-based `n_used`; the full `llr_trace` is returned as an
+  offline diagnostic (entries past `n_used` are counterfactual replay
+  values). Verified against R catIrt `termSPRT.R`/`logLik.brm.R`/`p.brm.R`
+  and Thompson (2007, doi:10.7275/fq3r-zz60); Reckase (1983), Eggen (1999),
+  and Wald (1947) are cited as historical origins via Thompson (not directly
+  read). Log-likelihood ratios are computed in stable log space (softplus /
+  log-sigmoid), so extreme-but-valid parameters that saturate the response
+  probability to numerical 0/1 yield finite LLRs instead of errors. Pinned
+  17-digit interior-crossing oracle, error-path and 500-rep Monte-Carlo
+  structural-invariant tests; 4 executed mutation kills (swapped boundaries,
+  dropped guessing floor, collapsed null hypothesis, off-by-one `n_used`).
 - **Owen-approximate posterior-predictive EPV item selection**
   (`fast_mlsirm.epv_select`; in `mlsirm_core::exposure`). Deliberately
   reduced scope of van der Linden's (1998, doi:10.1007/BF02294775) minimum

crates/fast-mlsirm-py/src/lib.rs

@@ -42,7 +42,8 @@ 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,
-    sympson_hetter as core_sympson_hetter, AStratifiedConfig, SympsonHetterConfig,
+    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::{
@@ -2430,6 +2431,57 @@ fn py_epv_select(
     Ok(out.into())
 }
 
+/// Single-cut binary-response Wald SPRT classification for CAT
+/// (`mlsirm_core::exposure::sprt_classify`). D = 1 logistic 3PL; point
+/// hypotheses at `theta_cut -/+ delta`; log Wald boundaries
+/// A = ln((1-beta)/alpha), B = ln(beta/(1-alpha)) with inclusive
+/// first-crossing decisions ("above"/"below"/"continue"). `llr_trace`
+/// entries past `n_used` are offline counterfactual replay values.
+///
+/// References (APA 7th; see the core module comment for read/not-read
+/// source status):
+/// 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 (READ)
+/// Nydick, S. W. (2014). catIrt (R package). (READ: termSPRT.R,
+/// logLik.brm.R, p.brm.R)
+/// Eggen, T. J. H. M. (1999). Applied Psychological Measurement, 23(3),
+/// 249-261. (NOT read; historical citation via Thompson)
+/// Reckase, M. D. (1983). A procedure for decision making using tailored
+/// testing. (NOT read; historical citation via Thompson)
+/// Wald, A. (1947). Sequential analysis. Wiley. (NOT read; boundary forms
+/// verified through the READ sources above)
+#[pyfunction]
+fn py_sprt_classify(
+    py: Python<'_>,
+    a: PyReadonlyArray1<'_, f64>,
+    b: PyReadonlyArray1<'_, f64>,
+    c: PyReadonlyArray1<'_, f64>,
+    responses: PyReadonlyArray1<'_, u8>,
+    theta_cut: f64,
+    delta: f64,
+    alpha: f64,
+    beta: f64,
+) -> PyResult<Py<pyo3::types::PyDict>> {
+    let res = core_sprt_classify(
+        a.as_slice()?,
+        b.as_slice()?,
+        c.as_slice()?,
+        responses.as_slice()?,
+        theta_cut,
+        delta,
+        alpha,
+        beta,
+    )
+    .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("llr", res.llr)?;
+    out.set_item("llr_trace", numpy::PyArray1::from_slice(py, &res.llr_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
@@ -6250,6 +6302,7 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(py_owen_cat, m)?)?;
     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!(guttman_lambdas, m)?)?;
     m.add_function(wrap_pyfunction!(tenberge_mu, m)?)?;
     m.add_function(wrap_pyfunction!(cronbach_alpha, m)?)?;

crates/mlsirm-core/src/exposure.rs

@@ -1403,3 +1403,192 @@ pub fn epv_select(
         predictive,
     })
 }
+
+// ===================== Wald SPRT classification for CAT =====================
+//
+// `sprt_classify` implements single-cut, binary-response SPRT classification
+// (Wald's sequential probability ratio test applied to IRT classification
+// testing). Two point hypotheses around the cut score,
+//   theta0 = theta_cut - delta,  theta1 = theta_cut + delta,
+// are compared through the cumulative binary log-likelihood ratio under the
+// D = 1 logistic 3PL
+//   P_i(theta) = c_i + (1 - c_i) / (1 + exp(-a_i (theta - b_i))),
+//   LLR_k = sum_{i<=k} [ u_i ln(P_i(theta1)/P_i(theta0))
+//                      + (1 - u_i) ln((1 - P_i(theta1))/(1 - P_i(theta0))) ],
+// against the log Wald boundaries
+//   A = ln((1 - beta) / alpha),  B = ln(beta / (1 - alpha)).
+// Responses are walked in order and the FIRST crossing decides (inclusive
+// comparisons, matching catIrt): LLR_k >= A -> "above" with n_used = k;
+// LLR_k <= B -> "below" with n_used = k; no crossing -> "continue" with
+// n_used = len(responses).
+//
+// CITATION GOVERNANCE / SCOPE (adversarial spec review, sprt_spec_review.md):
+// boundaries and the binary log-likelihood-ratio form were verified against
+// READ sources: catIrt R/termSPRT.R + R/logLik.brm.R + R/p.brm.R (GitHub
+// swnydick/catIrt) and Thompson (2007), p. 7. Reckase (1983) and Eggen
+// (1999) are historical citations via Thompson and were NOT directly read.
+// This function implements only a single-cut binary 3PL SPRT with D = 1
+// logistic-scale item parameters; it is not a multi-cut, polytomous, or
+// D = 1.7 compatibility layer (parameters calibrated on the D = 1.7 metric
+// must be rescaled by the caller, a_D1 = 1.7 * a_D17, before use).
+//
+// The returned decision/n_used are first-crossing SPRT results. llr_trace is
+// computed for ALL supplied responses as an offline diagnostic; entries after
+// n_used are counterfactual replay values - live CAT would terminate at
+// n_used and would not administer later items.
+//
+// References (APA 7th):
+// Wald, A. (1947). Sequential analysis. Wiley. (NOT read; boundary forms
+//     verified through the sources below)
+// 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 (READ: p. 7
+//     likelihood-ratio form and Wald decision points)
+// Nydick, S. W. (2014). catIrt: An R package for simulating IRT-based
+//     computerized adaptive tests. (READ: R/termSPRT.R boundary and
+//     inclusive-comparison conventions; R/logLik.brm.R binary log
+//     likelihood; R/p.brm.R D = 1 3PL)
+// Eggen, T. J. H. M. (1999). Item selection in adaptive testing with the
+//     sequential probability ratio test. Applied Psychological Measurement,
+//     23(3), 249-261. (NOT read; historical citation via Thompson)
+// Reckase, M. D. (1983). A procedure for decision making using tailored
+//     testing. (NOT read; historical citation via Thompson)
+
+/// Result of [`sprt_classify`]. `decision` is `"above"`, `"below"`, or
+/// `"continue"`; `n_used` is the 1-based count of responses consumed by the
+/// first boundary crossing (or all responses when no crossing occurs);
+/// `llr_trace` holds the cumulative log-likelihood ratio after every supplied
+/// response (entries past `n_used` are offline counterfactuals); `llr` is the
+/// final trace entry.
+#[derive(Debug, Clone)]
+pub struct SprtResult {
+    pub decision: &'static str,
+    pub n_used: usize,
+    pub llr_trace: Vec<f64>,
+    pub llr: f64,
+}
+
+/// Single-cut binary-response Wald SPRT classification (see module comment
+/// above for the exact verified contract and source status).
+pub fn sprt_classify(
+    a: &[f64],
+    b: &[f64],
+    c: &[f64],
+    responses: &[u8],
+    theta_cut: f64,
+    delta: f64,
+    alpha: f64,
+    beta: f64,
+) -> Result<SprtResult, String> {
+    let n = a.len();
+    if n == 0 {
+        return Err("sprt_classify: item pool is empty".into());
+    }
+    if b.len() != n || c.len() != n || responses.len() != n {
+        return Err(format!(
+            "sprt_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!("sprt_classify: a[{i}] must be finite and > 0"));
+        }
+        if !b[i].is_finite() {
+            return Err(format!("sprt_classify: b[{i}] must be finite"));
+        }
+        if !c[i].is_finite() || !(0.0..1.0).contains(&c[i]) {
+            return Err(format!(
+                "sprt_classify: c[{i}] must be finite and in [0, 1)"
+            ));
+        }
+        if responses[i] > 1 {
+            return Err(format!("sprt_classify: responses[{i}] must be 0 or 1"));
+        }
+    }
+    if !theta_cut.is_finite() {
+        return Err("sprt_classify: theta_cut must be finite".into());
+    }
+    if !delta.is_finite() || delta <= 0.0 {
+        return Err("sprt_classify: delta must be finite and > 0".into());
+    }
+    for (name, v) in [("alpha", alpha), ("beta", beta)] {
+        if !v.is_finite() || v <= 0.0 || v >= 1.0 {
+            return Err(format!(
+                "sprt_classify: {name} must be finite and in (0, 1)"
+            ));
+        }
+    }
+    if alpha + beta >= 1.0 {
+        return Err("sprt_classify: alpha + beta must be < 1".into());
+    }
+
+    let upper = ((1.0 - beta) / alpha).ln();
+    let lower = (beta / (1.0 - alpha)).ln();
+    let theta0 = theta_cut - delta;
+    let theta1 = theta_cut + delta;
+    // Stable softplus ln(1 + e^z): shift by max(z, 0) so exp never overflows.
+    let softplus = |z: f64| -> f64 {
+        if z > 0.0 {
+            z + (-z).exp().ln_1p()
+        } else {
+            z.exp().ln_1p()
+        }
+    };
+    // Stable log-probabilities under the D = 1 logistic 3PL
+    // P = c + (1 - c) sigmoid(z), z = a (theta - b) (crate CAT convention;
+    // catIrt p.brm.R). ln(1 - P) = ln(1 - c) - softplus(z) always; ln(P)
+    // needs the log-sigmoid branch -softplus(-z) only when c = 0 (for c > 0
+    // the direct form is bounded below by c and stays finite).
+    let ln_p = |z: f64, ci: f64| -> f64 {
+        if ci > 0.0 {
+            (ci + (1.0 - ci) / (1.0 + (-z).exp())).ln()
+        } else {
+            -softplus(-z)
+        }
+    };
+
+    let mut llr_trace = Vec::with_capacity(n);
+    let mut cum = 0.0_f64;
+    let mut decision = "continue";
+    let mut n_used = n;
+    for i in 0..n {
+        let z0 = a[i] * (theta0 - b[i]);
+        let z1 = a[i] * (theta1 - b[i]);
+        let inc = if responses[i] == 1 {
+            // ln(P(theta1)) - ln(P(theta0)), each log computed stably.
+            ln_p(z1, c[i]) - ln_p(z0, c[i])
+        } else {
+            // ln(1-P(theta1)) - ln(1-P(theta0)); the ln(1-c) terms cancel.
+            softplus(z0) - softplus(z1)
+        };
+        // Defensive: unreachable for validated inputs with the stable forms
+        // above (kept as a hard failure rather than silently propagating).
+        if !inc.is_finite() {
+            return Err(format!(
+                "sprt_classify: non-finite log-likelihood-ratio increment at item {i}"
+            ));
+        }
+        cum += inc;
+        llr_trace.push(cum);
+        // First crossing decides; inclusive comparisons (catIrt termSPRT.R).
+        if decision == "continue" {
+            if cum >= upper {
+                decision = "above";
+                n_used = i + 1;
+            } else if cum <= lower {
+                decision = "below";
+                n_used = i + 1;
+            }
+        }
+    }
+    Ok(SprtResult {
+        decision,
+        n_used,
+        llr: *llr_trace.last().unwrap(),
+        llr_trace,
+    })
+}

python/fast_mlsirm/__init__.py

@@ -97,6 +97,7 @@
     owen_cat as owen_cat,
     ccat_select as ccat_select,
     epv_select as epv_select,
+    sprt_classify as sprt_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
@@ -254,6 +255,7 @@
     "owen_cat",
     "ccat_select",
     "epv_select",
+    "sprt_classify",
     "AStratifiedResult",
     "omega_total_1f",
     "omega_total_1f_from_data",