Changes for version 0.292 - 2026-08-05
- fixed long-double bug https://www.cpantesters.org/cpan/report/506975f6-906a-11f1-8f30-a201c4f2440e
- `power_t_test` was cross-validated against R 4.6.1 `power.t.test` and against `scipy.stats.nct` driven by `scipy.optimize.brentq`, over a grid of 288 cases covering all five solved-for parameters, all three types, both alternatives and `strict`. Three fixes came out of it:
- The Simpson sum behind the noncentral *t* CDF put a fixed 30000-step grid on `u = w/(1+w)`, and the chi density it integrates defeats that at both ends. The density carries `w**(df-1)`, so unless `df` is a whole number some derivative of it is infinite at `w = 0` and Simpson's error bound does not hold: two good digits at `df = 1.2` with `sig_level = 1e-4`, five at `df = 1.2`, nine at `df = 1.8`. Substituting `w = z**m`, with `m` chosen so that `m*df - 1 >= 3`, restores the bounded derivatives and brings all of those to machine precision. It also puts the origin's contribution at zero, which subsumes a separate bug: the sum had been dropping its `u = 0` endpoint term, worth 7e-7 of absolute power at `df == 1`. `nu` is now also floored the way R floors it, per sample rather than in total.
- The same density has standard deviation `1/sqrt(2*df)` and so narrows without bound, while the grid did not. Past `df` of about 1e7 the steps went clean over the peak: `power_t_test(n => 4e7, delta => 0)` returned 0.138 where the answer can only be `sig_level/2`, and a large-cohort `n` solved 9% low. Above `df` of 1e3 the steps now go on `w` across +/- 12 standard deviations of the mode, with the chi normalisation taken from Stirling's series to keep the peak height from cancelling away; and above 4e5, where those log terms cancel too hard for any grid to help, the Abramowitz & Stegun 26.7.10 asymptotic form takes over -- the same formula, at the same cut-off, that R's `pnt.c` uses. That is also 25 times quicker than integrating.
- The power was formed as `1 - P(T <= t)`, which loses most of its digits to cancellation when the power is small. It is now integrated as the upper tail directly.
- The four inverse solvers were plain bisection stopped at the bracket width, which capped `n`, `delta`, `sd` and `sig_level` at R's own four or five significant figures. They now use regula falsi with the Illinois correction against a relative tolerance, so they match machine-precision `brentq` roots to ~1e-13 in fewer evaluations than the bisection took. The `tol` default moved from `1.22e-4` to `1e-12` to match.
- Nothing checked that the bracket held a root, so an unreachable target came back as a bracket endpoint wearing the requested power: solving for `sd` with `power => 0.01` returned `delta * 1e7`, and with a negative `delta` returned a negative standard deviation. Unreachable targets now croak and name the range searched. `sig_level` and `power` outside `[0, 1]`, an `n` below 2, a negative `sd`, and an unrecognised `type` or `alternative` are rejected as well -- `type => 'twosample'` used to be read silently as `'two.sample'`.
- New test file `t/power_t_test.R.scipy.t` carries the cross-validated grid.
Modules
Get basic statistical functions, like in R, but with Perl using XS for performance