fix: SE cumulative and EUR missing gamma(1/n) factor

[mwburgoyne]

SE._Nfn returns qi*tau/n * gammainc(1/n, (t/tau)^n). Since scipy.special.gammainc is the regularised lower incomplete gamma P(a,x), the closed-form integral of the SE rate qi*exp(-(t/tau)^n) needs the extra gamma(1/n) factor:

N(t) = qi*tau/n * gamma(1/n) * P(1/n, (t/tau)^n)

Without it, SE cumulative volume and EUR are off by a factor of Γ(1/n): understated for n < 0.5 (about 25% at n=0.4, 64% at n=0.3), exact at n=0.5, overstated by about 10% for n > 0.5.

Verified against adaptive quadrature of the rate (scipy.integrate.quad); a regression test test_SE_cum_matches_integral is added, since the previous check_model assertions only tested that cum was finite, so a constant-factor error passed silently.

For very small n, gamma(1/n) overflows (the closed-form EUR genuinely diverges there); the code guards on a finite coefficient and falls back to the existing numerical integrator.

🤖 Generated with Claude Code

[dsfulf]

Review — verified, the math is correct.

Traced the derivation and confirmed numerically against scipy.integrate.quad. With the substitution u = (t/tau)^n, the integral of qi*exp(-(t/tau)^n) is qi*tau/n * gamma_lower(1/n, (t/tau)^n), where gamma_lower is the unregularized lower incomplete gamma. Since scipy's gammainc returns the regularized P = gamma_lower / Gamma, the closed form requires the Gamma(1/n) factor. Empirical check (qi=1000, tau=30, t=5000):

n     new/quad   old/quad   Gamma(1/n)
0.3   1.00000    0.3600     2.7782
0.4   1.00000    0.7523     1.3293
0.5   1.00000    1.0000     1.0000   (Gamma(2)=1, unaffected)
0.6   1.00000    1.1077     0.9027
0.8   1.00000    1.1033     0.9064

Notes:

• The bug predates the v1.3.0 integration work — git blame puts the faulty line at 57b68ac (Aug 2021). The perf refactor never touched it (SE.cum is a closed form, not numerically integrated), so that work neither introduced nor caught it.
• It went undetected because check_model asserts only rate(0)==qi, cum(0)==0, and finiteness — every monotonicity/value check is commented out, so a constant-factor error is invisible to it. The new test_SE_cum_matches_integral (compare cum to quadrature of rate) closes exactly that gap.
• The overflow guard is sound: Gamma(1/n) overflows only for n < ~0.0058, and the fallback receives a validated ndarray from cum().
• Real-world impact: any SE fit with n != 0.5 had EUR wrong — understated up to ~64% (n=0.3), overstated ~10% (n>0.5). This is being called out as a breaking numerical change in the version history.

Nit: three trailing blank lines added in test/test_perf.py.

[dsfulf]

Verified the math: the closed-form SE cumulative requires the gamma(1/n) factor (scipy gammainc is the regularised P). Confirmed new/quad = 1.00000 across n in [0.3, 0.8]; the small-n overflow fallback to the numerical integrator matches quadrature to ~5e-9. Rebased onto main (post-#21) so the mypy unreachable-statement lint passes; full matrix green. Approving.