Marginal (MMLE-EM) estimation for the latent-space family: multigroup/multilevel/FIPC, QMC-EM, GPU E-step, Rust scoring & fit statistics, paper-batch features

CHANGELOG.md

@@ -2,8 +2,89 @@
 
 ## Unreleased
 
+### Added
+
+- **Marginal (MMLE-EM) estimation for the full latent-space family.**
+  `fit(estimator="mmle")` now fits `MIRT`/`MLS2PLM`/`MLSRM` (and `ULS2PLM`/
+  `ULSRM` under a population structure) by Bock-Aitkin-style marginal EM:
+  person latents `(theta, xi)` are integrated over Gauss-Hermite grids —
+  tractable via the simple-structure conditional factorization — with a
+  Fisher-preconditioned GEM M-step and the Jeon et al. (2021) LSIRM priors as
+  MAP penalties (`PenaltyConfig::lsirm_prior`). Rust core
+  (`mlsirm_core::marginal`) with a NumPy mirror
+  (`fast_mlsirm.estimators.marginal`) held to 1e-9 end-of-run parity
+  (`tests/test_marginal_parity.py`); design and paper basis in
+  `docs/mmle_marginal_lsirm_design.md`.
+- **Estimation-level multigroup and multilevel population structures** for the
+  marginal estimator: `fit(..., group_id=...)` (Bock-Zimowski group trait
+  means/SDs, common items, pinned reference group) and
+  `fit(..., cluster_id=...)` (Fox-Glas random intercept, `sigma_u`/ICC
+  estimated). Results surface on `FitResult.population` and persist through
+  `save_fit_result`; the CLI `fit` command gains `--estimator`, `--group-id`,
+  `--cluster-id`, `--q-theta`, `--q-xi`, `--q-u`, and `--tolerance`.
+- **wgpu E-step kernels for the marginal estimator**
+  (`mlsirm_core::gpu_marginal`): the E-step hot path runs in f32 on the GPU
+  with the same race-free slot-ownership reduction as the JML kernels, cutting
+  a 31k-person multilevel E-step iteration from ~110 s (CPU f64) to ~5 s on a
+  laptop RTX 3050 Ti; the M-step and final EAP pass stay on the CPU in f64,
+  and hosts without an adapter fall back to the CPU path unchanged.
+- **Likelihood-based fit statistics** (`fast_mlsirm.fitstats`): Orlando-Thissen
+  S-X² via the Lord-Wingersky recursion generalized to the joint `(theta, xi)`
+  grid (chi-square tail without SciPy), Benjamini-Hochberg FDR control,
+  Drasgow `l_z` and Snijders `l_z*` person fit with the MAP `r_0` correction,
+  and infit/outfit at the marginal EAPs.
+- **Item screening pipeline** (`fast_mlsirm.select_items`): iterative
+  fit → flag → remove → refit with sparse / S-X²-BH / mean-square band /
+  low-discrimination / map-isolation flags, an `l_z*` person screen, a
+  per-dimension item floor, and a full audit trail.
+- **Serving bundle + frozen-parameter scoring** (`fast_mlsirm.serving`):
+  schema-versioned JSON bundle of the calibrated item parameters and
+  population block, and `score_respondents()` EAP scoring of new response
+  payloads with items frozen — the fixed-parameter serving pattern used by
+  the downstream importance-assessment API. `fast-mlsirm score` scores a JSON
+  payload (or `.npy` matrix) against a bundle from the command line.
+
+- **QMC-EM and MC-EM integration rules** for the marginal estimator
+  (`FitConfig(xi_rule="qmc"|"mc", xi_points=..., xi_seed=...)`): the
+  latent-space integral runs on Halton low-discrepancy points (randomized-QMC
+  shift optional; Jank 2005) or seeded Monte Carlo draws (Wei & Tanner 1990;
+  Meng & Schilling 1996) instead of the tensor Gauss-Hermite grid — enabling
+  `latent_dim > 3` and better error scaling per node. Both constructions are
+  deterministic and bit-mirrored across the Rust/NumPy backends.
+- **Rust scoring module** (`mlsirm_core::scoring`, exposed via
+  `_core.score_bank_eap` / `score_bank_map` / `eapsum_tables`): EAP
+  (Bock & Mislevy 1982), MAP (posterior Newton with observed-information
+  SEs), and summed-score EAP conversion tables via the Lord-Wingersky
+  recursion (Thissen et al. 1995; Cai 2015), all under per-dimension
+  `N(mean_d, sd_d^2)` priors that cover single, multigroup
+  (`mu_g, sigma_g`) and multilevel populations (conditional
+  `N(u_hat_c, 1)` or marginal `N(0, sqrt(1 + sigma_u^2))`).
+  `score_respondents(..., method="eap"|"map"|"eapsum", prior=...)` and the
+  bundle's embedded `eapsum_tables` expose these to serving.
+- **Fit statistics moved to the Rust core** (`mlsirm_core::fitstats`): S-X²,
+  Benjamini-Hochberg, `l_z`/`l_z*`, infit/outfit now compute in Rust
+  (`fast_mlsirm.fitstats` delegates; the NumPy bodies remain the parity
+  reference/fallback). S-X² gains the `rms_residual` practical-significance
+  effect size (Sinharay & Haberman 2014) and `select_items` gates its flag on
+  `sx2_min_effect`; the mean-square gate now uses infit only (outfit is
+  reported, not gating — it explodes under very low pass rates); the person
+  screen threshold is configurable and the Snijders `r_0` correction is
+  centered on the population prior mean (cluster intercepts / group means).
+- **Fixed Item Parameter Calibration** (`fit(..., anchors=...)`): anchored
+  items stay frozen (optionally `tau` too) while new items and a freed
+  population mean/SD are estimated — the multiple-cycle prior-update (MWU-MEM
+  style) variant Kim (2006) found robust; latent-space orientation inherits
+  from the anchors (no PCA re-alignment). **Concurrent calibration** is the
+  existing multigroup path with structural missingness (Hanson & Béguin
+  2002), covered by a dedicated recovery test.
+
 ### Changed
 
+- `estimator="mmle"` with a spatial/multidimensional model now fits (routed to
+  the marginal estimator) instead of raising `NotImplementedError`; plain
+  `ULS2PLM`/`ULSRM` without a population structure keep the legacy
+  unidimensional fast path and its exact previous behavior.
+
 - Exposed the Rust MMLE-EM estimator (`mlsirm_core::mmle::fit_mmle_2pl`) through
   the PyO3 binding as `fast_mlsirm._core.fit_mmle_2pl`, so
   `fit(estimator="mmle")` now runs on the Rust core when the extension is built