feat(ksirt): kernel-smoothing nonparametric IRT (Ramsay 1991)

CHANGELOG.md

@@ -114,6 +114,52 @@
 
 ### Added
 
+- **Haberman subscore added-value analysis** (`fast_mlsirm.subscore_analysis`;
+  new `mlsirm_core::subscores`; Haberman, 2008, as cited in Sinharay, 2010).
+  For each subscale of a disjoint, exhaustive item partition computes the
+  PRMSEs of the three classical-test-theory true-subscore estimators — from
+  the observed subscore (`= Cronbach alpha`), from the observed total
+  (`rho^2(s_t, x_t) * alpha_x` with the true-score covariance row sum over
+  subscore columns only), and from both jointly (Wainer-style augmentation via
+  `tau`/`beta`/`gamma`) — plus per-person estimator matrices, the
+  `(K+1)^2` score correlation matrix, disattenuated subscore correlations, and
+  added-value decisions (Haberman's `PRMSE_s > PRMSE_x`; Sinharay's 2010
+  `+ 0.01` margin for augmentation, labeled — CRAN `CTTsub`'s relative rule is
+  documented but not implemented). Formulas verified against the Appendix of
+  Sinharay (2010, ETS RR-10-16) and the CRAN `subscore` R source read
+  line-by-line; degenerate samples (alpha outside `(0, 1]`, zero variance,
+  subscore collinear with the total) are rejected instead of propagating NaN.
+  For LLM-as-a-Judge item-quality management this decides whether per-domain
+  judge subscores add diagnostic value over the overall score. Rust-only
+  numerics; the Python wrapper validates and marshals. Tests pin every
+  reported statistic against literals from an independent NumPy transcription
+  of the R semantics on an asymmetric fixture with mixed added-value
+  outcomes, include rejection tests for the structural and degeneracy guards
+  (the defensive computed-PRMSE-range guard is not separately exercised), a
+  conditional dominance
+  sweep on guard-passing random data, and a 500-rep `#[ignore]` Monte Carlo
+  MSE comparison; three mutation spot-checks (dropped `m/(m-1)`, rowsum
+  including the total column, `tau` numerator sign flip) were run and killed.
+- **Kernel-smoothing nonparametric IRT** (`fast_mlsirm.ksirt_analysis`; new
+  `mlsirm_core::ksirt`; Ramsay, 1991, as cited in Mazza et al., 2014).
+  Estimates option characteristic curves by Nadaraya-Watson kernel regression
+  (gaussian/quadratic/uniform kernels) of option indicators on rank-based
+  ordinal ability estimates `qnorm(rank/(n+1))`, on an equally spaced
+  evaluation grid, with Silverman-rule default bandwidths, plus expected item
+  score and expected total score curves. Formulas verified against the
+  KernSmoothIRT JSS paper (Mazza et al., 2014, Sections 2-2.3) and the
+  KernSmoothIRT R/C++ package source read line-by-line; standard errors and
+  cross-validation bandwidth selection are deliberately out of scope (the R
+  implementation's SE accumulator is order-dependent and unverifiable from
+  read sources). For LLM-as-a-Judge item-quality management this reveals
+  non-monotone or poorly discriminating evaluation items without a parametric
+  model. Rust-only numerics; the Python wrapper validates and marshals. Tests
+  pin a hand-computed 4-person fixture (rank->theta qnorm literals, grid
+  endpoints, Silverman constant), enforce structural invariants
+  (row-sums-to-one with positive denominators, compact-support zeros,
+  zero-denominator fallback), and include a 500-replication Monte Carlo
+  recovery study (`#[ignore]`) under normal and skewed ability generation
+  using the rank-invariance composition oracle.
 - **Mokken scale analysis** (`fast_mlsirm.mokken_analysis`; new
   `mlsirm_core::mokken`; Mokken, 1971, as cited in van der Ark, 2007).
   Computes the Loevinger scalability coefficients `Hij`, `Hi`, `H` and their

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

@@ -60,6 +60,8 @@ use mlsirm_core::rasch_cml::{
     andersen_lr_test as core_andersen_lr, fit_rasch_cml as core_fit_rasch_cml,
 };
 use mlsirm_core::facets::fit_facets as core_fit_facets;
+use mlsirm_core::ksirt::{ksirt as core_ksirt, KsirtKernel};
+use mlsirm_core::subscores::subscores as core_subscores;
 use mlsirm_core::mokken::{aisp as core_mokken_aisp, coef_h as core_mokken_coef_h};
 use mlsirm_core::rsm::fit_rsm as core_fit_rsm;
 use mlsirm_core::rt::{
@@ -1513,6 +1515,124 @@ fn mokken_aisp(
     core_mokken_aisp(x.as_slice()?, n_persons, n_items, c, alpha).map_err(PyValueError::new_err)
 }
 
+/// Kernel-smoothing nonparametric option characteristic curves
+/// (`mlsirm_core::ksirt`; Ramsay, 1991, as cited in Mazza, Punzo, &
+/// McGuire, 2014, https://doi.org/10.18637/jss.v058.i06). `x` is a
+/// row-major complete `n_persons * n_items` pre-scored response matrix.
+/// Returns a dict with `theta` (`N`), `grid` (`Q`), `bandwidth` (`J`), and
+/// per-item lists `options`, `occ` (flattened `m_j * Q`, row-major by
+/// option), `expected` (`Q`), plus `expected_total` (`Q`).
+#[pyfunction]
+#[pyo3(signature = (x, n_persons, n_items, kernel = "gaussian", nevalpoints = 51, bandwidth = None))]
+fn ksirt_occ(
+    py: Python<'_>,
+    x: PyReadonlyArray1<'_, f64>,
+    n_persons: usize,
+    n_items: usize,
+    kernel: &str,
+    nevalpoints: usize,
+    bandwidth: Option<Vec<f64>>,
+) -> PyResult<Py<pyo3::types::PyDict>> {
+    let flat = x.as_slice()?;
+    if flat.len() != n_persons * n_items {
+        return Err(PyValueError::new_err(format!(
+            "x has {} entries, expected n_persons * n_items = {}",
+            flat.len(),
+            n_persons * n_items
+        )));
+    }
+    let kern = match kernel {
+        "gaussian" => KsirtKernel::Gaussian,
+        "quadratic" => KsirtKernel::Quadratic,
+        "uniform" => KsirtKernel::Uniform,
+        other => {
+            return Err(PyValueError::new_err(format!(
+                "unknown kernel '{other}' (expected gaussian, quadratic, or uniform)"
+            )))
+        }
+    };
+    let rows: Vec<&[f64]> = flat.chunks_exact(n_items).collect();
+    let res = core_ksirt(&rows, kern, nevalpoints, bandwidth.as_deref())
+        .map_err(PyValueError::new_err)?;
+    let out = pyo3::types::PyDict::new(py);
+    out.set_item("theta", res.theta)?;
+    out.set_item("grid", res.grid)?;
+    out.set_item("bandwidth", res.bandwidth)?;
+    out.set_item("expected_total", res.expected_total)?;
+    let options: Vec<Vec<f64>> = res.items.iter().map(|it| it.options.clone()).collect();
+    let occ: Vec<Vec<f64>> = res
+        .items
+        .iter()
+        .map(|it| it.occ.iter().flatten().copied().collect())
+        .collect();
+    let expected: Vec<Vec<f64>> = res.items.iter().map(|it| it.expected.clone()).collect();
+    out.set_item("options", options)?;
+    out.set_item("occ", occ)?;
+    out.set_item("expected", expected)?;
+    Ok(out.into())
+}
+
+/// Haberman subscore added-value analysis (`mlsirm_core::subscores`;
+/// Haberman, 2008, as cited in Sinharay, 2010,
+/// ETS RR-10-16). `x` is a row-major complete `n_persons * n_items`
+/// scored response matrix; `groups[j]` in `0..K` assigns item `j` to a
+/// subscale. Returns a dict with per-subscale `alpha`, `prmse_s`,
+/// `prmse_x`, `prmse_sx`, `tau`, `beta`, `gamma`, `added_value_s`,
+/// `added_value_sx`, `alpha_total`, the `(K+1)^2` flattened `corr`, the
+/// `K*K` flattened `disattenuated_corr` (NaN diagonal), and the `n*K`
+/// flattened estimator matrices `observed`, `subscore_s`, `subscore_x`,
+/// `subscore_sx` plus `total` (`n`).
+#[pyfunction]
+fn subscore_analysis(
+    py: Python<'_>,
+    x: PyReadonlyArray1<'_, f64>,
+    n_persons: usize,
+    n_items: usize,
+    groups: Vec<usize>,
+) -> PyResult<Py<pyo3::types::PyDict>> {
+    let flat = x.as_slice()?;
+    // Validate BEFORE allocating rows: unchecked n_persons * n_items can
+    // wrap on 64-bit (e.g. 2^63 * 2 == 0, matching an empty array) and then
+    // panic with capacity overflow inside the row allocation.
+    if n_persons < 3 || n_items < 4 || groups.len() != n_items {
+        return Err(PyValueError::new_err(
+            "need n_persons >= 3, n_items >= 4, and one group index per item",
+        ));
+    }
+    let expected = n_persons
+        .checked_mul(n_items)
+        .ok_or_else(|| PyValueError::new_err("n_persons * n_items overflows"))?;
+    if flat.len() != expected {
+        return Err(PyValueError::new_err(format!(
+            "x has {} entries, expected n_persons * n_items = {expected}",
+            flat.len(),
+        )));
+    }
+    let rows: Vec<Vec<f64>> = (0..n_persons)
+        .map(|i| flat[i * n_items..(i + 1) * n_items].to_vec())
+        .collect();
+    let res = core_subscores(&rows, &groups).map_err(PyValueError::new_err)?;
+    let out = pyo3::types::PyDict::new(py);
+    out.set_item("alpha", res.alpha)?;
+    out.set_item("alpha_total", res.alpha_total)?;
+    out.set_item("prmse_s", res.prmse_s)?;
+    out.set_item("prmse_x", res.prmse_x)?;
+    out.set_item("prmse_sx", res.prmse_sx)?;
+    out.set_item("tau", res.tau)?;
+    out.set_item("beta", res.beta)?;
+    out.set_item("gamma", res.gamma)?;
+    out.set_item("added_value_s", res.added_value_s)?;
+    out.set_item("added_value_sx", res.added_value_sx)?;
+    out.set_item("total", res.total)?;
+    let flatten = |m: Vec<Vec<f64>>| -> Vec<f64> { m.into_iter().flatten().collect() };
+    out.set_item("corr", flatten(res.corr))?;
+    out.set_item("disattenuated_corr", flatten(res.disattenuated_corr))?;
+    out.set_item("observed", flatten(res.observed))?;
+    out.set_item("subscore_s", flatten(res.subscore_s))?;
+    out.set_item("subscore_x", flatten(res.subscore_x))?;
+    out.set_item("subscore_sx", flatten(res.subscore_sx))?;
+    Ok(out.into())
+}
 
 /// Marginal-EM fit of a mixed Rasch / mixture-IRT model (`mlsirm_core::mixture`, Rost,
 /// 1990). `y`/`observed` are row-major `n_persons * n_items`; `model` is "rasch" or
@@ -5112,6 +5232,8 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(fit_facets, m)?)?;
     m.add_function(wrap_pyfunction!(mokken_coef_h, m)?)?;
     m.add_function(wrap_pyfunction!(mokken_aisp, m)?)?;
+    m.add_function(wrap_pyfunction!(ksirt_occ, m)?)?;
+    m.add_function(wrap_pyfunction!(subscore_analysis, m)?)?;
     m.add_function(wrap_pyfunction!(fit_mixture, m)?)?;
     m.add_function(wrap_pyfunction!(fit_lltm, m)?)?;
     m.add_function(wrap_pyfunction!(fit_testlet, m)?)?;