/ 
Stats-LikeR-0.301

Changes for version 0.301 - 2026-08-21

  • there are numerous additions of `restrict` keywords, which may or may not improve speed
  • kruskal_test
    • `kruskal_test` cross-validated against R 4.6.1's `stats::kruskal.test()` and SciPy 1.18.0's `TestKruskal`, driven by those suites' own cases rather than by cases invented here. Six bugs are fixed: five in how the arguments and the data are read before any ranking happens, and one in the chi-squared tail, which reaches every function that uses it. The test now needs a third of the memory and runs in under half the time, and it returns the same answer twice in a row, which it did not before.
    • Everything below is checked in the new `t/kruskal_test.R.scipy.t` (870 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 generator that produced them, `t/kruskal_test.R.scipy.R`, is committed beside it. There had been no `t/kruskal_test.t` at all — the only coverage was six assertions in `t/01.t` on the single Hollander & Wolfe example, and none of the six bugs would have shown up in it. The full suite is 125 files and 25,951 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.
    • NaN was ranked instead of dropped:
    • `looks_like_number` is true for `NaN`, so a `NaN` went into the ranking. R treats `NaN` as `NA` and `complete.cases(x, g)` removes it before `rank()` ever sees it:
    • | data | was | R 4.6.1 | |---|---|---| | `c(1,2,3,4,5,6)` with one `NaN`, n = 7 | `H = 4.5` | `H = 3.8571428571428577` | | `1:24` with one `NaN` | `H = 17.28` | `H = 16.5` |
    • It also handed `cmp_nv3` a comparison that is never true for any pair involving the `NaN`, which leaves `qsort` without the strict weak ordering it is entitled to — the same defect `wilcox_test` had fixed in 0.298, where the comment describing it is still in the file. `+Inf` and `-Inf` are neither `NA` nor `NaN` to R and a rank test has no trouble with them, so they are still kept and ranked; that is now pinned rather than incidental.
    • The bad value propagated: `table_one`'s `_t1_cont_p` hands its groups straight to `kruskal_test`, so a single `NaN` among nine observations was reported at `p = 0.0273` where dropping it gives `0.0439`.
    • A group with no data inflated the degrees of freedom:
    • The hash-of-arrays form counted every key in `k`, including a key whose array was empty or whose every element had been dropped. On `{a => [1,1,1], b => [2,2,2], c => []}` that gave `df = 2` and `p = 0.0820849986238988` for what is a two-group problem with `df = 1, p = 0.025347318677468304`. Such a group was already skipped when forming the statistic and when building `group_stats`; only `df` still counted it.
    • R refuses this case outright — `all groups must contain data` — and so does `kruskal_test` now, because the alternative is to test the groups that do have data under a `df` that counts one that does not. R's order of checks came with it: it filters each group, refuses an empty one, and only then counts what is left, so `{a => [], b => []}` is `all groups must contain data` and not `not enough observations`. SciPy takes the other side of this and returns `NaN` with a `SmallSampleWarning`; the divergence is recorded in the test file rather than papered over. The `x`/`g` form cannot reach any of this — it mints a group id the first time an observation survives the filter — so nothing changes there.
    • Group labels were truncated at a NUL and lost their UTF-8 flag:
    • The `x`/`g` path read the label with `SvPV_nolen` and then took `strlen` of it. Perl strings are counted, not NUL-terminated, so `"a\0X"` and `"a\0Y"` collapsed into one group: `kruskal_test([1..6], ["a\0X","a\0X","a\0Y","a\0Y","b","b"])` came back as two groups with `df = 1, H = 2.4` instead of three with `df = 2, H = 4.571428571428573`. Both paths also copied the label's bytes while dropping perl's UTF-8 flag when storing it into `group_stats`, so a label outside latin-1 came back as mojibake and the two input paths disagreed with each other about labels inside it. The length now travels with the string and carries the flag in its sign, which is `hv_store`'s own convention, so a label comes back `eq` to what went in. Dropping the `strlen` also drops a pass over every label.
    • A trailing named argument read past the argument stack:
    • The named-argument loop took `ST(arg_idx + 1)` without checking that there was one, so an odd argument list read one slot past the top of the stack — and what it found there changed which branch ran: `kruskal_test(\%h, 'x')` came back complaining that `'h'` cannot be mixed with `'x'`/`'g'`, because `x_sv` had been assigned whatever was past the end. `binom_test`, `chisq_test`, `fisher_test`, `wilcox_test`, `var_test` and `prcomp` all guard this; `kruskal_test` was the one that did not. It now croaks `odd number of named arguments`.
    • An infinite chi-squared statistic gave no p-value:
    • `get_p_value` short-circuits a statistic at or below zero and otherwise goes to `igamc`. `+Inf` is neither, so it reached the continued fraction, where the first `1/d` is `1/Inf = 0` and then `del = 0 * Inf` is `NaN` — an overwhelmingly significant result reported as no result at all. R's `pchisq(Inf, df, lower.tail = FALSE)` is `0`, and so is this now. `NaN` in gives `NaN` out, as R does, rather than running the continued fraction to its full 10,000-iteration safety bound first.
    • This is reachable from `kruskal_test`. When a sample has no variation at all the tie correction is `(n^3 - n)/(n^3 - n)`, and once `n^3` is past `2^53` the subtraction of `n` is lost from one side or the other, so an inexact zero is divided by an exact zero. R has the same problem and returns `+Inf`, `-Inf` or `NaN` depending on which way `n` rounded: `NaN` at `n = 250000`, `-Inf` at `300000`, `NaN` at `400000`, `+Inf` at `500000`, `NaN` at `750000`, `+Inf` at `1000000`, `NaN` at `1500000` and `+Inf` at `2000000`. `kruskal_test` now agrees with it on all eight. Below that the correction is exact and both give `NaN`, which the corpus pins. `get_p_value` is shared, so `chisq_test`, `prop_test`, `mcnemar_test`, `friedman_test`, `cmh_test`, `logrank_test` and `coxph` get the same fix.
    • Three times less memory, twice the speed:
    • The ranking no longer goes through `RankInfo` and `rank_and_count_ties`. `kruskal_test` wants per-group rank sums, not the ranks themselves, so it sorts a 16-byte `(value, group)` pair and adds each tie block's averaged rank straight into the group sums, instead of storing an `NV` rank per observation for a second pass to read.
    • It also no longer calls `qsort`. glibc's `qsort` is a mergesort that allocates a scratch buffer the size of the whole array — measured on glibc 2.39 as a `VmHWM` of 116 MB going to 230 MB across one sort of a 114 MB array — which was half of the function's peak memory, and its comparison goes through a function pointer that cannot be inlined. In its place is a median-of-three introsort that recurses on the smaller partition and loops on the larger, so the stack stays `O(log n)`, with a heapsort fallback past a depth of `2*floor(log2(n))` so an adversarial input cannot drive it to `O(n^2)`, and insertion sort for short runs. Sorting the same five-million-element array takes 0.375s against `qsort`'s 1.01s and allocates nothing. The third change is the group-label array on the `x`/`g` path, which was sized at one pointer per *observation* to hold one per *group* — 40 MB at `n = 5e6` to hold three pointers — and now grows on demand.
    • At `n = 5,000,000` over three groups, measured as `VmHWM` either side of the call:
    • | | 0.3 | 0.301 | |---|---|---| | peak memory | 228 MB (47.8 B/obs) | 76 MB (15.9 B/obs) | | `kruskal_test(\@x, \@g)` | 1.20 s | 0.557 s | | `kruskal_test(\%h)` | 1.11 s | 0.467 s |
    • Sorted, reversed, all-equal, organ-pipe and median-of-three-killer inputs all stay under 0.32s at `n = 2e6`, which is what the depth limit is there for. The sort is checked against an independent pure-Perl implementation of the whole test over 748 structured cases — those shapes at every n either side of the insertion-sort threshold — and 49,712 random ones.
    • The same input now gives the same answer:
    • `H` moved by up to `1.2e-14` between runs on identical data. Nothing was random: the sum of `R_i^2 / n_i` walked the groups by group id, and on the hash-of-arrays path an id is minted in `hv_iternext` order, which is perl's per-process hash order. Equal values were also left in whatever relative order the sort happened to leave them, which came from the same place.
    • The sort now orders by value and then by group, which makes it a total order, and the `k` terms of the sum are ordered before they are added — smallest first, which is the better-conditioned direction as well as a canonical one. `k` is the number of groups, not the number of observations, so it costs nothing next to the ranking. `H` is now bit-identical to R on all 37 corpus cases in all four call forms, and stays so across 60 runs under `PERL_PERTURB_KEYS=1`.
  • Documentation
    • `kruskal_test` gains two sections: what happens to non-numeric, undefined, `NaN` and infinite elements and to a group left with no data, and what the returned fields are — `statistic`, `parameter`, `method` and the p-value under both `p_value` and `p.value` from R's `htest`, plus the `size` and `mean` sub-hashes of `group_stats`, which are computed over the same observations the statistic used.

Documentation

CONTRIBUTING.md
README.md

Modules

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

Other files