feat(mokken): Mokken scale analysis with Loevinger H and AISP

CHANGELOG.md

@@ -114,6 +114,26 @@
 
 ### Added
 
+- **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
+  Mokken Z statistics, and partitions items into Mokken scales with the
+  automated item selection procedure (AISP, "search normal"), with sample
+  statistics and selection mechanics verified line-by-line against the mokken
+  R package source (van der Ark, 2007; Straat et al., 2013): `Hij =
+  S_ij/Smax_ij` with `Smax` from the comonotone (sorted-column) coupling,
+  `Hi`/`H` as ratios of pairwise sums, and per-scale Bonferroni-adjusted Z
+  gates. For LLM-as-a-Judge item-quality management this flags evaluation
+  items that fail to scale (label 0) and detects multidimensional item pools
+  before parametric calibration. Complete integer data required (dichotomous
+  or polytomous). Rust-only numerics; the Python wrapper validates and
+  marshals. Tests include a brute-force covariance oracle, an exact Guttman
+  `H = 1` anchor, a hand-computed Z fixture, a Z-gate anchor whose deletion
+  seeds a spurious scale (this test caught a real sign error in the normal
+  quantile during development), a hand-constructed Criterion-1 design whose
+  negative-`Hij` exclusion is the only active gate (mutation-verified), a
+  two-cluster AISP recovery, and an `#[ignore]` 500-replicate Monte Carlo
+  (normal + skewed traits).
 - **Many-Facet Rasch Model (MFRM) rater-severity calibration** (`fast_mlsirm.fit_facets`;
   new `mlsirm_core::facets`; Linacre, 1989; Eckes, 2015). Fits
   `ln[P(k)/P(k-1)] = theta_p - d_i - c_j - f_k` — the rating scale model

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

@@ -60,6 +60,7 @@ 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::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::{
     fit_rt_lognormal as core_fit_rt, rt_person_fit as core_rt_person_fit, RtConfig,
@@ -1470,6 +1471,49 @@ fn fit_facets(
     Ok(out.into())
 }
 
+/// Mokken scalability coefficients (`mlsirm_core::mokken::coef_h`).
+/// `x` is a row-major complete `n_persons * n_items` integer score matrix.
+/// Returns a dict with `hij`/`zij` (flattened `J*J`, NaN diagonal), `hi`,
+/// `zi` (`J`), and scalars `h`, `z`. Sample statistics follow the mokken R
+/// package (van der Ark, 2007, https://doi.org/10.18637/jss.v020.i11).
+#[pyfunction]
+#[pyo3(signature = (x, n_persons, n_items))]
+fn mokken_coef_h(
+    py: Python<'_>,
+    x: PyReadonlyArray1<'_, i64>,
+    n_persons: usize,
+    n_items: usize,
+) -> PyResult<Py<pyo3::types::PyDict>> {
+    let res = core_mokken_coef_h(x.as_slice()?, n_persons, n_items)
+        .map_err(PyValueError::new_err)?;
+    let out = pyo3::types::PyDict::new(py);
+    out.set_item("hij", res.hij)?;
+    out.set_item("hi", res.hi)?;
+    out.set_item("h", res.h)?;
+    out.set_item("zij", res.zij)?;
+    out.set_item("zi", res.zi)?;
+    out.set_item("z", res.z)?;
+    Ok(out.into())
+}
+
+/// Mokken automated item selection procedure (`mlsirm_core::mokken::aisp`,
+/// the "search normal" algorithm of the mokken R package). Returns per-item
+/// scale labels: 0 = unscalable, 1, 2, ... in formation order. `c` is the
+/// scalability lower bound (rule of thumb 0.3), `alpha` the nominal
+/// significance level.
+#[pyfunction]
+#[pyo3(signature = (x, n_persons, n_items, c = 0.3, alpha = 0.05))]
+fn mokken_aisp(
+    x: PyReadonlyArray1<'_, i64>,
+    n_persons: usize,
+    n_items: usize,
+    c: f64,
+    alpha: f64,
+) -> PyResult<Vec<u32>> {
+    core_mokken_aisp(x.as_slice()?, n_persons, n_items, c, alpha).map_err(PyValueError::new_err)
+}
+
+
 /// 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
 /// "2pl". `n_classes` latent classes each get their own item parameters. Returns a dict
@@ -5066,6 +5110,8 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(fit_crm, m)?)?;
     m.add_function(wrap_pyfunction!(fit_rsm, m)?)?;
     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!(fit_mixture, m)?)?;
     m.add_function(wrap_pyfunction!(fit_lltm, m)?)?;
     m.add_function(wrap_pyfunction!(fit_testlet, m)?)?;

crates/mlsirm-core/src/lib.rs

@@ -14,6 +14,7 @@ pub mod mhrm;
 pub mod mixed;
 pub mod mixture;
 pub mod mmle;
+pub mod mokken;
 pub mod nodes;
 pub mod nominal;
 pub mod oakes;