Marginal (MMLE-EM) estimation for the latent-space family: multigroup/multilevel/FIPC, QMC-EM, GPU E-step, Rust scoring & fit statistics, paper-batch features - #160
Closed
seonghobae wants to merge 210 commits into
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>
![opencode-agent[bot]](https://avatars.githubusercontent.com/in/1549082?s=60&v=4){kind=small}
Contributor
There was a problem hiding this comment.
Pull request overview
OpenCode cannot approve yet because required coverage evidence did not pass.
Review outcome
1. HIGH .github/workflows/opencode-review.yml:1 - Coverage evidence did not prove required test/docstring evidence
-
Problem: The required coverage-evidence job result was
failure, so OpenCode cannot establish approval sufficiency for this head. -
Root cause: Automated approval is only valid when the same-head coverage-evidence job proves supported repository test suites passed and configured docstring gates passed or were advisory, or reports not applicable because no supported source files or package manifests exist. Missing, failed, skipped, unavailable, or unsupported-tooling test evidence is a blocker.
-
Fix: Install or configure the repository test/docstring evidence tooling when source files or package manifests exist, rerun the current-head coverage-evidence job, and approve only after it reports
successwith required evidence or explicit no-source not-applicable evidence. -
Regression test: Keep the approval branch checking
needs.coverage-evidence.result == successbefore posting APPROVE, and publish REQUEST_CHANGES when coverage-evidence blocker states such as cancelled, skipped, failed, unsupported-tooling, or below-100 evidence are present. -
Result: REQUEST_CHANGES
-
Reason: coverage-evidence result was
failure, so required test/docstring evidence was not proven for current headd1311b2ccc0db7233f20e578e8af9eff15a6280d. -
Head SHA:
d1311b2ccc0db7233f20e578e8af9eff15a6280d -
Workflow run: 29311144089
-
Workflow attempt: 1
Coverage evidence
Coverage Decision
- Result: FAIL
- Test evidence: not proven passing
- Docstring evidence: not proven passing when configured
- Failure count: 2
Changed-File Evidence Map
flowchart LR
PR["PR changed files"] --> Evidence["OpenCode bounded evidence"]
Evidence --> S1["Changed file (25 files)"]
S1 --> I1["repository behavior"]
I1 --> R1["Review risk: Changed file (25 files)"]
R1 --> V1["required checks"]
Evidence --> S2["Docs (8 files)"]
S2 --> I2["operator or user guidance"]
I2 --> R2["Review risk: Docs (8 files)"]
R2 --> V2["docs review"]
Evidence --> S3["Test (7 files)"]
S3 --> I3["regression suite"]
I3 --> R3["Review risk: Test (7 files)"]
R3 --> V3["targeted test run"]
Contributor
OpenCode Review Overview
Pull request overviewOpenCode cannot approve yet because required coverage evidence did not pass. Review outcome1. HIGH .github/workflows/opencode-review.yml:1 - Coverage evidence did not prove required test/docstring evidence
Coverage evidenceCoverage Decision
Changed-File Evidence Mapflowchart LR
Evidence["OpenCode evidence"] --> Review["Current PR review path"]
Review --> Verify["Required checks"]
|
…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>
Contributor
There was a problem hiding this comment.
Pull request overview
OpenCode cannot approve yet because required coverage evidence did not pass.
Review outcome
1. HIGH .github/workflows/opencode-review.yml:1 - Coverage evidence did not prove required test/docstring evidence
-
Problem: The required coverage-evidence job result was
failure, so OpenCode cannot establish approval sufficiency for this head. -
Root cause: Automated approval is only valid when the same-head coverage-evidence job proves supported repository test suites passed and configured docstring gates passed or were advisory, or reports not applicable because no supported source files or package manifests exist. Missing, failed, skipped, unavailable, or unsupported-tooling test evidence is a blocker.
-
Fix: Install or configure the repository test/docstring evidence tooling when source files or package manifests exist, rerun the current-head coverage-evidence job, and approve only after it reports
successwith required evidence or explicit no-source not-applicable evidence. -
Regression test: Keep the approval branch checking
needs.coverage-evidence.result == successbefore posting APPROVE, and publish REQUEST_CHANGES when coverage-evidence blocker states such as cancelled, skipped, failed, unsupported-tooling, or below-100 evidence are present. -
Result: REQUEST_CHANGES
-
Reason: coverage-evidence result was
failure, so required test/docstring evidence was not proven for current head35d343f22785b15e90a729f0bbcd375fc642837a. -
Head SHA:
35d343f22785b15e90a729f0bbcd375fc642837a -
Workflow run: 29311601402
-
Workflow attempt: 1
Coverage evidence
Coverage Decision
- Result: FAIL
- Test evidence: not proven passing
- Docstring evidence: not proven passing when configured
- Failure count: 2
Changed-File Evidence Map
flowchart LR
PR["PR changed files"] --> Evidence["OpenCode bounded evidence"]
Evidence --> S1["Changed file (25 files)"]
S1 --> I1["repository behavior"]
I1 --> R1["Review risk: Changed file (25 files)"]
R1 --> V1["required checks"]
Evidence --> S2["Docs (9 files)"]
S2 --> I2["operator or user guidance"]
I2 --> R2["Review risk: Docs (9 files)"]
R2 --> V2["docs review"]
Evidence --> S3["Test (7 files)"]
S3 --> I3["regression suite"]
I3 --> R3["Review risk: Test (7 files)"]
R3 --> V3["targeted test run"]
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>
Contributor
There was a problem hiding this comment.
Pull request overview
OpenCode cannot approve yet because required coverage evidence did not pass.
Review outcome
1. HIGH .github/workflows/opencode-review.yml:1 - Coverage evidence did not prove required test/docstring evidence
-
Problem: The required coverage-evidence job result was
failure, so OpenCode cannot establish approval sufficiency for this head. -
Root cause: Automated approval is only valid when the same-head coverage-evidence job proves supported repository test suites passed and configured docstring gates passed or were advisory, or reports not applicable because no supported source files or package manifests exist. Missing, failed, skipped, unavailable, or unsupported-tooling test evidence is a blocker.
-
Fix: Install or configure the repository test/docstring evidence tooling when source files or package manifests exist, rerun the current-head coverage-evidence job, and approve only after it reports
successwith required evidence or explicit no-source not-applicable evidence. -
Regression test: Keep the approval branch checking
needs.coverage-evidence.result == successbefore posting APPROVE, and publish REQUEST_CHANGES when coverage-evidence blocker states such as cancelled, skipped, failed, unsupported-tooling, or below-100 evidence are present. -
Result: REQUEST_CHANGES
-
Reason: coverage-evidence result was
failure, so required test/docstring evidence was not proven for current head4e8a143e70c8e21e7cb99c40387ce5479ebe1594. -
Head SHA:
4e8a143e70c8e21e7cb99c40387ce5479ebe1594 -
Workflow run: 29312112905
-
Workflow attempt: 1
Coverage evidence
Coverage Decision
- Result: FAIL
- Test evidence: not proven passing
- Docstring evidence: not proven passing when configured
- Failure count: 2
Changed-File Evidence Map
flowchart LR
PR["PR changed files"] --> Evidence["OpenCode bounded evidence"]
Evidence --> S1["Changed file (25 files)"]
S1 --> I1["repository behavior"]
I1 --> R1["Review risk: Changed file (25 files)"]
R1 --> V1["required checks"]
Evidence --> S2["Docs (9 files)"]
S2 --> I2["operator or user guidance"]
I2 --> R2["Review risk: Docs (9 files)"]
R2 --> V2["docs review"]
Evidence --> S3["Test (7 files)"]
S3 --> I3["regression suite"]
I3 --> R3["Review risk: Test (7 files)"]
R3 --> V3["targeted test run"]
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>
Contributor
There was a problem hiding this comment.
Pull request overview
OpenCode cannot approve yet because required coverage evidence did not pass.
Review outcome
1. HIGH .github/workflows/opencode-review.yml:1 - Coverage evidence did not prove required test/docstring evidence
-
Problem: The required coverage-evidence job result was
failure, so OpenCode cannot establish approval sufficiency for this head. -
Root cause: Automated approval is only valid when the same-head coverage-evidence job proves supported repository test suites passed and configured docstring gates passed or were advisory, or reports not applicable because no supported source files or package manifests exist. Missing, failed, skipped, unavailable, or unsupported-tooling test evidence is a blocker.
-
Fix: Install or configure the repository test/docstring evidence tooling when source files or package manifests exist, rerun the current-head coverage-evidence job, and approve only after it reports
successwith required evidence or explicit no-source not-applicable evidence. -
Regression test: Keep the approval branch checking
needs.coverage-evidence.result == successbefore posting APPROVE, and publish REQUEST_CHANGES when coverage-evidence blocker states such as cancelled, skipped, failed, unsupported-tooling, or below-100 evidence are present. -
Result: REQUEST_CHANGES
-
Reason: coverage-evidence result was
failure, so required test/docstring evidence was not proven for current head7130c3de86611ce91b8ce4734fc9ac70a8a19abe. -
Head SHA:
7130c3de86611ce91b8ce4734fc9ac70a8a19abe -
Workflow run: 29314322214
-
Workflow attempt: 1
Coverage evidence
Coverage Decision
- Result: FAIL
- Test evidence: not proven passing
- Docstring evidence: not proven passing when configured
- Failure count: 2
Changed-File Evidence Map
flowchart LR
PR["PR changed files"] --> Evidence["OpenCode bounded evidence"]
Evidence --> S1["Changed file (27 files)"]
S1 --> I1["repository behavior"]
I1 --> R1["Review risk: Changed file (27 files)"]
R1 --> V1["required checks"]
Evidence --> S2["Docs (10 files)"]
S2 --> I2["operator or user guidance"]
I2 --> R2["Review risk: Docs (10 files)"]
R2 --> V2["docs review"]
Evidence --> S3["Test (7 files)"]
S3 --> I3["regression suite"]
I3 --> R3["Review risk: Test (7 files)"]
R3 --> V3["targeted test run"]
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>
Contributor
There was a problem hiding this comment.
Pull request overview
OpenCode cannot approve yet because required coverage evidence did not pass.
Review outcome
1. HIGH .github/workflows/opencode-review.yml:1 - Coverage evidence did not prove required test/docstring evidence
-
Problem: The required coverage-evidence job result was
failure, so OpenCode cannot establish approval sufficiency for this head. -
Root cause: Automated approval is only valid when the same-head coverage-evidence job proves supported repository test suites passed and configured docstring gates passed or were advisory, or reports not applicable because no supported source files or package manifests exist. Missing, failed, skipped, unavailable, or unsupported-tooling test evidence is a blocker.
-
Fix: Install or configure the repository test/docstring evidence tooling when source files or package manifests exist, rerun the current-head coverage-evidence job, and approve only after it reports
successwith required evidence or explicit no-source not-applicable evidence. -
Regression test: Keep the approval branch checking
needs.coverage-evidence.result == successbefore posting APPROVE, and publish REQUEST_CHANGES when coverage-evidence blocker states such as cancelled, skipped, failed, unsupported-tooling, or below-100 evidence are present. -
Result: REQUEST_CHANGES
-
Reason: coverage-evidence result was
failure, so required test/docstring evidence was not proven for current head77946ddc1af0a3b91dd4612980d485bf3ffe0454. -
Head SHA:
77946ddc1af0a3b91dd4612980d485bf3ffe0454 -
Workflow run: 29316547720
-
Workflow attempt: 1
Coverage evidence
Coverage Decision
- Result: FAIL
- Test evidence: not proven passing
- Docstring evidence: not proven passing when configured
- Failure count: 2
Changed-File Evidence Map
flowchart LR
PR["PR changed files"] --> Evidence["OpenCode bounded evidence"]
Evidence --> S1["Changed file (27 files)"]
S1 --> I1["repository behavior"]
I1 --> R1["Review risk: Changed file (27 files)"]
R1 --> V1["required checks"]
Evidence --> S2["Docs (10 files)"]
S2 --> I2["operator or user guidance"]
I2 --> R2["Review risk: Docs (10 files)"]
R2 --> V2["docs review"]
Evidence --> S3["Test (8 files)"]
S3 --> I3["regression suite"]
I3 --> R3["Review risk: Test (8 files)"]
R3 --> V3["targeted test run"]
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>
Contributor
There was a problem hiding this comment.
Pull request overview
OpenCode cannot approve yet because required coverage evidence did not pass.
Review outcome
1. HIGH .github/workflows/opencode-review.yml:1 - Coverage evidence did not prove required test/docstring evidence
-
Problem: The required coverage-evidence job result was
failure, so OpenCode cannot establish approval sufficiency for this head. -
Root cause: Automated approval is only valid when the same-head coverage-evidence job proves supported repository test suites passed and configured docstring gates passed or were advisory, or reports not applicable because no supported source files or package manifests exist. Missing, failed, skipped, unavailable, or unsupported-tooling test evidence is a blocker.
-
Fix: Install or configure the repository test/docstring evidence tooling when source files or package manifests exist, rerun the current-head coverage-evidence job, and approve only after it reports
successwith required evidence or explicit no-source not-applicable evidence. -
Regression test: Keep the approval branch checking
needs.coverage-evidence.result == successbefore posting APPROVE, and publish REQUEST_CHANGES when coverage-evidence blocker states such as cancelled, skipped, failed, unsupported-tooling, or below-100 evidence are present. -
Result: REQUEST_CHANGES
-
Reason: coverage-evidence result was
failure, so required test/docstring evidence was not proven for current headc2048efe9d6bb600709aabfc8343927325405acf. -
Head SHA:
c2048efe9d6bb600709aabfc8343927325405acf -
Workflow run: 29318123391
-
Workflow attempt: 1
Coverage evidence
Coverage Decision
- Result: FAIL
- Test evidence: not proven passing
- Docstring evidence: not proven passing when configured
- Failure count: 2
Changed-File Evidence Map
flowchart LR
PR["PR changed files"] --> Evidence["OpenCode bounded evidence"]
Evidence --> S1["Changed file (27 files)"]
S1 --> I1["repository behavior"]
I1 --> R1["Review risk: Changed file (27 files)"]
R1 --> V1["required checks"]
Evidence --> S2["Docs (10 files)"]
S2 --> I2["operator or user guidance"]
I2 --> R2["Review risk: Docs (10 files)"]
R2 --> V2["docs review"]
Evidence --> S3["Test (8 files)"]
S3 --> I3["regression suite"]
I3 --> R3["Review risk: Test (8 files)"]
R3 --> V3["targeted test run"]
…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>
Contributor
There was a problem hiding this comment.
Pull request overview
OpenCode cannot approve yet because required coverage evidence did not pass.
Review outcome
1. HIGH .github/workflows/opencode-review.yml:1 - Coverage evidence did not prove required test/docstring evidence
-
Problem: The required coverage-evidence job result was
failure, so OpenCode cannot establish approval sufficiency for this head. -
Root cause: Automated approval is only valid when the same-head coverage-evidence job proves supported repository test suites passed and configured docstring gates passed or were advisory, or reports not applicable because no supported source files or package manifests exist. Missing, failed, skipped, unavailable, or unsupported-tooling test evidence is a blocker.
-
Fix: Install or configure the repository test/docstring evidence tooling when source files or package manifests exist, rerun the current-head coverage-evidence job, and approve only after it reports
successwith required evidence or explicit no-source not-applicable evidence. -
Regression test: Keep the approval branch checking
needs.coverage-evidence.result == successbefore posting APPROVE, and publish REQUEST_CHANGES when coverage-evidence blocker states such as cancelled, skipped, failed, unsupported-tooling, or below-100 evidence are present. -
Result: REQUEST_CHANGES
-
Reason: coverage-evidence result was
failure, so required test/docstring evidence was not proven for current headb5d9d90c0a397c5bbeb4aa864f4fc6b291cd8936. -
Head SHA:
b5d9d90c0a397c5bbeb4aa864f4fc6b291cd8936 -
Workflow run: 29318985676
-
Workflow attempt: 1
Coverage evidence
Coverage Decision
- Result: FAIL
- Test evidence: not proven passing
- Docstring evidence: not proven passing when configured
- Failure count: 2
Changed-File Evidence Map
flowchart LR
PR["PR changed files"] --> Evidence["OpenCode bounded evidence"]
Evidence --> S1["Changed file (27 files)"]
S1 --> I1["repository behavior"]
I1 --> R1["Review risk: Changed file (27 files)"]
R1 --> V1["required checks"]
Evidence --> S2["Docs (10 files)"]
S2 --> I2["operator or user guidance"]
I2 --> R2["Review risk: Docs (10 files)"]
R2 --> V2["docs review"]
Evidence --> S3["Test (8 files)"]
S3 --> I3["regression suite"]
I3 --> R3["Review risk: Test (8 files)"]
R3 --> V3["targeted test run"]
…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>
Problem The Rust EAP GPU kernel was implemented but the core convenience API, PyO3 binding, and Python serving API all defaulted to CPU. An explicit GPU request also fell back silently, making it possible to mistake CPU fallback for an actual GPU run. Reproduction/Evidence CodeGraph traced score_respondents -> score_bank_eap -> score_eap_device -> try_score_eap_gpu and showed the public defaults were CPU. WGPU_BACKEND=metal selected an Apple Metal adapter; CPU versus GPU maximum absolute differences were loglik 5.391e-7, theta 4.912e-7, theta_sd 7.951e-7, and xi 2.330e-7 against a 2e-3 tolerance. Root cause GPU scoring was introduced as opt-in and later public layers retained the original CPU defaults. The GPU parity test treated a missing adapter as a soft skip even when Metal had been explicitly requested. Change Default EAP execution to Device::Auto in Rust and device=auto in PyO3 and serving, preserve explicit CPU f64 execution, warn when explicit GPU requests fall back, and require an actual adapter when WGPU_BACKEND=metal is set. Add public-default and fallback regressions plus changelog documentation. Validation WGPU_BACKEND=metal cargo test -p mlsirm-core --release scoring:: -- --nocapture: 24 passed, 0 failed, 1 unrelated ignored; max GPU deltas below 8e-7. cargo test -p mlsirm-core --no-default-features default_eap_policy_matches_auto_device: 1 passed. cargo test --release --manifest-path crates/fast-mlsirm-py/Cargo.toml: 3 passed. WGPU_BACKEND=metal uv run pytest tests/test_serving.py -ra: 6 passed. uv run pytest --collect-only -q: 697 collected. cargo test --workspace --release -- --list: 436 tests; source audit found 39 pre-existing ignores. Focused Ruff and git diff --check passed. Sources Bock, R. D., & Mislevy, R. J. (1982). Adaptive EAP estimation of ability in a microcomputer environment. Applied Psychological Measurement, 6(4), 431-444. https://doi.org/10.1177/014662168200600405 Existing Zotero record and attached PDF metadata were verified; no Zotero item was changed.
Problem score_respondents accepted a device policy, but when the compiled extension was unavailable its EAP path silently ran the NumPy reference implementation. Even device=gpu therefore succeeded on CPU outside Rust, contradicting the GPU-first contract and hiding an unavailable accelerator. Reproduction/Evidence Before this change, monkeypatching fast_mlsirm.serving._core_module to return None and calling score_respondents on a two-item MIRT bundle with device=gpu returned theta=[2.0816681711721685e-17], theta_sd=[0.83473056368705], xi=[0.0], and loglik=-1.5767249033611719 instead of rejecting the ignored device request. Root cause The serving EAP branch used score_bank_eap only when the extension loaded, then imported and invoked estimators.marginal.score_eap as a transparent fallback. That fallback had no device parameter and executed NumPy CPU code. Change Require the compiled Rust core for serving EAP, remove the transparent NumPy fallback/import, and add a regression test proving an explicit GPU request cannot escape to Python CPU. Rust score_eap_device retains GPU-preferred auto execution and Rust CPU f64 fallback. Validation - uv run ruff check python/fast_mlsirm/serving.py: passed - uv run ruff check --ignore E741,F841 tests/test_serving.py: passed - uv run pytest tests/test_serving.py -ra: 7 passed in 47.70s - uv run pytest --collect-only -q: 698 tests collected - CARGO_TARGET_DIR=/tmp/fast-mlsirm-pr160-gpu-7644999 cargo test --workspace --release -- --list: 436 tests, 0 benchmarks - WGPU_BACKEND=metal CARGO_TARGET_DIR=/tmp/fast-mlsirm-pr160-gpu-7644999 cargo test -p mlsirm-core --release gpu_eap_matches_cpu_reduction -- --nocapture: 1 passed; max absolute differences loglik=5.391e-7, theta=4.912e-7, theta_sd=7.951e-7, xi=2.330e-7, tolerance=2e-3 - git diff --check: passed Sources Bock, R. D., & Mislevy, R. J. (1982). Adaptive EAP estimation of ability in a microcomputer environment. Applied Psychological Measurement, 6(4), 431-444. https://doi.org/10.1177/014662168200600405 This correction changes execution policy only; it does not alter the cited EAP mathematics.
Problem: CAT EAP scoring already preferred the Rust wgpu path, but fixed-bank item and test information remained Rust CPU-only. The complete adaptive selection step therefore did not satisfy the GPU-first execution contract, and a first f32 kernel could leak invalid information for values that remain representable in the f64 public contract. Reproduction/Evidence: On actual Metal, ordinary fixed-bank information matched CPU f64 with maximum absolute differences 1.187e-7 for item information and 1.954e-7 for test information. A valid alpha=100 public serving case returned CPU [[0.0]] but the direct f32 GPU path produced NaN before the readback guard. Root cause: bank_information had no device dispatch or GPU kernel. WGSL f32 exp(alpha) also has a smaller finite range than the Rust f64 reference, and the initial readback did not reject non-finite or negative information. Change: Add cached wgpu information kernels for MIRT, distance LSIRM, and inner-product bifactor banks. Make bank information and CAT default to Device::Auto, expose the device through PyO3 and serving, retain explicit Rust CPU fallback, reject invalid GPU readbacks, and harden shape/finite-input validation. Validation: - WGPU_BACKEND=metal cargo test -p mlsirm-core --release gpu_bank_information_matches_cpu_reduction -- --nocapture: 1 passed; actual Metal; item/test max abs 1.187e-7/1.954e-7 - cargo test -p mlsirm-core --release scoring::cat_pv_tests -- --nocapture: 3 passed - cargo test overflow regressions with default and no-default features: passed - uv run pytest tests/test_serving.py -q -ra: 8 passed in 24.07s; 0 skipped/xfail/xpass/deselected - uv run pytest --collect-only -q: 699 collected - focused Ruff, rustfmt, PyO3 checks/tests, release extension build, and git diff --check: passed Sources: No statistical formula or source claim changed. Existing Bock and Mislevy (1982) EAP and Magis (2013) information/CAT citation scope is retained; GPU dispatch and f32 fallback are repository implementation policy.
Problem Plausible-value posterior reduction and seeded sampling remained a scalar Rust CPU-only path. The public PyO3 and serving APIs exposed no device policy, and the CPU implementation did not bound parallel work by hardware concurrency. Reproduction/Evidence CodeGraph traced serving.plausible_values through the PyO3 plausible_values binding to scoring::plausible_values, where the posterior grid reduction and categorical draws ran only in a single CPU loop. An explicit Metal parity run now exercises the spatial MLS2PLM xi and theta selections over 2,048 outputs. Root cause The original sampler coupled its LCG stream to the scalar loop and had no GPU dispatcher or reusable posterior buffers, so it could neither preserve deterministic draws across devices nor fall back from bounded GPU failures. Change Add two-pass wgpu kernels for posterior log-sum-exp and seeded xi/theta categorical selection. Default Rust, PyO3, and Python APIs to Device::Auto; validate f32 outputs and adapter/buffer limits before using them. Retain an explicit deterministic f64 CPU path using fixed contiguous person shards capped by available_parallelism to minimize scheduling and context switching. Validation WGPU_BACKEND=metal cargo test -p mlsirm-core plausible_values --all-features -- --nocapture: 2 passed, max abs CPU/GPU difference 7.652738309715801e-8, mean difference 2.44e-9. cargo test --workspace: 400 passed, 39 ignored, 0 failed. cargo check --workspace --all-features and cargo check -p mlsirm-core --no-default-features passed. uv run --no-sync pytest -q -ra: 700 passed with no skips/xfails/xpasses. Final rebuilt-extension selection: 9 passed, 259 deselected; public Metal max abs difference 7.652738309716e-8. Focused rustfmt, Ruff, and git diff checks passed; whole cargo fmt remains blocked only by pre-existing tests/unit/dif_tests.rs drift. Sources Marsman, M., Maris, G., Bechger, T., & Glas, C. (2016). What can we learn from plausible values? Psychometrika, 81(2), 274-289. https://doi.org/10.1007/s11336-016-9497-x. The GPU and CPU scheduling policy is repository-specific and is not attributed to this source.
Problem EAPsum conversion-table construction ran only in Rust CPU, while serving summed responses and looked up tables in Python. This bypassed the repository Rust GPU-first execution contract and left the fallback outside the required Rust-only numerical path. Reproduction/Evidence CodeGraph traced score_respondents(method=eapsum) to a Python score loop and eapsum_tables to a scalar Rust CPU reduction with no Device dispatch. Explicit WGPU_BACKEND=metal parity coverage was absent. Root cause The original EAPsum feature predated the shared auto|gpu|cpu dispatch policy. Its Lord-Wingersky recursion and posterior moments were implemented as one CPU-only helper, and PyO3 exposed tables but no Rust lookup endpoint. Change Add bounded wgpu kernels for Lord-Wingersky table construction, posterior moments, and respondent lookup. Add Device-aware Rust table and lookup APIs, PyO3 bindings, and Python delegation. Retain f64 CPU fallbacks using fixed contiguous available_parallelism shards to limit scheduling and context switching. Validate response completeness, table shapes, finite values, GPU bounds, and f32 outputs. Document the verified APA 7 source. Validation - WGPU_BACKEND=metal cargo test -p mlsirm-core eapsum -- --nocapture: 3 passed; actual Metal max errors score_prob 1.622e-7, EAP 2.706e-6, SD 3.742e-6, lookup EAP 5.068e-9, lookup SD 2.481e-8. - WGPU_BACKEND=metal cargo test --workspace -- --nocapture: 402 passed, 39 ignored, 0 failed. - WGPU_BACKEND=metal .venv/bin/pytest -q -ra: 701 passed. - Focused Python serving/scoring/security selection: 24 passed, 268 deselected. - cargo check --workspace --all-features; cargo check -p mlsirm-core --no-default-features; PyO3 check/build: passed. - Focused rustfmt, Ruff, and git diff --check: passed. - Strict Clippy remains blocked by 148 pre-existing repository-wide warnings; no final diagnostic points at gpu_eapsum.rs. Sources Thissen, D., Pommerich, M., Billeaud, K., & Williams, V. S. L. (1995). Item response theory for scores on tests including polytomous items with ordered responses. Applied Psychological Measurement, 19(1), 39-49. https://doi.org/10.1177/014662169501900105
The library already had the full polytomous model family and polytomous EAP scoring, but its only bias-reduced ML ability estimator was dichotomous-only. EAP shrinks toward the population mean, which is exactly what individual score reporting must not do, so this closes a gap rather than adding a third way to do the same thing. Solves dlnL/dtheta + J/(2I) = 0 with I = sum_k P'_k^2/P_k and J = sum_k P'_k P''_k/P_k over the person's observed items -- the exact generalization of the shipped dichotomous sum_i P'P''/(PQ), which is its two-category case. J is computed DIRECTLY, never as a derivative of I. GRM and GPCM; PCM is the GPCM path at slope = 1. RSM is deliberately not supported and the code says why: its fitted (delta, shared tau) form is not convertible through any exposed API. Per-category quantities are formed division-free from the sigmoids, so no category probability ever appears in a denominator and no probability floor is needed. Both polytomous log-likelihoods are log-concave, yet the weighted objective is genuinely multimodal because the Warm weight is not: a 3-item GPCM bank in the suite has stationary points at +0.0988, +0.3774 and +1.3314 while max lnL'' < 0, so a solver bracketing the first sign change errs by 1.23 logits (2.36 on the GRM fixture). The global grid scan of score_wle is therefore reused unchanged, with the grid demand scaled by (n_cat - 1). VERIFICATION STATUS, stated in the code because it bounds what may be claimed. That the polytomous Warm correction is J/(2I) with J = sum_k P'P''/P is confirmed from the catR package's SOURCE, not from a primary paper -- and catR keeps its Jeffreys-prior branch as a separate expression, so the two estimators are kept distinct here too. Penfield and Bergeron (2005) treat the GPCM but their equations were not obtainable, so nothing here rests on them and they are not cited as a source. Separately, and proved in-repository: J = I' holds exactly for both shipped families, because J - I' = -E[l' l''] vanishes (the GPCM's l'' is category-free; the GRM telescopes through v_0 = v_K = 0). That identity is used ONLY as a test oracle, never as a shortcut, since it fails for per-boundary slopes and for the 3PL -- both shipped as negative controls. One coverage limit is stated rather than papered over: because J = I' is exact here, replacing J with a numerical derivative of I is behaviour-preserving for these two families and NO polytomous test can detect it. The discriminating anchors live in the dichotomous suite, where a lower asymptote breaks the identity. A family added later must re-derive J. Also corrects the dichotomous WLE documentation, which claimed J coincides with I'/2 for the 2PL/Rasch. The correct statement is J = I' exactly there, from I' = 2J - T with T = sum_i P'^3 (1 - 2P)/(PQ)^2, and T = J only when c = 0, d = 1. The identity is now pinned by a test, because the first attempt at this correction was itself wrong -- it dropped the (1 - 2P) factor, giving a value ~5x off at the 2PL -- and nothing caught it. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Problem: Empirical EAP reliability was a CPU-only scalar reduction with no execution-device contract, contrary to the GPU-first Rust policy. Python callers could not request a CPU reference or verify an actual accelerator path. Reproduction/Evidence: CodeGraph traced python/fast_mlsirm/fitstats.py through the PyO3 empirical_reliability wrapper to scoring.rs, where every dimension and person was reduced serially. On the fixed 513-person fixture, CPU f64 returns [0.9153743829939073, 0.7196718700963158] and Metal f32 returns [0.9153743982315063, 0.7196718454360962], a maximum absolute difference of 2.4660219644090375e-08. Root cause: The reliability API predated the repository Device dispatch convention, and gpu_scoring.rs only exposed item/test information kernels. Change: Add an actual wgpu empirical-reliability kernel with finite f32 conversion, buffer/workgroup bounds, and validated readback. Default Rust and Python calls now use Device::Auto, explicit cpu/gpu/auto remains available, and unavailable/invalid GPU work falls back to a fixed contiguous-shard f64 CPU reduction sized by available_parallelism to avoid oversubscription. Validation: - WGPU_BACKEND=metal cargo test -p mlsirm-core empirical_reliability -- --nocapture: 3 passed, including a strict actual-adapter test - cargo test -p mlsirm-core --no-default-features empirical_reliability -- --nocapture: 2 passed - cargo check --workspace --all-targets: passed - WGPU_BACKEND=metal uv run pytest tests/test_fitstats.py tests/test_paper_features.py -k "empirical_reliability or wle" -ra: 6 passed, 135 deselected - WGPU_BACKEND=metal cargo test -p mlsirm-core wle -- --nocapture: 19 passed, 1 existing literature-grade test ignored - Pre-arrival parent 9b4b297 full suites: Python 702 passed; Rust 404 passed, 39 ignored - Ruff, focused rustfmt, CodeGraph re-index, and git diff --check: passed Sources: Bechger, T. M., Maris, G., Verstralen, H. H. F. M., & Beguin, A. A. (2003). Using classical test theory in combination with item response theory. Applied Psychological Measurement, 27(5), 319-334. https://doi.org/10.1177/0146621603257518 Stanley, L. M., & Edwards, M. C. (2016). Reliability and model fit. Educational and Psychological Measurement, 76(6), 976-985. https://doi.org/10.1177/0013164416638900
Problem: The current-head Rust CI ran the actual-GPU empirical-reliability parity test on an Ubuntu runner with no usable adapter. The test unconditionally unwrapped the optional GPU result, turning an environmental constraint into a product failure. Reproduction/Evidence: GitHub Actions run 29929333811, job 88954572052 failed at scoring::reliability_tests::empirical_reliability_explicit_gpu_matches_cpu_reference with 414 passed, 1 failed and 39 ignored. The panic was "this test requires a usable GPU adapter" at tests/unit/scoring_reliability_tests.rs:70. Root cause: Unlike the existing EAP, EAPsum and bank-information GPU parity tests, the reliability regression did not distinguish an unavailable adapter from a requested Metal backend. Change: Treat an absent adapter as an explicitly reported environment skip in ordinary runs. Keep WGPU_BACKEND=metal strict: an explicit Metal request still fails if no Metal adapter is selected, and available adapters still execute the CPU-f64 versus GPU-f32 parity assertions. Validation: - WGPU_BACKEND=metal cargo test -p mlsirm-core scoring::reliability_tests::empirical_reliability_explicit_gpu_matches_cpu_reference -- --nocapture: 1 passed - WGPU_BACKEND=metal cargo test --workspace -- --nocapture: 415 passed, 0 failed, 39 ignored; QMC correlated D4 converged at 110/200 with 9.492470326222247e-6 < 1e-5 - uv run pytest tests/test_fitstats.py -k empirical_reliability -ra: 3 passed, 22 deselected - pytest collect-only: 703 tests; Rust list: 454 tests - rustfmt --edition 2021 --check tests/unit/scoring_reliability_tests.rs and git diff --check: passed Sources: This changes only environment-sensitive test dispatch, not the empirical-reliability formula or documentation. Existing verified sources remain Bechger et al. (2003), https://doi.org/10.1177/0146621603257518, and Stanley and Edwards (2016), https://doi.org/10.1177/0013164416638900.
Contributor
There was a problem hiding this comment.
Pull request overview
OpenCode cannot approve yet because required coverage evidence did not pass.
Review outcome
1. HIGH .github/workflows/opencode-review.yml:1 - Coverage evidence did not prove required test/docstring evidence
-
Problem: The required coverage-evidence job result was
failure, so OpenCode cannot establish approval sufficiency for this head. -
Root cause: Automated approval is only valid when the same-head coverage-evidence job proves supported repository test suites passed and configured docstring gates passed or were advisory, or reports not applicable because no supported source files or package manifests exist. Missing, failed, skipped, unavailable, or unsupported-tooling test evidence is a blocker.
-
Fix: Install or configure the repository test/docstring evidence tooling when source files or package manifests exist, rerun the current-head coverage-evidence job, and approve only after it reports
successwith required evidence or explicit no-source not-applicable evidence. -
Regression test: Keep the approval branch checking
needs.coverage-evidence.result == successbefore posting APPROVE, and publish REQUEST_CHANGES when coverage-evidence blocker states such as cancelled, skipped, failed, unsupported-tooling, or below-100 evidence are present. -
Result: REQUEST_CHANGES
-
Reason: coverage-evidence result was
failure, so required test/docstring evidence was not proven for current head1243031dc5f4c85b9a965270d6f131cb309730f2. -
Head SHA:
1243031dc5f4c85b9a965270d6f131cb309730f2 -
Workflow run: 29937067070
-
Workflow attempt: 1
Coverage evidence
Coverage Decision
- Result: FAIL
- Test evidence: not proven passing
- Docstring evidence: not proven passing when configured
- Failure count: 2
Changed-File Evidence Map
flowchart LR
Evidence["OpenCode evidence"] --> Review["Current PR review path"]
Review --> Verify["Required checks"]
Contributor
Author
|
Superseded by #213 for the current validation-and-hardening stream. |
seonghobae
added a commit
that referenced
this pull request
Jul 25, 2026
* feat(marginal): MMLE-EM estimator for the latent-space family with multigroup/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>
* feat(fitstats): S-X2, l_z/l_z* person fit, item screening pipeline, and 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>
* docs(changelog) + feat(cli): document the marginal-estimator feature 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>
* feat: QMC/MC-EM rules, Rust scoring (EAP/MAP/EAPsum) + fitstats core, 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>
* docs: phase-2 design notes (QMC/MC-EM, scoring, FIPC, numerical caveats)
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* feat: paper batch 1 — zero-inflated mixture, position covariate, validation 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>
* feat(oakes): Oakes-identity standard errors for marginal fits (Pritikin 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>
* docs(papers): implemented-literature map + group A/B/C implementation specs
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* feat: paper batch 3 — information/CAT/plausible values, residual & pairwise 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>
* feat(fitstats): Chen-Thissen (1997) local-dependence indices; batch-4 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>
* feat(scoring): empirical EAP reliability (Stanley-Edwards 2016; Milanzi 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>
* feat(bifactor): BIFAC2PLM via an inner-product interaction kind
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>
* feat(fitstats): M2 limited-information GOF with RMSEA2 CI and SRMSR
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>
* feat(linking): moment + Haebara + Stocking-Lord IRT scale linking
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>
* test(coverage): cover linking + M2 branches; docstrings on new helpers
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>
* fix(security): harden untrusted input boundaries (Strix VULN-0001..0007)
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>
* test(coverage): cover node/quadrature/parse/validate defensive branches
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>
* fix(security): finite/bounded numeric config + n_draws/n_dims guards
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>
* fix(config): keep validation messages backward-compatible after hardening
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>
* feat(gpu): GPU EAP scoring kernel (score_pass), opt-in via Device
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>
* security: harden 11 second-pass Strix findings on PR #160 paper-feature 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>
* security: chase Strix batch-3 (11 findings) + proactive boundary audit (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>
* security: Rust-core backstops from the boundary audit (OOB/panic guards)
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>
* docs: GPCM/nominal polytomous kernel design spec (for the next model-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>
* gpcm: NumPy parity oracle for the unified softmax category cell
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>
* gpcm: analytic M-step gradient reference (FD-verified)
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>
* gpcm: NumPy GPCM MMLE-EM reference fitter (recovery-verified)
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>
* docs: resolve the polytomous-LSIRM design fork from the literature
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>
* gpcm: polytomous response cells + gradients in Rust (the compute path)
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>
* docs: keep the polytomous design spec model-generic
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>
* gpcm: unidimensional polytomous marginal-EM fitter in Rust
mlsirm-core::poly::fit_poly_unidim — the Rust compute path for polytomous
fitting (no latent-space interaction yet): Bock-Aitkin marginal EM over the
crate's Gauss-Hermite grid (bit-identical to the NumPy hermegauss reference),
reusing the poly cells and gradients, with a per-item Newton M-step
(finite-difference Hessian on the analytic gradient, small dense solve).
Supports both PolyModel::Grm (default) and ::Gpcm.
Rust recovery tests (deterministic LCG simulation) for both models: slope
correlation > 0.9 and GRM threshold MAE < 0.25. cargo -p mlsirm-core --lib
75 pass. Next: fold the poly cell into the latent-space marginal E-step
(build_tables/person_pass) and expose via PyO3 with cross-language parity.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: expose the Rust polytomous fitter via PyO3 with cross-language parity
PyO3 wrappers for the polytomous compute path: _core.gpcm_cell_logprobs,
_core.grm_cell_logprobs, and _core.fit_poly_unidim(y, ..., model="grm"|"gpcm").
The Python-side category_logprobs/gpcm_node_gradient/fit_gpcm_numpy stay as
parity references only.
Parity test: the Rust GPCM cell matches the NumPy category_logprobs reference
to <=1e-12, and the Rust unidimensional fitter agrees with fit_gpcm_numpy on
recovered slopes (same EM/Newton, bit-identical GH grid). Full Python suite
352 pass; cargo -p mlsirm-core --lib 75 pass.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: public fit_polytomous API (GRM/GPCM, compute in Rust)
fit_polytomous(responses, n_cat, model="grm"|"gpcm") — a thin orchestration
wrapper over mlsirm_core::poly::fit_poly_unidim; validates inputs and returns
slopes, category parameters, log-likelihood, and (GPCM) Muraki step thresholds.
GRM is the default. Exported from the package top level.
Test: GRM recovery through the public API (slope corr > 0.9) plus input-
validation rejections. Full suite 353 pass. The latent-space polytomous
extension (same cell inside the marginal quadrature) remains the next milestone.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: polytomous EAP scoring (fit -> score), compute in Rust
mlsirm_core::poly::score_poly_eap computes EAP trait scores + posterior SDs
from polytomous responses given fitted item parameters, over the Gauss-Hermite
grid; exposed via PyO3 and wrapped by the public score_polytomous(responses,
fit). Completes the fit -> score pair for GRM/GPCM.
Tests: Rust score_poly_eap recovers true theta (corr > 0.8 at the true item
params); Python fit_polytomous -> score_polytomous round-trip recovers theta.
cargo -p mlsirm-core --lib 76 pass; full pytest 354 pass.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: NumPy GRM cell parity mirror (completes the cell parity contract)
grm_category_logprobs is the NumPy parity reference for the Rust GRM
cumulative-logit cell (mlsirm_core::poly::grm_logprobs); previously only the
GPCM softmax cell had a cross-language mirror. Parity test: Rust
grm_cell_logprobs matches the NumPy reference to 1e-12, plus NumPy
self-consistency (normalization + binary reduction).
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: polytomous item/test information curves (compute in Rust)
mlsirm_core::poly::poly_item_information + poly_information_curves compute
Fisher information I(theta) = sum_k (dP_k/dtheta)^2 / P_k for GRM and GPCM
(GPCM reduces to a^2*Var_P(scores); GRM to the cumulative-slope form), exposed
via PyO3 and wrapped by information_polytomous(fit, theta) -> item/test info.
Tests: analytic information matches a central finite-difference of the cell to
1e-4 (Rust); the Python API returns positive curves with test = item-info row
sum. cargo -p mlsirm-core --lib 77 pass; full pytest 356 pass.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* docs: record polytomous scoring and information in the changelog
The polytomous entry now describes the full fit -> score -> information
subsystem (fit_polytomous / score_polytomous / information_polytomous), all
compute in the Rust core.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: missing-data support for the polytomous fitter and scorer
fit_poly_unidim and score_poly_eap take an optional observed mask; missing
responses are marginalized out of each person's likelihood/posterior (the same
treatment the binary marginal path already applies). Exposed through PyO3 and
the Python API, where NaN in responses marks a missing cell (a shared
_poly_int_and_mask validator).
Tests: Rust GPCM recovery under ~25% MCAR missingness (slope corr > 0.9);
Python fit_polytomous/score_polytomous round-trip on NaN-holed data.
cargo -p mlsirm-core --lib 78 pass; full pytest 357 pass.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* docs: note missing-data (NaN) handling in the polytomous changelog entry
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: latent-space polytomous LSIRM estimator + rigorous validation
mlsirm_core::poly_marginal::fit_poly_lsirm — the polytomous (GRM/GPCM) cell
embedded in an interaction map, fit by marginal EM. Unidimensional trait theta;
the person latent position xi is integrated over a tensor Gauss-Hermite grid and
the item position zeta_i is estimated. Fixed gamma=1 (Go et al. 2024 lsirm12pl
identification). Fully additive: reuses the poly cells/gradients, quadrature,
and the exact binary chain rule (per-node poly g_base multiplies d base/d a =
a*theta and d base/d zeta_k = (xi_k - zeta_k)/dist) — touches neither the binary
estimator nor the GPU. Design confirmed by a marginal.rs mapping workflow.
Validation strengthened per review — absolute agreement, not correlation:
- Recovery test uses RMSE: item-slope RMSE and item-item distance-matrix RMSE
(the distance matrix is exactly invariant to the position rotation/
reflection/translation ambiguity, and gamma=1 fixes its scale).
- Cross-validation against an already-validated reference (not self-recovery):
at K=2 the GPCM cell IS the 2PL, so fit_poly_unidim(K=2) must reproduce the
repo's binary MMLE-EM (mmle::fit_mmle_2pl) item parameters on the same data —
slope and intercept RMSE < 0.1.
cargo -p mlsirm-core --lib 80 pass.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: extend K=2 cross-validation to GRM + RMSE under missingness
The K=2-vs-trusted-binary-MMLE cross-check now runs for both cells (GRM, the
default, also reduces to the 2PL at K=2), so both match mmle::fit_mmle_2pl to
RMSE < 0.1. The missing-data recovery test gains an absolute slope RMSE
assertion alongside the correlation check.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: person scores + Python API for the latent-space polytomous LSIRM
fit_poly_lsirm now returns per-person EAP trait scores (theta_eap/theta_sd) and
EAP latent positions (xi_eap) from a post-convergence posterior pass over the
(theta, xi) grid. Exposed via PyO3 and the public fit_lsirm_polytomous(
responses, n_cat, latent_dim, model=...) -> PolyLsirmFit.
Tests: the Rust recovery test adds person trait-EAP correlation (EAP is
prior-shrunk, so correlation is the right metric there); a Python end-to-end
test recovers item positions (distance-matrix RMSE) and slopes (RMSE) and
checks the returned person scores. cargo -p mlsirm-core --lib 80 pass; full
pytest 358 pass.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: polytomous model-selection indices (Kang, Cohen & Sung 2009)
polytomous_information_criteria(fit, n_persons) returns AIC/BIC/CAIC/AICc/SABIC
and the free-parameter count for a fitted GRM/GPCM or latent-space polytomous
model, from its marginal log-likelihood and parameter count (Kang, Cohen &
Sung 2009, Model Selection Indices for Polytomous Items). Exported at the
package top level.
Test: correct parameter count (n_items*n_cat for the unidimensional cell),
finite indices, and the expected AIC < BIC < CAIC ordering at N=1500.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: generalized S-X2 polytomous item fit (Kang & Chen 2008/2011)
Extend the binary Orlando-Thissen S-X2 to ordered polytomous items with
the generalized Lord-Wingersky summed-score recursion (Thissen, Pommerich,
Billeaud & Williams 1995): group persons by summed score, form the model-
expected category proportions
E_ikz = INT P_i(z|t) f*_i(k-z|t) phi(t) dt / INT f(k|t) phi(t) dt
from the leave-one-out distribution f*_i, merge boundary score groups and
collapse adjacent categories to a minimum expected frequency, and report
per-item chi-square with df = sum_g (cells_g - 1) - m.
All numerical work in Rust (poly.rs), reusing the GRM/GPCM cells.
Validation:
- reduces EXACTLY (< 1e-8, df equal) to the trusted binary fitstats::s_x2
at n_cat=2 for both GRM and GPCM on a shared quadrature grid;
- at the true generating parameters E[S-X2] tracks the retained cell count
(ratio in [0.85, 1.15]) with a sub-15% flag rate, for GPCM (2008) and
GRM (2011) alike -- the calibration a mis-scaled index (Yen Q1 / G2)
fails.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: Python API for generalized S-X2 polytomous item fit + APA 7th refs
Expose poly_s_x2 through PyO3 (poly_item_fit_sx2) and the public
item_fit_polytomous(responses, fit) wrapper: per-item statistic, df,
p_value, n_cells for a fitted GRM/GPCM, with NaN = missing marginalized to
a complete-case summed-score table. Python test covers calibration at the
fitted model (statistic ~ df, low flag rate), the df = n_cells - m identity,
missing data, and input validation.
Cite the methodological basis in APA 7th ed. form in the Rust doc comments
and Python docstrings (Kang & Chen 2008/2011; Orlando & Thissen 2000;
Thissen, Pommerich, Billeaud & Williams 1995; Lord & Wingersky 1984), per
the standing documentation directive.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: Monte-Carlo recovery study (normal vs skew ability), RMSE + bias
Add a parameter-recovery Monte-Carlo for the GPCM fitter generating from the
published Kang & Chen (2008, p. 397) item scheme (slope ~ lognormal(0, .5^2),
four step difficulties ~ N(+-1.5, +-0.5; SD .5)), across two ability
conditions:
- normal theta ~ N(0,1) (matched prior): slope RMSE .08 / |bias| .01,
intercept RMSE .17 / |bias| .04 -- tight, near-unbiased recovery;
- skew theta = Exp(1)-1 (mean 0, var 1, skew 2), a prior
misspecification: slope RMSE .15 / |bias| .11, intercept RMSE .21 /
|bias| .10 -- recovery holds but degrades with measurable bias.
Reports absolute-agreement RMSE and signed bias per parameter over
replications (not correlation), printed under --nocapture. APA 7th refs in
the doc comment.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: split MC recovery into fast CI guard + 500-replication study
Refactor the GPCM recovery Monte-Carlo into a reusable helper with two
entry points: a fast regression guard (20 reps, run in CI) and an
#[ignore]-d literature-grade study (500 replications, N=2000, run with
`cargo test --release -- --ignored --nocapture`), matching the >=500-rep
norm of the IRT Monte-Carlo literature.
500-replication results (Kang & Chen 2008 generating scheme):
- normal theta ~ N(0,1): slope RMSE .072 / |bias| .005, intercept
RMSE .137 / |bias| .010 -- near-unbiased under the matched prior;
- skew theta = Exp(1)-1: slope RMSE .146 / |bias| .107, intercept
RMSE .185 / |bias| .096 -- systematic bias from prior misspecification
(stable across 500 reps, not Monte-Carlo noise).
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: polytomous M2 limited-information GOF (Maydeu-Olivares & Joe 2014)
Add fitstats::poly_m2 — the ordered-category generalization of the binary
m2_rmsea2. Uses the cumulative marginal form (univariate P(Y_i>=c),
bivariate P(Y_i>=c,Y_j>=d)), provably the same statistic as the paper's
category-equality form and reducing EXACTLY to m2_rmsea2 at K=2. Model
moments factor over the N(0,1) grid by local independence; the multinomial
covariance Xi uses the nesting rule 1{Y>=c}1{Y>=c'}=1{Y>=max} (max-threshold
collapse); reuses the one-Cholesky residual-projection solve, RMSEA2 + CI,
and SRMSR. df = Q - P, Q = n(K-1) + C(n,2)(K-1)^2, P = n*K.
Spec derived and adversarially verified by a subagent workflow against the
paper (cumulative==equality confirmed numerically to 7e-11) and the existing
binary code; guards for lengths, category range, and degenerate margins added
per the verification.
Validation:
- K=2 reduces EXACTLY to the trusted binary m2_rmsea2 (GRM and GPCM: M2, df,
p-value, RMSEA2, moment/param counts) to <1e-4;
- Monte-Carlo calibration (GPCM, 500 reps, N=2000): under a matched N(0,1)
ability mean(M2)/df=0.994 with Type I 0.050 (nominal); under a right-skewed
population (Exp(1)-1) mean(M2)/df=4.15 with power 1.00 -- exactly the null
calibration + misfit power M2 is meant to show. Fast CI guard + #[ignore]
500-rep study. APA 7th references in the doc comment.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: Python API for polytomous M2 goodness-of-fit + APA 7th refs
Expose poly_m2 through PyO3 and the public m2_polytomous(responses, fit)
wrapper: m2, df, p_value, rmsea2 (+90% CI), srmsr and the moment/parameter/
complete-case counts for a fitted GRM/GPCM, with NaN = missing marginalized to
complete cases. Python test verifies the df = n(K-1)+C(n,2)(K-1)^2 - nK
bookkeeping, a good fit for correctly-specified data, rejection under a
misspecified fit, and the too-few-items guard. APA 7th references in the Rust
doc comment and Python docstring.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: polytomous item-pair local dependence (Chen & Thissen 1997)
Add fitstats::poly_local_dependence -- per-pair LD diagnostics for a fitted
GRM/GPCM, the ordered-category generalization of the binary pairwise
chi-square. For each item pair it compares the observed K x K contingency
table to the model-implied joint under local independence
(E_ab = N * sum_t w_t P_i(a|t) P_j(b|t)) and reports the Pearson X2, the
likelihood-ratio G2, df = (K-1)^2 (Chen-Thissen / mirt convention), the
chi-square p-value, Cramer's V, and the largest standardized cell residual.
Pairwise-complete cases; degenerate cells guarded.
Statistic form, G2, Cramer's V and the df were adversarially verified by a
subagent against both source papers: df = cells - independence-params =
K^2 - (2K-1) = (K-1)^2 (binary -> 1), and the reference is heuristic /
conservative (the null is stochastically smaller than chi-square), so it reads
as a diagnostic screen for residual association.
Validation:
- deterministic anchor: at K=2 the X2 equals a from-scratch 2x2 Pearson
chi-square vs the same model joint (<1e-8);
- Monte-Carlo at FITTED parameters (500 reps, N=2000): null clean pairs
X2/df=0.84 with Type I 0.027 (conservative, as the papers predict); a
2-item testlet is localized to that pair at X2/df=10.9 with power 1.00
(clean pairs stay ~1.3); a skewed population inflates all pairs (a
detectable distribution misfit). Fast CI guard + #[ignore] 500-rep study.
APA 7th references in the doc comment.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: Python API for polytomous local dependence + APA 7th refs
Expose poly_local_dependence through PyO3 and the public
local_dependence_polytomous(responses, fit) wrapper: per-pair item_i/item_j,
x2, g2, df, p_value, cramers_v, max_abs_std_resid, n_pair for a fitted
GRM/GPCM, NaN = missing. Python test covers the C(n,2) pair bookkeeping, the
df=(K-1)^2 identity, calibration under local independence, testlet detection,
and the item-count guard. APA 7th references in the Rust doc comment and
Python docstring.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: nominal categories model estimator (Thissen, Cai & Bock 2010; Bock 1972)
Add poly::fit_nominal -- unidimensional nominal-model marginal MLE with a free
scoring function a_k and intercept c_k per category,
P(Y=k|theta) = softmax_k(a_k*theta + c_k), identified by a_0=c_0=0 with
theta~N(0,1). Reuses the softmax cell (gpcm_logprobs) and its residual gradient
(gpcm_node_gradient: d/dc_k = r_k - n*P_k, d/da_k = (r_k - n*P_k)*theta) in a
Bock-Aitkin EM with a Newton M-step; the GPCM is the special case a_k = a*k.
Parameterization and identification adversarially verified by a subagent
against Thissen, Cai & Bock (2010): the a_0=c_0=0 / theta~N(0,1) baseline is
exactly the identification the authors adopt (p.45) and needs no extra scoring
contrast; only a global reflection (a_k,theta)->(-a_k,-theta) remains.
Validation:
- GPCM nesting: fitting nominal to GPCM data reaches loglik >= the GPCM fit and
recovers linear scores (a_2/a_1 ~ 2);
- recovery Monte-Carlo (500 reps, N=2000, per-item sign alignment): matched
N(0,1) ability -> score RMSE .148 / |bias| .011, intercept RMSE .104 /
|bias| .007 (near-unbiased); skewed ability -> RMSE .435 / |bias| .388
(prior-misspecification bias). Fast CI guard + #[ignore] 500-rep study.
APA 7th references in the doc comment.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: Python API for the nominal categories model + APA 7th refs
Expose fit_nominal through PyO3 and the public fit_nominal_polytomous(
responses, n_cat) wrapper returning a NominalFit (scores, intercepts, loglik,
n_iter), NaN = missing. Python test covers the (n_items, n_cat-1) shapes, GPCM
nesting (loglik >= GPCM fit, recovered scores linear a_2/a_1 ~ 2), and input
validation. APA 7th references in the Rust doc comment and Python docstring.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: polytomous person fit l_z / l_z* (Drasgow-Levine-Williams 1985; Snijders 2001)
Add poly::poly_person_fit -- the standardized log-likelihood person-fit
statistic l_z (Drasgow, Levine & Williams 1985) and its estimated-trait
correction l_z* (Snijders 2001) for a fitted GRM/GPCM, at the EAP trait. With
l_0 = sum_i log P_i(y_i), E = sum_i sum_k P_ik log P_ik, V = sum_i Var_k(log
P_ik): l_z = (l_0 - E)/sqrt(V); l_z* subtracts the log-likelihood/score
covariance (c = sumCov/sumI, tau2 = V - sumCov^2/sumI) and adds the MAP prior
score. The trait-score derivative is a central difference, so it is model
agnostic. Reuses the poly cells and score_poly_eap.
Validation:
- reduces EXACTLY to the trusted binary fitstats::person_fit l_z at n_cat=2 on
a shared EAP trait (<1e-6; l_z* to finite-difference tolerance <5e-3);
- Monte-Carlo (GPCM, 500 reps): with model respondents l_z* is ~N(0,1)
(mean -0.15, sd 1.04, Type I 0.081 -- slightly high at a 20-item test, a
documented finite-length effect), and inconsistent responders (implied trait
alternating +-1.6) are flagged with power 0.86. Fast CI guard + #[ignore]
500-rep study. APA 7th references in the doc comment.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: Python API for polytomous person fit (l_z / l_z*) + APA 7th refs
Expose poly_person_fit through PyO3 and the public person_fit_polytomous(
responses, fit) wrapper returning per-person lz, lz_star, theta_eap, and
flagged, NaN = missing. Python test covers calibration on a clean fit
(l_z* ~ N(0,1), Type I near nominal) and flagging of inconsistent responders
evaluated against that fit, plus the item-count guard. APA 7th references in
the Rust doc comment and Python docstring.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: polytomous computerized adaptive testing (Dodd, De Ayala & Koch 1995)
Add poly::poly_cat_next_item (maximum-Fisher-information selection at the
running trait estimate) and poly::poly_cat_simulate (full adaptive
administration): pick by max information, generate the response at the true
trait, re-estimate the EAP trait and posterior SD after each item via
score_poly_eap, stop at an SE threshold (or a fixed length). Composes the
already-validated poly_item_information and score_poly_eap.
Validation (procedure -> trait-recovery efficiency), 500-simulee Monte-Carlo:
- variable-length CAT (stop at SE < 0.30) recovers theta with RMSE 0.29
(normal) / 0.33 (skew) using only ~9.7 of 40 bank items;
- at a fixed length of 12, maximum-information selection beats random selection
(RMSE 0.266 vs 0.334 normal; 0.301 vs 0.401 skew). Fast CI guard + #[ignore]
500-simulee study. APA 7th references in the doc comment.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* gpcm: Python API for polytomous CAT simulation + APA 7th refs
Expose poly_cat_simulate through PyO3 and the public cat_simulate_polytomous(
true_theta, fit) wrapper returning per-simulee theta_eap, theta_sd, n_used.
Python test covers efficient trait recovery (RMSE < 0.4 using < 75% of the
bank), max-information beating random selection at a fixed length, and config
validation. APA 7th references in the Rust doc comment and Python docstring.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* Add polytomous IRT likelihood-ratio DIF (two-group marginal EM)
Detect differential item functioning for GRM/GPCM items with the
parametric IRT-LR method (Thissen, Steinberg & Wainer, 1993; Woehr &
Meriac, 2010). Fit a compact model (all items group-invariant) once,
then per studied item an augmented model (that item freed per group)
with every other item as the anchor; LR = 2*(ll_aug - ll_compact) is
referred to chi2((n_groups-1)*n_cat). Each non-reference group's latent
N(mu_g, sigma_g^2) is estimated in BOTH models (reference group pinned to
N(0,1)) so group ability differences (impact) are absorbed, not misread
as DIF.
Compute in Rust: poly::fit_poly_multigroup is a Bock-Zimowski multi-group
marginal EM whose per-item M-step reuses the single-group Newton step on
each group's nodes/expected-counts stacked (the concatenation is exactly
the Bock-Zimowski pooling); poly_dif_sweep drives the compact/augmented
fits and BH-FDR flagging. PyO3 poly_dif exposes it as the public
dif_polytomous(responses, group_id, n_cat).
Validation (500-rep Monte-Carlo, GPCM, K=3, impact focal N(0.5, 1.2^2)):
under no DIF the test is calibrated (Type I 0.042, mean LR 2.92 ~ df=3),
uniform difficulty DIF is detected with power 0.996 and non-uniform slope
DIF with 0.920, and a skewed focal population inflates Type I only mildly
(0.057). A structural test confirms the augmented fit never falls below
the compact one and recovers the focal mu/sigma.
Robustness (from an adversarial review of untested paths):
- Surface a non-finite fit as NaN (compact -> error, augmented item ->
NaN row, unflagged) instead of letting (2*(NaN-ll)).max(0.0) silently
report a DIF item as clean; can arise for GRM when a focal group has a
rarely-used category and thresholds disorder (GPCM recommended there).
- Densify group labels in the Python wrapper and reject empty declared
groups in the core so df = (n_groups-1)*n_cat counts only groups backed
by data (non-contiguous / 1-based codes no longer inflate df).
- Document effect_size as an unsigned across-group range (magnitude).
Fast CI guards plus an #[ignore] 500-rep study in Rust; Python tests cover
the sweep, label densification, and the no-silent-false-negative guard.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* Add nonparametric polytomous person-fit U3poly (Emons 2008)
Generalize van der Flier's (1982) U3 person-fit statistic to ordered
polytomous items (Emons, 2008): a model-free index that needs no fitted
IRT model. Each item-step response function P(Y_i >= m) is estimated by
its sample proportion and turned into a logit weight; a person's observed
weighted score is compared to the largest and smallest weighted scores
attainable at that person's total score (the conditioning group), giving
U3 in [0, 1] with 1 = maximally popularity-inconsistent (misfit). The
attainable min/max bounds come from exact min-plus / max-plus DP rather
than the flat "sum of the top-k step weights" shortcut, which over-counts
once an unused category breaks within-item monotonicity.
Compute in Rust: poly::u3_poly_person_fit (+ u3_poly_bootstrap_cutoff,
a parametric-bootstrap simulated critical value, since U3poly has no
usable analytic null and the normal reference is documented as
miscalibrated). PyO3 u3_person_fit / u3_bootstrap_cutoff expose the public
u3_person_fit_polytomous / u3_cutoff_polytomous.
Validation:
- exact n_cat=2 reduction to a from-scratch van der Flier U3 (max abs diff
< 1e-10), the trusted-binary correctness anchor;
- 500-rep Monte-Carlo (GPCM, K=5, n=600): the simulated cutoff calibrates
the marginal flag rate under a matched population (Type I 0.052 normal /
0.054 skew, cutoff estimated under the matching latent shape) and detects
careless responders with power ~1.00. The per-total-score-group flag-rate
deviation (0.066 / 0.083) is reported to make transparent that a single
pooled cutoff cannot fully condition on the total score.
Robustness (from an adversarial review of the missing-data paths):
- an all-missing respondent (no conditioning group) returns NaN, not a
silent 0.0 "perfect fit" that could never be flagged;
- document that the bootstrap cutoff is calibrated for complete-length
patterns and should not flag persons with substantial missingness.
Complements the parametric l_z/l_z* (person_fit_polytomous) with a
distribution-free screen. Fast CI guard + #[ignore] 500-rep study in Rust;
Python test covers the sweep, careless detection, K=2, and the all-missing
guard.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* Add observed-score equating (Kolen & Brennan)
Add a new mlsirm_core::equating module and the public
equate_observed_scores / equate_neat — the raw-score-to-raw-score
complement to the IRT scale linking in irt_link. Covers the
equivalent-groups design (mean, linear, and equipercentile equating with
the Kolen-Brennan uniform-kernel continuization, equated scores kept
real-valued) and the common-item non-equivalent-groups (NEAT) design via
chained equipercentile and frequency-estimation (post-stratification)
equipercentile.
The attainable percentile-rank inverse interpolates within the bracketing
score's interval at every point including the low boundary (needed for an
exact self-equating identity). Frequency-estimation builds synthetic
densities from the two score-by-anchor tables and renormalizes the
surviving mass, so a partially non-overlapping anchor degrades toward each
group's own marginal; a fully disjoint anchor is rejected rather than
returning a silently-collapsed table.
Compute in Rust; exposed via PyO3 and a Python equating.py (EquateResult).
Validation:
- three exact identities as correctness anchors: equipercentile
self-equate is the identity to < 1e-9 (including the low boundary at
x = 0); mean/linear recover a known integer-affine transform to < 1e-9;
both NEAT methods collapse to EG equipercentile under equal anchor
distributions to < 1e-9;
- a 500-replication Monte-Carlo against a deterministic Lord-Wingersky
population equating: the empirical equipercentile converges at the
expected rate (interior RMSE 0.53 at N=1000 -> 0.26 at N=4000, ratio
1.99 ~ sqrt(4); max bias 0.068 -> 0.031).
Robustness (from an adversarial review): FE rejects fully disjoint anchor
support; strict raw-score range validation (bin by the containing
[c-0.5, c+0.5) interval); the k=None inference path validates empty /
non-finite input; documented that FE reports synthetic-population moments.
Deferred behind the density/table interface: Tucker/Levine linear NEAT,
log-linear presmoothing, and Gaussian-kernel equating.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* Add kernel equating and log-linear presmoothing
Two enhancements to the equating module, behind its density interface.
loglinear_smooth(counts, degree) presmooths a score-frequency
distribution by Poisson-ML log-linear fitting (Holland & Thayer, 2000):
an orthonormal polynomial design on a centered/scaled score, Newton with
step-halving, reusing poly::solve_small. It preserves the first `degree`
sample moments exactly while damping sampling noise, returns AIC/BIC for
degree selection, and reproduces the raw relative frequencies when
saturated (degree = k).
equate_eg_ext adds a Gaussian-kernel continuization (von Davier, Holland
& Thayer, 2004) and optional per-form presmoothing to the equipercentile
family, behind one extended entry point whose uniform-kernel path
reproduces the existing equipercentile bit-for-bit. The Gaussian cdf
F_h(x) = sum_j r_j Phi((x - a x_j - (1-a) mu)/(a h)) is inverted by
safeguarded Newton; the bandwidth is chosen by the penalty method.
fitstats::erfc is promoted to pub(crate) for the normal cdf.
Compute in Rust; exposed via PyO3 (loglinear_smooth,
equate_observed_scores_ext) and Python (loglinear_smooth,
equate_observed_scores_kernel). EquateResult gains h_x/h_y.
Validation:
- exact-identity anchors: uniform-kernel equating equals the
equipercentile to < 1e-12; presmoothing preserves the first T moments
to < 1e-8 and, saturated, reproduces rel_freq; the Gaussian self-equate
is the identity, a large bandwidth drives kernel equating to linear to
< 1e-4, and the continuized density preserves the discrete mean/variance;
- a 500-rep Monte-Carlo against the population Gaussian-kernel transform:
interior RMSE 0.53 -> 0.26 from N=1000 to 4000 (ratio 2.03 ~ sqrt(4)),
max bias 0.049 -> 0.020.
Robustness (from an adversarial review): scale-relative Newton tolerances
(so large-N fits converge instead of spuriously reporting non-convergence);
optimal_bandwidth never returns a bandwidth worse than its grid best on a
non-unimodal penalty; the default degree clamps to k for short forms;
bandwidth validation only under the Gaussian kernel; smooth-degree input
validation.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* Add Tucker & Levine linear NEAT equating
Add equate_neat_linear — the linear observed-score equating methods for
the common-item non-equivalent-groups (NEAT) design (Kolen & Brennan,
2014, Sec. 4.3-4.4; Brennan, 2006), alongside the existing chained and
frequency-estimation equipercentile NEAT methods.
Each method forms synthetic-population moments of the two forms (weighted
by w1) from a group total-on-anchor slope gamma, then equates linearly.
Tucker uses the regression slope Cov(total, V) / Var(V); Levine uses the
congeneric effective-length ratio, which differs for an internal anchor
(Var(total)/Cov) versus an external one ((Var(total)+Cov)/(Var(V)+Cov)) —
encoded as NeatLinearMethod {Tucker, LevineObserved} x AnchorKind
{Internal, External}.
Compute in Rust (equating::equate_neat_linear); exposed via PyO3 and a
Python equate_neat_linear.
Validation:
- exact reduction to equivalent-groups linear equating under equal anchor
moments — all four Tucker/Levine x internal/external variants, for w1 in
{0, 0.5, 1}, match equate_eg(Linear) to < 1e-9;
- a hand-computed check that pins the internal-vs-external Levine gamma
(the historically error-prone formula) against an independent NumPy
oracle: the three gamma branches give three distinct conversions;
- a 500-replication Monte-Carlo under a common-regression generative model
(which satisfies the Tucker assumption by construction): the equated
table converges to the large-N reference at the 1/sqrt(N) rate (interior
RMSE 0.39 -> 0.19 from N=1000 to 4000, ratio 2.02 ~ sqrt(4); max bias
0.051 -> 0.034).
The gamma formulae and synthetic moments were independently reviewed
line-by-line against the reference equations and the oracle. Deferred:
Levine true-score equating, Braun-Holland linear.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* Add standard errors of equating (bootstrap + delta-method)
Add equating_standard_errors, the per-score-point sampling standard error
of the equated score for the equivalent-groups design (Kolen & Brennan,
2014, ch. 7), by two routes.
bootstrap_see (nonparametric bootstrap, all EG methods): resamples
examinees per group independently with replacement at the observed sample
sizes (the two forms are given to separate random samples, so their
sampling errors are independent), re-equates each of n_boot replicates
through equate_eg unchanged, and returns the per-score bootstrap SD
(divisor n_boot-1) and a percentile confidence interval. This is the only
SEE route that covers equipercentile.
analytic_see (delta-method, normal-theory, mean/linear): closed forms
Var[e_Y(x)] = sigma_x^2/n_x + sigma_y^2/n_y (mean, constant in x) and
sigma_y^2 (1+z^2/2)(1/n_x+1/n_y) (linear, z=(x-mu_x)/sigma_x), with an
asymptotic-normal interval; errors on equipercentile.
Compute in Rust (equating::bootstrap_see / analytic_see, reusing equate_eg
and nodes::inv_normal_cdf, no new dependency); exposed via PyO3 and Python.
Validation:
- the analytic-Linear SEE agrees with the bootstrap-Linear SEE within
Monte-Carlo tolerance (interior relative gap < 0.15);
- the Mean SEE is constant in x and equals the closed form;
- bootstrap SE is positive, the CI brackets the estimate, it shrinks
~1/sqrt(N), and it is seed-deterministic;
- a 500-replication Monte-Carlo confirms the bootstrap SE recovers the
TRUE sampling SD of e_Y(x) (from an outer fresh-sample Monte-Carlo) --
interior boot/true ratio in [0.95, 1.08] for equipercentile.
Both closed forms were independently re-derived and the implementation
reviewed line-by-line. Deferred: NEAT bootstrap SEE, analytic
equipercentile/kernel SEE.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
* Add lognormal response-time model (van der Linden, 2007)
Add a new mlsirm_core::rt module and the public fit_response_times: the
speed-side analogue of the 2PL for item response times, opening a
response-time modality alongside the accuracy models. For person j
(latent speed tau_j) and item i (time intensity beta_i, time
discrimination alpha_i),
ln(T_ij) | tau_j ~ Normal(beta_i - tau_j, 1/alpha_i^2),
estimated by marginal-ML EM with tau ~ Normal(0, sigma_tau^2), speed
scored by EAP.
Because the model is conditionally Gaussian with a unit loading on tau,
the speed posterior, the marginal likelihood, and the EAP are all exact
closed forms (matrix-determinant / Sherman-Morrison), so the estimator
needs neither Gauss-Hermite quadrature nor a line search -- the EM is
exact O(nnz) coordinate ascent, with the alpha M-step carrying the +v_j
posterior-variance correction (dropping it biases alpha high). The
log-time metric identifies the speed scale, so sigma_tau is estimated and
only the location is pinned (mu_tau = 0); alpha multiplies the residual,
not beta - tau, so there is no alpha<->sigma_tau trade-off. Compute in
Rust; PyO3 + Python wrapper; missing/non-positive times are marginalized
per person.
Validation:
- an exact identity anchor: the closed-form (Woodbury) marginal
log-likelihood equals a dense multivariate-normal log-pdf to < 1e-9,
certifying ln|Sigma|, the quadratic form, and all sign conventions;
- a reduction anchor: sigma_tau -> 0 collapses to the per-item lognormal
MLE (beta = mean log-time, 1/alpha^2 = var log-time);
- a 500-replication Monte-Carlo (N=800, 20 items, ~30% missing) under
normal AND a misspecified skew speed population: the item parameters
stay essentially unbiased (RMSE alpha 0.067 / beta 0.027, bias beta
-0.0001 under skew) with speed recovered at corr 0.92 -- the level-1 RT
item pa…
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.
Summary
Model-design PR implementing marginal maximum likelihood (Bock-Aitkin EM) for the full latent-space model family, estimation-level population structures, a wgpu GPU E-step, all-Rust scoring/fit-statistic compute, and the implementable core of four supplied literature batches. Design doc:
docs/mmle_marginal_lsirm_design.md; per-equation paper basis:docs/papers/mmle-lsirm-formula-compilation.md(+group_{a,b,c}_specs.md,implemented-literature-map.md,corpus-triage-batch{3,4}.md).Estimator (
mlsirm_core::marginal)MIRT/MLS2PLM/MLSRM/ULS2PLM/ULSRMunder the simple-structure contract; deterministic Gauss-Hermite E-step made tractable by the conditional factorization (Q_xi^K * sum_d Q_theta, notQ^(1+D+K)); Fisher-preconditioned GEM M-step; Jeon et al. (2021) LSIRM priors as MAP penalties (PenaltyConfig::lsirm_prior).sigma_u/ICC) populations; FIPC anchors (Kim 2006, MWU-MEM-style) withPopulationSpec::SingleFree; concurrent calibration = multigroup + structural missingness (Hanson-Beguin), with a recovery test.xi_rule = gh | qmc | mc): Halton (+ Cranley-Patterson shift) and seeded-MC node sets (Jank 2005; Wei-Tanner 1990; Meng-Schilling 1996), liftinglatent_dim <= 3for the grid rule.zero_inflation=True; Perumean-Chaney 2013 guidance) and a context-varying item covariate with estimated coefficient (Debeer-Janssen 2013 position effects).gpu_marginal.rs): ~110 s -> ~5 s per multilevel E-step iteration on a 31k x 57 dataset (RTX 3050 Ti); race-free slot-ownership reductions; M-step and final EAP stay CPU f64; 1 GiB buffer guard with CPU fallback.Scoring & serving (
scoring.rs,serving.py)fast-mlsirm scoreCLI.Fit statistics & validation (
fitstats.rs,agreement.rs,oakes.rs)rms_residualpractical-significance guard (Sinharay-Haberman 2014),l_z/l_z*(Snijders, prior-mean-centeredr_0), infit/outfit, BH-FDR, information criteria (Kang-Cohen-Sung 2009), Vuong tests (Schneider et al. 2019), Q3/GDDM (Svetina-Levy 2014), residual item fit (Haberman et al. 2013), adjusted pairwise chi2/df (Tay-Drasgow 2012), Chen-Thissen (1997) LD indices, resampling person fit (Sinharay 2016), TCC drift (Guo et al. 2015).validate_judge): QWK/r/SMD/degradation/subgroup-SMD with the paper's thresholds.select_items()screening pipeline,dif_analysis()(virtual-item LR; Jeon-Rijmen-Rabe-Hesketh 2013 + Makransky-Glas 2013),irtree_expand()(Jeon-De Boeck 2016).Contracts
tests/test_marginal_parity.py).mlsirm-core; Python is orchestration + parity references.cargo test --workspace+ binding-crate manifest run).estimator="mmle"fast path for plainULS2PLM/ULSRMunchanged; spatial models now fit instead of raising.Real-data validation (31,356 x 57, 307 teams)
Zero-inflated multilevel MLS2PLM selected by BIC (Δ=310 vs plain multilevel; Δ=2,336 vs single); ICC = 0.80; π = 0.21; corrected screening retains 21/57 items with a full audit; serving bundle + EAPsum tables exported and smoke-tested end to end (EAP/MAP/EAPsum).
Explicitly-scoped follow-ups (docs/papers/corpus-triage-batch4.md, priority order)
BIFAC2PLMbifactor family (Gibbons-Hedeker reduction slots into the existing conditional-factorization E-step)🤖 Generated with Claude Code