MMLE marginal LSIRM stack with GPU-first scoring and audit fixes - #213
Merged
Conversation
…ltigroup/multilevel structures and a wgpu E-step Model-design PR (paper basis in docs/mmle_marginal_lsirm_design.md and docs/papers/mmle-lsirm-formula-compilation.md): - crates/mlsirm-core/src/marginal.rs — deterministic Bock-Aitkin-style marginal EM for MIRT/MLS2PLM/MLSRM/ULS2PLM/ULSRM. Person latents (theta, xi) integrated over Gauss-Hermite grids, feasible via the simple-structure conditional factorization (Q_xi^K * sum_d Q_theta, not Q^(1+D+K)). Fisher-preconditioned GEM M-step with Armijo line search; MAP penalties default to the Jeon et al. (2021) priors (PenaltyConfig::lsirm_prior). Multigroup (Bock-Zimowski) group means/SDs with pinned reference group; multilevel (Fox-Glas) random intercept with estimated sigma_u. PCA-aligned zeta/xi for the rotation/reflection invariance. - crates/mlsirm-core/src/quadrature.rs — embedded hermegauss tables (7..41 nodes), bit-identical to the NumPy reference. - crates/mlsirm-core/src/gpu_marginal.rs — wgpu f32 E-step kernels (lp/nbar/item passes, race-free slot ownership); M-step and final EAP stay CPU f64. Measured 110s -> ~5s per multilevel E-step iteration on a 31k x 57 dataset (RTX 3050 Ti). - python/fast_mlsirm/estimators/marginal.py — NumPy mirror; parity with the Rust core at 1e-9 after full EM runs (tests/test_marginal_parity.py). - fit(group_id=..., cluster_id=...) drives the population structures; FitResult.population carries mu/sigma/sigma_u/u_eap/icc/theta_sd and save_fit_result persists them. CLI fit gains --estimator/--group-id/ --cluster-id/--q-theta/--q-xi/--q-u/--tolerance. - Plain ULS2PLM/ULSRM keep the legacy fast path (unchanged behavior); spatial models under estimator="mmle" now fit instead of raising. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…nd serving bundle - python/fast_mlsirm/fitstats.py — Orlando-Thissen S-X2 with the Lord-Wingersky recursion generalized to the joint (theta, xi) quadrature grid (per-dimension summed scores, expected-count collapsing, chi-square p via a no-SciPy regularized upper incomplete gamma), Benjamini-Hochberg FDR, Drasgow l_z and Snijders l_z* with the MAP r_0 correction at EAP estimates, infit/outfit at the marginal EAPs, and select_items(): the literature-grounded fit -> flag -> remove -> refit loop (sparse / S-X2-BH / MSQ band / low discrimination / map isolation flags, person-fit screen, per-dimension item floor, full audit trail). - python/fast_mlsirm/serving.py — schema-versioned JSON serving bundle (frozen item parameters + population block + screening audit) and score_respondents(): EAP scoring of new response payloads (dict or dense) against the frozen bundle, the same fixed-parameter pattern as the downstream importance-assessment API. - estimators/marginal.py gains score_eap() (one E-step pass, no updates); public API re-exports; tests for chi2_sf/BH/Lord-Wingersky against enumeration, S-X2 flagging a scrambled item, l_z* calibration and aberrant-person detection, screening pipeline behavior, and bundle round-trip/scoring monotonicity. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…set; add `fast-mlsirm score` `fast-mlsirm score --bundle b.json --responses r.json` scores new respondents against a frozen serving bundle (JSON code->0/1 payloads or a .npy matrix), mirroring the downstream fixed-parameter serving flow. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… FIPC and concurrent calibration All numeric compute now lives in the Rust core; Python keeps thin wrappers plus the NumPy parity references. Paper basis: Part II of docs/papers/mmle-lsirm-formula-compilation.md (Wei-Tanner 1990, Booth-Hobert 1999, Jank 2005, Meng-Schilling 1996, Bock-Mislevy 1982, Thissen et al. 1995, Lord-Wingersky 1984 via Cai 2015, Kim-Cohen 1998, Hanson-Beguin 2002, Kim 2006, Sinharay-Haberman 2014). - mlsirm-core/nodes.rs: shared latent-space node sets — tensor GH, Halton QMC (+ Cranley-Patterson shift), seeded MC; Acklam inverse normal CDF. Bit- mirrored in the NumPy reference (parity ~1e-15 after full EM runs). - mlsirm-core/scoring.rs: ItemBank scoring — EAP, damped-Newton MAP with observed-information SEs, EAPsum tables via Lord-Wingersky; per-dimension N(mean, sd^2) priors cover single/multigroup/multilevel serving. - mlsirm-core/fitstats.rs: S-X2 (+ rms_residual effect size), BH, lz/lz*, infit/outfit; chi-square tail via regularized upper incomplete gamma. - marginal.rs: xi_rule/xi_points/xi_seed; Anchors (FIPC) with frozen items and optional frozen tau; PopulationSpec::SingleFree (free mu/sigma, anchor-identified); GPU guard for oversized QMC node sets. - select_items methodology hardening after the first real-data run over-pruned 45/57 items: S-X2 flag now requires a practical effect size, the MSQ gate uses infit only, the person screen is threshold-configurable and prior-mean centered. - fast-mlsirm score/serving: method="eap"|"map"|"eapsum", prior override for known-team/known-group conditioning, eapsum_tables embedded in the bundle. - Tests: 263 pytest + 43 cargo, incl. QMC/MC backend parity at 1e-9, FIPC anchor freezing + population recovery, concurrent-calibration recovery under structural missingness, EAPsum monotonicity, MAP SEs. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…dation gates, DIF, Vuong, Q3/GDDM, IRTree expansion, information criteria
Implements the implementable core of the supplied literature set:
- Perumean-Chaney et al. (2013) -> zero-inflated marginal mixture
(FitConfig(zero_inflation=True)): structural-zero class with EM-estimated
pi, responsibilities folded into the E-step weights (GPU kernels
untouched); works under single/multigroup/multilevel populations.
- Debeer & Janssen (2013) -> context-varying item covariate with one
estimated coefficient (fit(covariate={"w", "init_delta"})): the linear
item-position effect with booklet groups as contexts; Newton coordinate in
the M-step; identification guards (single-context covariate rejected).
- Williamson, Xi & Breyer (2012) -> mlsirm_core::agreement + validate_judge:
QWK/kappa, Pearson r, SMD, degradation-vs-human-human and subgroup-SMD
conjunctive acceptance gates with the paper's thresholds.
- Jeon, Rijmen & Rabe-Hesketh (2013) + Makransky & Glas (2013) ->
dif_analysis(): LR DIF screen via group-specific virtual items with all
other items anchored, BH-FDR over studied items, logit effect sizes.
- Schneider, Chalmers, Debelak & Merkle (2019) -> vuong_nonnested(): Vuong
z from casewise logliks with optional Schwarz correction (erfc-based
normal tail in the core).
- Svetina & Levy (2014) -> dimensionality_residuals(): Yen Q3 and GDDM from
EAP residuals in the Rust core.
- Jeon & De Boeck (2016) -> irtree_expand(): mapping-matrix pseudo-item
expansion (their Eq. 9 equivalence makes IRTrees ordinary binary IRT on
the expanded matrix; off-path nodes reuse NaN missingness).
- Kang, Cohen & Sung (2009) -> information_criteria in the core
(AIC/BIC/AICc/SABIC/CAIC + free-parameter counting incl. anchors/ZI/
covariate), surfaced as FitResult.ic and in fit_summary.json.
ZI and covariate are mirrored in the NumPy reference (parity at 1e-9);
Wolkowitz & Skorupski (2013) documented as superseded by marginal-ML MAR
handling; Ferrando et al. (2009) and Joubert et al. (2015) documented as
out-of-scope for binary judge data (specs in docs-research/group_*_specs.md).
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…in 2017) mlsirm_core::oakes::observed_information_oakes assembles the observed (penalized) information at the MML solution as the M-step Hessian at the fixed posterior (central FD of the analytic Q-gradient, no E-steps) plus the Oakes cross term (one CPU f64 E-step per parameter) — the estimator Pritikin (2017) recommends over the supplemented-EM family. SEs cover the item-side parameters and tau, conditional on the population parameters; exposed as fast_mlsirm.oakes_standard_errors(result, responses, factor_id, config). Internal-consistency test: the Oakes assembly matches the full central difference of the marginal score on probe coordinates. Serving bundles now carry pi_zero / covariate_delta in the population block. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… specs Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…irwise fit, resampling person fit, TCC drift; corpus triage map Implemented from the third supplied corpus (triage in docs/papers/corpus-triage-batch3.md): - Magis (2013): 4PL item information (reduces to 2PL at c=0/d=1); bank_information() item/test information at arbitrary trait points. - Bock & Mislevy (1982) + Wang, Kuo & Chao (2010): cat_next_item() adaptive EAP step over a frozen serving bundle — targets the trait dimension with the largest posterior SD, ranks unadministered items by information. - Marsman et al. (2016): plausible_values() seeded posterior draws from the scoring grid for secondary analyses. - Haberman, Sinharay & Chon (2013): residual_item_fit() standardized EAP-bin residual fit (long-test regime documented; S-X2 for short tests). - Tay & Drasgow (2012): adjusted_chi2_pairs() N=3000-adjusted pairwise chi2/df ratios for local-dependence screening. - Sinharay (2016): person_fit_resampling() parametric-bootstrap empirical p-values for l_z* at the EAP estimates. - Guo, Zheng & Chang (2015): tcc_drift() stepwise TCC drift detection between two same-scale calibrations. All compute in mlsirm-core (scoring.rs / fitstats.rs) with PyO3 wrappers and bundle-level Python APIs; 3PL/4PL estimation and the BIFAC2PLM bifactor variant are scoped as the next model-design PRs in the triage map, and the polytomous/CDM/Bayesian clusters are dispositioned there with reasons. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… corpus triage ld_indices() computes the signed standardized pairwise LD X2/G2 against the model-implied joint probabilities on the scoring grid — the classic local-dependence screen, complementing Q3/GDDM and the Tay-Drasgow adjusted ratios. docs/papers/corpus-triage-batch4.md dispositions the fourth reading set and consolidates the explicitly-requested roadmap (BIFAC2PLM bifactor, M2/RMSEA2, 3PL/4PL estimation, polytomous kernels, linking utilities, response times) in priority order. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…zi et al. 2015); batch-5 corpus triage empirical_reliability() reports the marginal Var(EAP)/(Var(EAP)+mean(SE^2)) per trait dimension for marginal fits, with the sources' caveat that the coefficient presumes a well-fitting model. corpus-triage-batch5.md dispositions the fifth reading set and settles the BIFAC2PLM design as an inner-product interaction kind on the existing eta plumbing (tables/GPU kernels unchanged) — the next model-design PR, now requested across three batches, followed by M2, general mixtures, 3PL/4PL, and polytomous kernels. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Implements the full-information dichotomous bifactor model (Gibbons &
Hedeker 1992; Cai, Yang & Hansen 2011) — requested across three literature
batches — as a third interaction kind on the marginal engine:
- InteractionKind::{None, Distance, Inner} parameterizes every eta site
(tables, M-step gradients with d eta/d zeta_k = x_k, scoring EAP/MAP/
information, fit statistics, Oakes SEs); tau is inert for the inner kind
and excluded from the parameter vector/counts.
- BIFAC2PLM = inner kind at latent_dim >= 1: the general factor rides the
existing conditional-factorization E-step (the Gibbons-Hedeker dimension
reduction), so tables and the wgpu kernels carry over unchanged; at
latent_dim = 1, zeta_i is the general-factor loading lambda_i.
- Positive-manifold loading initialization (a mixed-sign circle init locks
sign-split local optima); marginal estimator only — the JML path guards.
- Recovery tests (Rust + Python) and 1e-9 Rust/NumPy parity; empirical
reliability, EAPsum, CAT, plausible values, and the fit-statistic stack
all operate on BIFAC2PLM banks through the shared kind dispatch.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add M2 (Maydeu-Olivares & Joe 2005/2006; Cai & Hansen 2013) on the univariate + bivariate residual margins: statistic, df, chi-square p-value, the RMSEA2 approximate-fit index with a 90% noncentral-chi-square confidence interval, and the bivariate SRMSR (Maydeu-Olivares 2013). Every model-implied margin and the up-to-4th-order entries of the multinomial residual covariance Xi_2 are exact via the local-independence factorization over the (theta, xi) node set -- pi_S = sum_c w_c prod P_i(c) -- the same factorization the E-step uses (Cai-Hansen dimension reduction). Delta_2 is central-differenced from the node moments; the quadratic form is evaluated through one Cholesky of Xi_2 (never an explicit inverse): M2 = N ( e'Xi^-1 e - g'(D'Xi^-1 D)^-1 g ), g = D'Xi^-1 e. Kind-aware Rust core (mlsirm_core::fitstats::m2_rmsea2) is the compute path; a NumPy reference (fast_mlsirm.fitstats.m2 / _m2_numpy) is held to 1e-6 parity. Calibration tests in both suites contrast a well-specified fit (RMSEA2 < 0.03) against injected local dependence (RMSEA2 > 0.08). Batch-6 corpus triage; M2/RMSEA2 moves from the roadmap to done. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add mlsirm_core::linking (+ fast_mlsirm.irt_link): put a separately- calibrated new form onto the reference scale from common items, returning theta_old = A*theta_new + B. Covers the moment methods (mean/mean, mean/sigma) and the characteristic-curve methods of Haebara (1980) and Stocking & Lord (1983) for the unidimensional 2PL/1PL common-item case those procedures reduce to (Kolen & Brennan 2014). New-form items transform onto the old scale in the engine's eta = a*theta+b form as a* = a_new/A, b* = b_new - (a_new/A)*B; the characteristic-curve loss is minimized by a self-contained Nelder-Mead from the mean/sigma start, integrated over a standard-normal Gauss-Hermite grid. Rust compute path; four-method recovery tests (Rust + Python) recover a known transform. Motivated by the corpus linking papers (Kim & Lee 2006; Yao & Boughton 2009; Brossman & Lee 2013). Complements link_fixed_item_parameters and the FIPC serving path. Batch-6 triage updated. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add branch tests for the new code so its error/edge paths are exercised under the lib test suite (not only the integration recovery test): - linking: LinkMethod::parse (all arms + unknown), the moment zero-spread fallback, the cc-objective non-positive/NaN guard, Nelder-Mead contraction/shrink on a non-smooth objective, and the irt_link input guards (bad slopes, empty/mismatched grid). - M2: the too-few-items / length-mismatch / non-positive-df / too-few-complete-cases guards, plus a small hand-built-bank run that exercises the Cholesky, Delta, Xi, CI, and SRMSR body directly. Docstrings added to the new Python helpers (M2Result, _MutBank, _m2_numpy, _ncchi2_cdf, _nc_lambda_for, n_dims_of, IrtLinkResult) for the review gate's docstring-evidence check. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Address the Strix scan findings on PR #160 — DoS / data-poisoning guards for serving/fitting entry points that may take untrusted data: - VULN-0001: compact + validate population labels (group_id/cluster_id) in fit.py and inference.py. n_groups/n_clusters is now the distinct-label count (<= n_persons), not max+1, so sparse ids like [0, 1e9] can't force billion-row allocations; negative/non-integer/non-finite/wrong-length ids are rejected. - VULN-0004/0005: FitConfig.validate bounds latent_dim (<= 8) and xi_points (<= 1_000_000). - VULN-0006/0007: load_serving_bundle parses JSON strictly (no NaN/Infinity) and runs _validate_bundle (consistent, bounded dimensions; in-range factor_id; finite alpha/b/zeta/tau/eps_distance; supported quadrature); score_respondents and plausible_values validate the bundle at entry. - VULN-0002: plausible_values enforces the 0/1 finite response domain that score_respondents already required. - VULN-0003: validate_judge validates labels (1-D, equal length, finite, integer, 0<=label<k) before the uint32 conversion. Regression tests in tests/test_security_hardening.py (19) cover every finding; existing serving/scoring/estimator suites unaffected. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add unit tests for previously-uncovered defensive branches toward the 100%-line coverage gate: gh_rule None (unsupported size) and the Halton high-latent-dim / shift paths in nodes; XiRuleKind::parse (all arms + unknown) in marginal; validate_bank's four error branches and the y/observed length guard via score_eap in scoring. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Second Strix pass on PR #160. Close the remaining input-validation gaps (the earlier pass bounded only latent_dim/xi_points): - VULN-0004: FitConfig.validate now rejects non-finite learning_rate, init_gamma, eps_distance, tolerance, and gradient_clip (a bare `x <= 0` lets NaN/Inf through), and bounds max_iter (<=100_000), n_restarts (<=1_000), and m_steps (<=1_000) so oversized loops/allocations and NaN-poisoned fits are rejected up front. - VULN-0005: plausible_values bounds n_draws (1..100_000) before forwarding to the core; serving_prior bounds n_dims (1..64) for direct callers. VULN-0001/0002 (non-finite params, malformed bundle) were already closed by the prior commit's strict JSON parse + _validate_bundle; VULN-0003 (missing-module import failure) is a false positive from the scanner's PR-diff-only checkout (backend.py/diagnostics.py/report.py all exist and the package imports). Regression tests extended to 32. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…ning
The security-hardening bounds changed several FitConfig.validate messages,
breaking test_config.py's message-substring assertions (the values were
still correctly rejected). Reword the new messages as supersets of the old
text ("... >= 1 and <= N", "... > 0 and finite") so existing tests pass and
the stricter semantics remain.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Offload Bock-Mislevy (1982) EAP scoring to the wgpu path — the 31k-person serving hot path — per the all-math-in-Rust / maximize-GPU policy. - WGSL score_pass: one thread per person, race-free (each person owns its theta_eap/theta_sd/xi_eap/loglik slots; no atomics/slot-ownership unlike the E-step). Reuses the cell_l binary-sparsity decomposition and the same logp0/logp1/c0 tables the CPU scoring builds. f32, so ~1e-4 vs the f64 CPU reduction. Separate 19-binding bind-group layout + pipeline. - score_eap now delegates to a new score_eap_device(..., device): Cpu keeps the exact f64 reduction (the default -- all precision-sensitive callers and serving parity unchanged); Gpu/Auto try score_eap_gpu and fall back to CPU when no adapter or n_dims/latent_dim > 8. The CPU reduction was extracted to score_eap_cpu_reduce (the parity reference). - PyO3 score_bank_eap gains device="cpu"; serving.score_respondents threads device through so serving can request the GPU. On-device parity test gpu_eap_matches_cpu_reduction (<=2e-3) and a PyO3 device=gpu-vs-cpu smoke both pass on the local RTX; lib suite 67 green, serving/scoring pytest green. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…re surface Extends the input-validation / DoS guards to the code added in this PR: - preprocessing.irtree_expand: bound dense expansion (persons*items*nodes <= 50M) before allocating; validate node_dims (finite, non-negative, integer) before int64 cast (VULN-0001/0002) - validation: reject labels above uint32 max before the narrowing cast; require human_human baseline to match the paired sample size (VULN-0003/0011) - inference.observed_information: cap finite-difference Hessian at 5000 params (O(n^2) memory and objective calls); oakes_standard_errors validates factor_id (1-D, per-item, finite, non-negative, integer) before deriving n_dims (VULN-0004/0005) - serving._validate_bundle + estimators.marginal._xi_grid: reject tensor grids with q_xi ** latent_dim > 1M points (VULN-0006) - linking.link_fixed_item_parameters: reject duplicate/fractional/ negative/non-finite anchors, non-2-D theta, non-finite item params, and non-finite computed linking coefficients (VULN-0007..0010) Regression tests in tests/test_security_hardening.py (53 pass); full suite 325 pass. All 11 Strix PoCs verified blocked. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…t (6 more) Strix re-scan of b5d9d90 reported 12 findings; #1 ("incomplete package release") is a scanner artifact of its PR-scope-only checkout (every named module exists; `import fast_mlsirm` succeeds). The 11 real findings: - serving.score_respondents/plausible_values: bound the dense respondent matrix (rows x n_items) before np.full (VULN-0002) - linking: range-check anchor indices on the float BEFORE the int64 cast (uint64 max wrapped to -1 -> last-item index) + same for factor_id (0003) - validation.validate_judge: bound category count k (dense k x k core matrix) (0004) - preprocessing.irtree_expand: 50M-element ceiling (400 MB, inclusive) -> 64 MiB byte budget (0005) - config.MLS2PLMConfig.validate: bound sim dims + n_persons*n_items cells (0006); FitConfig.validate: bound aggregate max_iter*n_restarts (0008) - estimators.marginal.fit_marginal_numpy: bound population counts (n_groups/n_clusters <= n_persons) (0007) + EM working set (0012) - inference.observed_information: reject non-finite step (0009); oakes_standard_errors: validate factor_id (0010, extended below) - fitstats: shared _validate_factor_id bounds n_dims for s_x2/person_fit/ infit_outfit and all public entries (0010) Proactive boundary-audit workflow found 6 more Strix had not surfaced: - serving._validate_bundle: bound scoring-table product (max(items,dims) x q_theta x q_xi**latent_dim; 55+ GB otherwise) and validate the population block (serving_prior read an unvalidated, attacker-controlled sigma_u -> TypeError/OverflowError crash or silent Inf/NaN score poisoning) - linking.irt_link: validate slope/intercept finiteness + slope positivity before the Nelder-Mead core (NaN would panic it); link_fixed_item_parameters requires a positive linking scale - validation.validate_judge: compact sparse subgroup labels (core loops 0..max(label)+1 -> O(4e9) CPU-DoS); oakes: reject n_dims > n_items Regression tests in tests/test_security_hardening.py (all 17+ new cases, suite 76 file-local); full suite 348 pass; every PoC verified blocked. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Defense-in-depth for the confirmed audit findings whose reachable paths are now guarded at the Python layer, but which panic if the core is called directly (PyO3 / Rust callers): - fitstats::s_x2 rejects non-dichotomous observed responses — a non-0/1 value made the summed-score index (sized n_d+1) go out of bounds - fitstats::infit_outfit validates theta/xi lengths before indexing theta[p*n_dims+d] / xi[p*latent_dim+k] - scoring::validate_prior rejects non-finite prior mean/sd (a NaN sd passed the bare `sd <= 0` check, poisoning quadrature with Inf/NaN) Rust regression tests (sx2_rejects_non_dichotomous_responses, infit_outfit_rejects_wrong_theta_length); cargo test -p mlsirm-core green (69 lib tests); full Python suite 348 pass against the rebuilt extension. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…design PR) Synthesized by the gpcm-design workflow (study -> 3 independent designs -> judge panel -> merge). Unified scoring-function softmax cell nesting binary 2PL / GPCM (Muraki) / nominal (Bock) in one kernel that reuses eta_at_kind, with a ResponseModel::Bernoulli branch preserving binary bit-parity. Records the forced modeling consequence (the latent-space term enters category-scaled for identification, so the polytomous-LSIRM space axis is ordered, not nominal) for maintainer awareness. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
First step of the GPCM/nominal model-design PR (docs/papers/ gpcm-nominal-design-spec.md). `category_logprobs(base, scores, intercepts)` is the stable log-softmax of psi_k = scores[k]*base + intercepts[k] with the baseline category pinned — the parity reference the Rust cell will be held to. Nests binary 2PL exactly (logP_1 = log_sigmoid(base+b), the design's free bit-parity check), GPCM (fixed integer scores), and nominal (free scores). Design-neutral and pure Python: does not touch the working E-step core. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Second GPCM foundation piece. gpcm_node_gradient returns the category residual (r_m - n*P_m) for intercepts, the score-weighted residual R = sum_k s_k*resid_k for the base coordinate (which reuses the binary deta terms verbatim), and the nominal-score gradient resid_m*base. A central finite-difference test confirms all three components — de-risking the Rust M-step math before it is written. Still pure Python; the working E-step core is untouched. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Third GPCM foundation piece: fit_gpcm_numpy, a unidimensional GPCM (Muraki)
marginal MLE via Bock-Aitkin EM, exercising the verified cell
(category_logprobs) and residual gradient (gpcm_node_gradient) in a full EM
loop with a Newton M-step. Recovers known slopes (corr > 0.9) and category
intercepts from simulated data. Reports Muraki step thresholds
b_{i,k}=c_{i,k-1}-c_{i,k}. This is the parity oracle the Rust polytomous
kernel will mirror; still pure Python, working core untouched.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Literature search (alphaXiv + Consensus) on how the latent-space interaction combines with polytomous categories. Jeon et al. 2021 (base LSIRM) and Go et al. 2024/2025 (lsirm12pl, the authors' own package) both add the interaction as a single scalar on the linear predictor and route non-binary responses through "a suitable link function, as in GLMs" — not a softmax. No published ordinal/polytomous LSIRM exists yet. Resolution: the cumulative-logit Graded Response Model (Samejima) is the identification-clean default (single shared -gamma*d shift across thresholds, no forced category-scaling); GPCM/nominal softmax kept as documented options. Appends the resolution to the design spec. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Per the all-math-in-Rust directive, the polytomous compute lives in mlsirm-core::poly, not Python. Implements both response families over the shared linear predictor base = a*theta + interaction(x): - GRM (Samejima, cumulative logit) — the identification-clean default for the LSIRM family (single shared -gamma*d shift across thresholds, no forced category scaling; see the design spec's literature resolution): grm_logprobs + grm_node_gradient. - GPCM/nominal (Muraki/Bock, softmax) — partial-credit option: gpcm_logprobs + gpcm_node_gradient (category residual / score-weighted / nominal-score grads). Rust unit tests: cell normalization, binary bit-parity (logP_1 = log_sigmoid), GPCM monotonicity, and central finite-difference checks of every gradient component (4 pass). The NumPy category_logprobs/gpcm_node_gradient/fit_gpcm_numpy are demoted to parity references (their proper role, like fit_marginal_numpy). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Replace dataset-specific references in the spec with generic response-format wording (the model default is determined by the data's response format: ordinal/rubric -> GRM, partial-credit -> GPCM). No behavior change. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Apply observed-aware category validation to the remaining polytomous multigroup, S-X2, and U3 paths so masked sentinel values are ignored consistently. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
GitHub dependency-graph API returned HTTP 502 on last run; no code changes -- re-running to get a clean result. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
Implement a reduced-scope response-bias adaptation for binary diagnostics using Rust compute kernels and thin Python wiring. - add core leniency residual statistics in mlsirm-core fitstats - expose via PyO3 as leniency_residuals_stat - publish person/model-level leniency metrics in fit_diagnostics - add mutation-oriented mask/sign regression tests (Rust + Python) - update implemented-literature-map status to adapted Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
Set person-level leniency outputs to finite defaults when a respondent has no observed cells, while keeping summary statistics computed only on valid rows. This restores the fit diagnostics finite-value contract used by stability tests. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
…cre, 1989) (#218) * feat(facets): many-facet Rasch model rater-severity calibration Add mlsirm_core::facets implementing the MFRM (Linacre, 1989; Eckes, 2015): adjacent-category log-odds theta - d_i - c_j - f_k, i.e. the rating scale model (Andrich, 1978) with a rater facet, estimated by marginal-ML EM on a Gauss-Hermite grid (Bock & Aitkin, 1981). This is MMLE, not Linacre's JMLE, and the docs state the Facets-comparability caveat. Identification: theta ~ N(0,1), severities and thresholds centered (n_parameters = I + (J-1) + (K-2)); each EM cycle re-absorbs centering shifts into item difficulty, which is likelihood-invariant. Reports Linacre's connectedness requirement via union-find over the person-mediated item/rater co-observation graph; connected=false means cross-component comparisons rest solely on the trait prior. Rust-only numerics reusing rsm_logprobs and solve_small; PyO3 fit_facets binding; thin validating Python wrapper fast_mlsirm.fit_facets over a persons x items x raters NaN-missing array. Tests: FD gradient anchors for locations and thresholds, J=1 reduction to fit_rsm, asymmetric severity recovery, sparse-design recovery, disconnected + bridged connectivity, input rejection, monotone loglik trace, and an #[ignore] 500-rep Monte Carlo (normal and skew-normal traits) bounding severity bias/RMSE. A gradient sign-flip mutant was verified to fail 4 tests. Adversarial spec review and implementation review completed; the one confirmed defect (a false mutation-kill claim in a test docstring) is fixed by a positive connectivity assert. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> * Treat negative sentinels as missing in fit_facets (review feedback) The repo-wide missing-response convention is NaN, -1 (negative sentinels), or an explicit mask; the facets wrapper previously accepted only NaN and rejected negatives with a ValueError. The observed mask now excludes negative cells before marshaling to the Rust core (which already honors the mask), docstrings document the convention, and a regression test asserts -1-coded and NaN-coded missing cells produce identical fits. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --------- Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
seonghobae
added a commit
that referenced
this pull request
Jul 25, 2026
Audit fixes identified during adversarial review of PR #213: - poly.rs: fix degenerate GRM threshold handling - fitstats.rs: correctness improvements - poly_marginal.rs: updated tests - polytomous.py: minor corrections Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
seonghobae
added a commit
that referenced
this pull request
Jul 25, 2026
…240) Audit fixes identified during adversarial review of PR #213: - poly.rs: fix degenerate GRM threshold handling - fitstats.rs: correctness improvements - poly_marginal.rs: updated tests - polytomous.py: minor corrections Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Why
This branch brings the full MMLE-marginal LSIRM feature stack into a single reviewable PR and adds a final audit pass for literature traceability and runtime performance. The audit was needed to verify commit #160 claims against local Zotero evidence and to tighten CPU fallback throughput where GPU is unavailable.
What changed
Device::Auto/device="auto") with explicit CPU fallback behavior.crates/mlsirm-core/src/scoring.rsby replacing serial reductions with fixed-shard parallel reductions for:score_eap_cpu_reducebank_information_cpu_reduceThis reduces context-switch overhead by using bounded worker sharding instead of fine-grained task dispatch.
docs/papers/implemented-literature-map.mdso it points to actual in-repo paper spec files underdocs/papers/.Zotero local API audit notes
http://localhost:23119/apithat referenced psychometric papers with local attachments are resolvable, including PDF-backed entries (for example, Drasgow et al. 1985 and generalized S-X2 material)./items/{key}/fileredirection to local PDF storage for available attachments.Validation
cargo test --manifest-path crates\mlsirm-core\Cargo.toml scoring::