[Timefreq] Port correct_mc random-symbol fitting

[suraj-ranganath]

Port the EEGLAB Ramberg-Schmeiser helpers used by correct_mc and replace the temporary neighboring-bin p-value estimator with standalone rsfit. Bundle correctfit with the same distribution-fitting family, without expanding bootstat distfit. Part of epic #131 and stacked under PR #152.

Closes #147

[claude]

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

---

PR Review: [Timefreq] Port correct_mc random-symbol fitting

• Gather context (read changed files, EEGLAB reference)
• Review parity of ported helpers (rsfit, rsadjust, rsget, rspdfsolv, rspfunc, correctfit)
• Review correct_mc integration changes
• Review tests and docs
• Post review feedback

Note: the vendored EEGLAB checkout (src/eegprep/eeglab/) isn't present in this CI environment, so the parity findings below are grounded in the Ramberg–Schmeiser moment definitions and in internal consistency between the files in this PR, rather than a line diff against the MATLAB source. The @pytest.mark.matlab parity tests are skipped here.

Code review

• Overall assessment: needs changes — one likely sign bug in rsadjust, plus a test gap that lets it slip through.
• Highest-risk area: rsadjust λ2 sign selection (affects fitted distribution orientation in correct_mc).
• Merge recommendation: needs changes.

Blocking

None.

Important

rsadjust.py:30 — third-moment numerator uses -2*a*b instead of -3*a*b.

The sign of l2 is chosen from the standardized third moment of the RS distribution, whose numerator is C − 3AB + 2A³. This PR's own rspdfsolv.py:52 computes that exact quantity as c - 3.0 * a * b + 2.0 * a**3, but rsadjust uses c - 2.0 * a * b + 2.0 * a**3:

# rsadjust.py:30
if skew * (c - 2.0 * a * b + 2.0 * a**3) < 0:
    l2 = -l2

When a ≠ 0 and the two expressions differ in sign, l2 (and therefore l1 = mean - a/l2) flips, mirroring the fitted distribution and producing wrong tail p-values. correct_mc fits neighboring-bin correlations that are typically positive/skewed, so this path is exercised in practice. Why it matters: silently wrong multiple-comparison p-values, with no error raised.

Suggested fix (align with rspdfsolv):

if skew * (c - 3.0 * a * b + 2.0 * a**3) < 0:
    l2 = -l2

Please confirm against functions/timefreqfunc/rsadjust.m in the EEGLAB source before merging. Fix this →

Nits

None.

Test gaps

rsadjust is only tested on the symmetric a = 0 case, which cannot catch the bug above. Both test_ramberg_schmeiser_helpers_cover_analytic_cases and the MATLAB parity test call rsadjust(1.0, 1.0, 0.0, 1/12, 0.0). For l3 == l4 we get a = 0, so the c - …a·b… + 2a³ term reduces to c and the -2ab/-3ab distinction vanishes — the sign branch is never discriminated.

Suggested test: add an asymmetric case with l3 != l4 and nonzero skewness (e.g. rsadjust(0.2, 1.5, mean, var, skew)), asserting the sign of l2. Ideally cross-check the full rsfit round-trip on a skewed sample so correct_mc's real path is covered, not just the symmetric uniform case in test_ramberg_schmeiser_helpers_cover_analytic_cases (lambdas[[0,2,3]] only).

EEGLAB parity notes

• rspfunc, rsget, rspdfsolv (incl. the l3*l4 < 0 ×2 penalty), and rsfit's moment/shape/adjust flow match the expected EEGLAB structure. The added _fit_shape_parameters fallback initial guess [-0.1, -0.1] and the variance_term <= 0 guard in rspdfsolv are reasonable robustness additions that only engage when the primary path fails to converge.
• correctfit.py reworks EEGLAB's (pval, fitstruct) signature into keyword args (allpval/gamparams/zeromode). It's self-consistent and tested, but it's not on the correct_mc rsfit path this PR targets — worth a one-line note in the PR description that it's bundled as part of the correction family rather than required by the issue surface.
• correct_mc integration (lag-1 neighbor correlation → rsfit(corrs, 0.0) → Bonferroni threshold) is internally consistent and matches the new tests.

[suraj-ranganath]

🤖 Specification

Problem: correct_mc previously returned standalone neighboring-bin p-values that did not use EEGLAB's Ramberg-Schmeiser fitting path, leaving the issue #147 rows partial/deferred.

Approach: add standalone ports of rsfit, rsadjust, rsget, rspdfsolv, rspfunc, and correctfit under functions/timefreqfunc; wire correct_mc to rsfit(correlations, 0) while preserving its existing input/output workflow; update exports, docs, matrix rows, and implementation notes.

Key code: rsfit computes sample moments, solves lambda3/lambda4 with rspdfsolv, adjusts lambda1/lambda2 through rsadjust, then calls rsget for the fitted cumulative probability. rsadjust intentionally preserves EEGLAB rsadjust.m's C - 2AB + 2A^3 lambda2 sign convention.

Tests: focused tests cover correct_mc delegation, analytic RS cases, asymmetric rsadjust parity behavior, correctfit gamma-parameter correction, the full time-frequency statistics file, matrix validation, ruff, ty, and PR-scope pre-commit. MATLAB parity coverage is present and skip-aware when matlab.engine is unavailable locally.