/ 
Stats-LikeR-0.3

Changes for version 0.3 - 2026-08-16

  • shapiro_test
    • `shapiro_test` rebuilt against R 4.6.1's `src/library/stats/src/swilk.c` — AS R94, Royston (1995) — driven by R's and SciPy's own test suites rather than by cases invented here. Four bugs are fixed, one of them a case R's regression suite tests for by name, and the statistic is now more accurate than R's own on a sample whose values dwarf its spread.
    • Everything below is checked in the new `t/shapiro_test.R.scipy.t` (146 tests), whose expected values are frozen literals with their provenance recorded in the file header; it needs no R and no Python to run. The generators that produced them, `t/shapiro_test.R.scipy.R` and `t/shapiro_test.R.scipy.py`, are committed beside it. The full suite is 124 files and 25,081 tests, and `./test.all.perls.pl` passes on all five local perls — `5.10.1`, `5.12.5` (long double), `5.42.3`, `5.44.0` and `5.44.0-quadmath` — with no warnings on any of them.
    • The p-value could come back negative:
    • R's `tests/reg-tests-1b.R` contains exactly one `shapiro.test` assertion, and it is this:
      • stopifnot(shapiro.test(c(0,0,1))$p.value >= 0)
    • `shapiro_test([0,0,1])` returned `-4.6648135328131477e-15`. At `n = 3` the p-value is `6/pi * (asin(sqrt(W)) - asin(sqrt(3/4)))` and `W` has an exact floor of `3/4` that `c(0,0,1)` sits on, so the subtraction lands on zero from whichever side the constants round to; R clamps the result at 0 and this module did not. It is the same defect SciPy fixed as gh-18322. The clamp is in, and `asin(sqrt(3/4)) = pi/3` is now carried to NV width rather than R's 15 digits, so the same case comes out at `+4.2e-16` before clamping instead of below zero.
    • W and the p-value were only good to nine digits:
    • The expected normal order statistics that AS R94 weights the sample with came from `inverse_normal_cdf()`, which is Moro's approximation and good to about `1e-9`. They go straight into `W`, so nine digits there is nine digits in the answer — where R reports sixteen. Against the values SciPy pins in `TestShapiro`, every one of them annotated upstream as *"reference values generated using R shapiro.test"*:
    • | SciPy case | W was | W is | R 4.6.1 | |---|---|---|---| | `test_basic` x1 | `0.900472879324135` | `0.900472879317561` | `0.90047287931756` | | `test_basic` x2 | `0.959026945965277` | `0.959026946032345` | `0.95902694603234` | | `test_basic2` x4 | `0.834666275331324` | `0.834666275318169` | `0.83466627531817` |
    • and the p-values with them:
    • | SciPy case | p was | p is | R 4.6.1 | |---|---|---|---| | `test_basic` x1 | `0.0420895752342124` | `0.0420895752222577` | `0.04208957522226` | | `test_basic` x2 | `0.524597929157127` | `0.524597930470668` | `0.5245979304707` | | `test_basic2` x4 | `0.000913490482316994` | `0.000913490481812984` | `0.000913490481813` |
    • Moro's value is still the starting point, but Newton against the `erfc`-based normal CDF finishes it. The loop stops as soon as another pass could not move the answer, so a `double` build pays for one refinement and only the wider NVs pay for a second. Across a 180-sample sweep over normal, uniform, exponential, log-normal, Cauchy, tied, tiny-scale and grid data at every n from 3 to 5000, the worst remaining disagreement with R is `2.0e-15` in `W` and `2.6e-12` in the p-value — the latter is not sloppier arithmetic but the same last bit amplified, since the p-value is a function of `log(1 - W)` over a sigma of about `0.6`.
    • 1 - W was formed by subtracting from 1:
    • The p-value depends on `log(1 - W)`, and `W` runs to within `1e-5` of 1 on a large normal sample, so computing `W = b^2/ssq` and then `1 - W` throws away exactly the digits the p-value is made of. R does not do this — its `swilk.c` forms `w1 = (ssassx - sax) * (ssassx + sax) / (ssa * ssx)` directly and says so in a comment — and now neither does this module.
    • More accurate than R when the values dwarf their own spread:
    • R's `swilk.c` divides the sample by its range but never centres it, so `1e9 + noise` loses most of its significant digits before `W` is ever formed. SciPy filed the same complaint from the other end as gh-14462 and works around it by subtracting the median; this module now does that too, which costs one subtraction per value.
    • Measured against a 60-digit `mpmath` evaluation of AS R94 on the identical doubles, over `1e6 + noise` and `1e9 + noise` at every n from 3 to 5000:
    • | | worst relative error in W | worst in the p-value | |---|---|---| | R 4.6.1 | `1.9e-8` | `2.6e-7` | | `shapiro_test` | `1.0e-15` | `1.4e-13` |
    • On well-conditioned samples the two still agree to the last few ulp, so this is a divergence only where R has already lost the digits. `t/shapiro_test.R.scipy.t` asserts the invariance rather than R's number there, and records why at that section.
    • Faster as well:
    • The sort now goes through the module's own introsort rather than `qsort()`, whose comparator the compiler cannot inline; the order statistics are generated for half the sample and mirrored, since the weights are antisymmetric; and ten `pow()` calls became Horner evaluations. `pow()` is a `__float128` call on a quadmath perl.
    • | n | was | is | |---|---|---| | 10 | 0.83 µs | 0.78 µs | | 100 | 4.78 µs | 5.05 µs | | 1000 | 79.3 µs | 49.8 µs | | 5000 | 578 µs | 459 µs |
    • `n = 100` is the one size that got slower: at that length the accurate quantiles are most of the work and there is not enough sorting to pay for them. That trade was taken deliberately.
    • One documented value was wrong:
    • The hash printed under `shapiro_test` in this README, in `read.me.pod` and in the module's own POD showed `statistic 0.960870680168535` and `p.value 0.589650577093106` — the pre-fix numbers, and for `[1..19]` while the example above them calls `shapiro_test([1..5])`. It now shows what `[1..5]` actually returns, `0.986762155447719` and `0.96717393596804`, which is R's `shapiro.test(1:5)` to the last digit R prints.
  • quantile
    • Interpolation ran between order statistics that were equal:
    • R's type 7 interpolates only when the index falls strictly between two order statistics *that differ* — `index > lo & x[hi] != qs` in `quantile.default`. This module always evaluated `(1 - g) * x[j] + g * x[j+1]`, which does not return `v` when both sides are `v`. On a two-valued sample at `n = 999` it reported `0.99999999999994` for `1`, and on the 602 identical values R's PR#16672 was filed about it failed to return that value at every prob — which is the monotonicity failure the PR is about. Both now match R exactly.
    • Probabilities a hair outside [0, 1] were refused:
    • A probability arrived at by arithmetic rather than written down can land just outside the interval. R allows `100 * .Machine$double.eps` of overshoot and clamps to the endpoint — its PR#17891, `quantile(0:1, 1+1e-14) == 1` — where this module raised an error. It now clamps within the same allowance and still errors on anything further out. R's constant is used rather than `NV_EPSILON` on purpose: it is part of what the function *accepts*, so a long-double or `__float128` build must not reject a `probs` vector R takes.
    • Faster:
    • The sort was `qsort()` with a function-pointer comparator; it is now the same introsort `shapiro_test` uses. Ordering 5000 NVs costs about 61,000 comparisons, and paying for an indirect call on every one of them is most of what a sort of that size costs.
    • | n | was | is | |---|---|---| | 100 | 3.2 µs | 2.7 µs | | 1000 | 52.5 µs | 22.8 µs | | 10,000 | 1.07 ms | 0.72 ms | | 100,000 | 13.7 ms | 8.8 ms |
    • Both fixes and the sort are covered by the new `t/quantile.R.t` (197 tests, or 205 under `EXTENDED_TESTING`), built on the two assertions R's own suite makes about `quantile` — that `quantile(x, ((1:n)-1)/(n-1))` recovers `sort(x)`, and that it equals the type-7 interpolation computed by hand off the sorted sample — run over seven input shapes chosen to break a quicksort (sorted, reversed, organ pipe, two-valued, tie ladder, sawtooth) at every n either side of the insertion-sort threshold and the recursion depth limit, plus PR#16672 and PR#17891 verbatim and 79 frozen R value tables. Its generator, `t/quantile.R.R`, is committed beside it.
  • Documentation
    • Illustrations for three more functions, drawn by `t.test.plots.pl` and `skew.kurtosis.plots.pl`, both committed:
    • **`t_test`** gains six: what the estimate, the standard error and the null distribution are and which area of it the p-value is; how `conf_int` is the estimate plus or minus a t quantile and how `conf_level` sets that quantile; the three `alternative`s side by side with the region each counts and the interval that goes with it; `p_value` as a function of `mu`, crossing `1 - conf_level` exactly at the two bounds of `conf_int`; paired, `var_equal` and Welch on the same data, with the Welch degrees of freedom as the two spreads separate; and two distributions separating with the interval retreating from `mu` as the p-value falls.
    • **`skew`** gains a left-tailed, a symmetric and a right-tailed sample against the same `N(0, 1)` curve, with the mean and median of each, which is what the sign of the statistic is reporting.
    • **`kurtosis`** gains a flat-shouldered, a normal and a heavy-tailed sample with the tails behind each drawn out, since it is the tails and not the peak that the statistic is measuring.

Documentation

CONTRIBUTING.md
README.md

Modules

Stats::LikeR
Get basic statistical functions, like in R, but with Perl using XS for performance

Other files