[Parity] Implement statistics package

[suraj-ranganath]

Add an EEGLAB-style statistics package with NumPy/SciPy helpers for statcond, t-tests, ANOVA, FDR, surrogate p-values and confidence intervals. Includes focused tests, MATLAB Engine parity coverage, API docs, implementation notes, and parity matrix updates.

Fixes #134
Part of #131

[suraj-ranganath]

🤖 Specification

Problem: EEGPrep had no runtime equivalent for EEGLAB functions/statistics, leaving docs/parity/eeglab_core_parity_matrix.json Phase 3 rows open and forcing Phase 4/6 work to risk private statistics copies. The missing surface included FDR, surrogate p-values/CIs, statcond, cell-array t-tests/ANOVA, correlation, and concatdata.

Approach: Added src/eegprep/functions/statistics/_core.py with NumPy/SciPy implementations and same-name wrapper modules for EEGLAB discoverability. Public helpers use explicit final-axis case/subject/surrogate contracts; statcond returns named fields while preserving tuple unpacking. Two-way ANOVA uses named row/column/interaction effects to avoid MATLAB output-order ambiguity.

Key code: FDR and surrogate helpers live in _core.py lines 104-218; surrogdistrib and statcond live in _core.py lines 504-631. The same-name modules are thin exports so future STUDY/time-frequency code can import stable package APIs without copying formulas.

Tests: tests/test_statistics_package.py covers FDR, surrogate p-values/CIs, t-tests, ANOVA, statcond dispatch, deterministic resampling, invalid inputs, and teststat. MATLAB Engine parity tests compare deterministic helpers against the vendored EEGLAB statistics folder with ANOVA tolerance rtol=1e-4/atol=1e-5 for MATLAB single-precision intermediates.

[claude]

Claude finished @suraj-ranganath's task in 3m 59s —— View job

---

Code review

• Overall assessment: Needs changes (one correctness bug + one repo-hygiene issue; otherwise solid)
• Highest-risk area: stat_surrogate_pvals two-tailed p-value
• Merge recommendation: Needs changes

Reviewed the new functions/statistics package: _core.py numerics (FDR, surrogate p-values/CIs, t-tests, ANOVAs, statcond, surrogdistrib), the same-name wrapper modules, tests, docs, and the parity-matrix updates. The architecture is clean — dataclasses with __iter__ to preserve EEGLAB-style tuple unpacking, explicit last-axis case/surrogate contracts, and named two-way effects to dodge MATLAB output-order ambiguity. Two items below should be addressed before merge.

Blocking

None.

Important

1. stat_surrogate_pvals two-tailed branch can return p > 1 — _core.py:177
When the observed value ties with surrogate samples, p_right and p_left can both exceed 0.5, so 2 * np.minimum(p_right, p_left) exceeds 1. Verified: stat_surrogate_pvals([[1,2,2,3]], [2], "both") → 1.5. EEGLAB clamps the two-tailed value to 1, and a probability > 1 is invalid regardless. The unit test only uses inputs that land exactly on 1.0, so it misses this. Fix: return np.minimum(2 * np.minimum(p_right, p_left), 1.0) and add a tie case. (Inline suggestion posted.) Fix this →

2. .notes/implementation-notes.html is committed but .gitignore excludes it — .gitignore has .notes/* and whitelists only .gitkeep and eeglab-migration-gap-audit.md. This HTML file is not whitelisted, so it appears to have been force-added by accident. Remove it from the PR: git rm --cached .notes/implementation-notes.html.

Nits

• fdr (_core.py:139-141) computes correction before validating fdr_type, so an invalid type does the harmonic-sum work and only then raises. Reorder the validation above the correction line for clarity. Low priority.
• statcond (_core.py:622-624): passing alpha with method="param" raises because no surrogate exists. That's a reasonable constraint, but the docstring for alpha doesn't mention it requires a nonparametric method — worth a one-line note.

Test gaps

• No test covers the two-tailed p > 1 tie case (Important refactor for packaging and test distribution #1) — add one alongside the clamp fix.
• The nonparametric two-way statcond surrogate path is only checked for shape/seeding (test_nonparametric_statcond_and_surrogdistrib_are_seeded), not numerical correctness. The deterministic helpers have MATLAB parity, but the surrogate two-way ANOVA assembly itself is unverified. Acceptable for this PR; flagging so it isn't assumed covered.

EEGLAB parity notes

• The MATLAB parity tests (@pytest.mark.matlab) addpath to src/eegprep/eeglab/functions/statistics, which isn't present in this checkout — they skip locally and only run where the vendored EEGLAB reference + MATLAB Engine exist. The two-way ANOVA output-order mapping (anova2_cell first output → columns, second → rows) is encoded only in that gated test. Make sure that MATLAB lane actually runs somewhere before relying on the mapping, since a swapped row/column assignment would otherwise pass all the runnable tests.
• teststat.py intentionally ports EEGLAB's teststat demo/test script as a Python smoke helper rather than a numeric function — the parity matrix documents this, and it's a reasonable call.
  · branch phase/3-statistics-package

[claude]

Severity: important — the two-tailed branch can return p-values greater than 1.

When the observed value is tied with surrogate samples, p_right and p_left can both exceed 0.5, so 2 * min(p_right, p_left) exceeds 1. Example: stat_surrogate_pvals([[1,2,2,3]], [2], "both") gives p_right=0.75, p_left=0.75 → 1.5. EEGLAB clamps the two-tailed p-value to 1; a probability above 1 is also invalid on its own terms. The existing unit test (test_surrogate_pvals_and_ci_use_last_axis) only exercises inputs where the result lands exactly on 1.0, so it misses this.

Suggested change

	return 2 * np.minimum(p_right, p_left)
	if tail_name == "both":
	return np.minimum(2 * np.minimum(p_right, p_left), 1.0)

Please also add a tie case to the test.