| English | 日本語 |
rStats.js puts accuracy first: distribution functions stay within 1e-14 of the true value even in the far tails (see accuracy.md). That makes them slower than R’s: about 4 to 30 times for the distribution functions (most d/p functions about 10 times), 3 to 25 times for linear algebra. A single call still takes about a microsecond (a quantile function a few to 15 µs), and a test of a typical sample well under a millisecond.
µs per call, one call at a time. The R column is R 4.5.2 on the same machine (one vectorised call over the same kind of inputs, divided by their number).
| Function | rStats.js | R |
|---|---|---|
pbeta |
1.6 | 0.16 |
dbeta |
0.94 | 0.11 |
pt |
1.2 | 0.12 |
pgamma |
0.66 | 0.17 |
pf |
1.5 | |
pchisq |
0.52 | |
dt / dgamma |
0.62 / 0.47 | |
dbinom / pbinom |
0.54 / 1.4 | |
pnbinom |
12.2 |
| Quantile function | rStats.js, typical | rStats.js, tail | R, typical | R, tail |
|---|---|---|---|---|
qbeta |
15.2 | 14.9 | 0.91 | 2.7 |
qt |
8.8 | 9.8 | 0.33 | |
qgamma |
4.8 | 4.0 | ||
qchisq |
4.3 | |||
qf |
11.4 |
“Typical” inputs: shape parameters log-uniform on 0.5 to 50, p uniform on 0 to 1. “Tail”: shapes on 0.1 to 1e4, p from 1e-5 to 1e-300. The results are the same as R’s (see accuracy.md); the gap is the cost of the double-double arithmetic (see the reference section below).
Large parameters take constant time. pbeta(0.5, a, a) takes 1.0 µs for every a from 1000 to 1e6 (2.4 µs at a
= 10), pbinom(n/2, n, 0.5) 1.0 to 1.1 µs for n from 1e3 to 1e6, pt with df = 1e6 1.1 to 1.3 µs, pf with 1e5
denominator degrees of freedom 1.2 to 2.0 µs; qt with df = 1e6 takes 4.0 µs, qf(p, 5, 1e5) 7.5 µs,
qbeta(p, 1e6, 2.5) 17.1 µs, qbinom(p, 1e6, 1e-5) 3.0 µs.
A test: tTest (Welch, 30 and 25 observations) takes 17.3 µs; most of it is qt for the interval.
Matrices of standard normal values (chol: AᵀA + nI; symmetric eigen: (A + Aᵀ)/2), ms per call. R uses OpenBLAS
0.3.26 with 4 threads (blocked LAPACK, SIMD); rStats.js uses plain unblocked loops.
| Operation | n = 100: rStats.js | R | n = 300: rStats.js | R |
|---|---|---|---|---|
qr |
0.76 | 0.25 | 18.4 | 4.9 |
solve (one right-hand side) |
0.48 | 0.10 | 11.0 | 1.0 |
chol |
0.22 | < 0.05 | 4.7 | 0.20 |
svd |
6.5 | 1.3 | 162 | 15.3 |
eigen (symmetric) |
4.0 | 0.95 | 96 | 7.7 |
eigen (general) |
10.7 | 3.5 | 263 | 32.9 |
For the matrix sizes typical in statistics (up to a few dozen columns) a call takes well under a millisecond to a few ms.
One fit of lmer (REML) from the data, structure included, the median of 3 fits; lme4 2.0.6 with its default
optimiser on the same data (written out from scripts/lmer-timing.mjs’s seeded generator), the median of 3. Both
use a sparse Cholesky factorisation; rStats.js’s criterion at its estimate was never larger than lme4’s.
| Model | n | q | rStats.js | lme4 |
|---|---|---|---|---|
y ~ x + (1 \| s) + (1 \| i), 200 subjects × 100 items crossed |
20000 | 300 | 0.59 s | 0.20 s |
| the same, 500 × 200 | 100000 | 700 | 4.13 s | 1.51 s |
y ~ x + (x \| s), 1000 subjects × 10 |
10000 | 2000 | 0.20 s | 0.07 s |
| the same, 5000 subjects × 10 | 50000 | 10000 | 0.87 s | 0.35 s |
With two or more random-effects terms, a fit takes about twice as long as the local minimisation alone (0.32 s and
2.09 s for the crossed models in 0.14.0): after it, the fit is compared with the boundary points where one term is
removed, each minimised again (see “Local minima” in differences-from-r.md). Small models (the
fixtures: up to a few hundred rows) take 0.1 to 7 ms. node scripts/lmer-timing.mjs measures the
fits and single evaluations of the criterion.
lmerSummary (Satterthwaite’s degrees of freedom: the analytic Hessian of the deviance and the Jacobians of the
covariance of the fixed effects) adds, after the fit: 0.03 s (crossed 200 × 100), 0.25 s (500 × 200), 0.01 s
(1000 subjects) and 0.04 s (5000 subjects), Node 24 on arm64.
One fit of glmer (the Laplace approximation unless stated; data from a seeded generator, random intercepts of
standard deviation 0.5 and one uniform predictor), Apple M1 (Mac mini), Node.js 24, the code of 2026-11-05 (after
0.17.0). The deviance’s gradient is analytic, so the cost grows about linearly with the rows; lme4 fits these in
well under a second, since it stops its optimiser and penalised IRLS early (rStats.js’s estimate is the exact
optimum). There is no size limit: 100000 rows take about 13 s.
| Model | n | rStats.js |
|---|---|---|
poisson y ~ x + (1 \| g), 100 groups |
20000 | 2.2 s |
binomial y ~ x + (1 \| g), 200 groups |
20000 | 2.4 s |
binomial y ~ x + (1 \| g), 500 groups, nAGQ = 100 |
5000 | 4.7 s |
poisson y ~ x + (1 \| g) + (1 \| h), 200 × 37 crossed |
4000 | 1.6 s |
poisson y ~ x + (x \| g), 200 groups |
4000 | 1.6 s |
poisson with 9 predictors + (1 \| g), 200 groups |
4000 | 1.5 s |
poisson y ~ x + (1 \| g), 1000 groups |
100000 | 12.9 s |
npm run size (scripts/size.mjs) bundles each import from the built dist/ with esbuild, then minifies and gzips
it. Tree-shaking keeps only the imported functions and their dependencies.
| Import | Minified | gzip |
|---|---|---|
pnorm |
3.1 KB | 1.6 KB |
setSeed, rnorm |
4.0 KB | 2.0 KB |
solve / svd / eigen |
5.6 / 5.2 / 12.7 KB | 2.4 / 2.4 / 5.1 KB |
wilcoxTest |
46.7 KB | 18.9 KB |
fisherTest |
62.8 KB | 27.9 KB |
tTest |
70.9 KB | 31.0 KB |
qbeta |
98.4 KB | 43.9 KB |
pt, qt |
89.7 KB | 40.4 KB |
glm, glmSummary / glm with all its methods |
85.2 / 122.3 KB | 35.0 / 51.4 KB |
everything (except /mixed) |
767.6 KB | 289.5 KB |
@numero-uno-jpn/rstats/distributions (all of it) |
241.7 KB | 95.6 KB |
@numero-uno-jpn/rstats/tests (all of it) |
591.7 KB | 229.0 KB |
@numero-uno-jpn/rstats/linalg (all of it) |
65.6 KB | 24.1 KB |
@numero-uno-jpn/rstats/descriptive (all of it) |
58.8 KB | 21.7 KB |
@numero-uno-jpn/rstats/numerics (all of it) |
64.6 KB | 26.3 KB |
@numero-uno-jpn/rstats/special (all of it) |
67.6 KB | 30.4 KB |
@numero-uno-jpn/rstats/mixed: lmer alone / with lmerTest’s tests |
50.1 / 164.6 KB | 19.4 / 70.0 KB |
@numero-uno-jpn/rstats/mixed: glmer alone / with glmerSummary, glmerAnova, glmerPredict |
116.9 / 160.4 KB | 47.2 / 66.1 KB |
@numero-uno-jpn/rstats/mixed (all of it) |
225.5 KB | 91.8 KB |
The subpath rows are an area bundled on its own. The per-area CDN files (dist/<area>.js) share code between areas
through dist/chunks/, so one area loaded from them transfers more: 39 KB (special) to 250 KB (tests) gzipped,
chunks included. pt is large because the non-central t and the double-double building blocks come with it; the
tests with a formula form bring the formula and model-frame code.
A call that would allocate more than 2^24 values in one array, or run for more than about 10 seconds on the
development machine, throws a RangeError before it allocates or computes. The exact tests name the way out in the
message (use exact: false). R computes most of these however long they take, or runs out of memory. Where an
allocation could be avoided without changing the values or the random stream, it was (rwilcox with m + n above
2^24, the goodness-of-fit simulation of chisqTest).
| Function | Throws when |
|---|---|
fisherTest (2 × 2), mantelhaenTest(exact: true) |
the support of the count has more than 2^24 values; mantelhaenTest also when the total reaches 2^53 |
oddsRatio |
"midp" and "fisher": the support of the count has more than 2^23 values (about 6 s), and for "fisher" an entry is above 2^31 - 1; every method: the hypergeometric p-values would take more than 2^26 terms of phyper (counts of about 1e13) |
boschlooTest |
(n1 + 1)(n2 + 1) is above 2^24, or the grid over the nuisance parameter would take more than 3.4e8 steps (about 8 s; max(1000, 32 sqrt(n1 + n2)) points times min(n + 1, 40 sqrt(n)) for both groups: a small group against one of 270000) |
ksTest(exact: true) (one sample) |
the matrix power needs more than 2^33 multiply-adds |
ansariTest(exact: true), jonckheereTest(exact: true) |
the distribution needs more than 2^24 values or 2^32 steps |
corTest(method: "kendall", exact: true) |
n is above 1290 |
pAdjust(method: "hommel") |
n is above 65536 |
r2dtable; simulated p-values of chisqTest (tables) and fisherTest |
the grand total is above 2^24 |
chisqTest goodness-of-fit simulation |
B times the total count is above 2^27 |
dwilcox, pwilcox, qwilcox (and the exact two-sample wilcoxTest) |
the table would take more than 2^29 steps |
wilcoxTest(exact: true) with ties or zero differences (the conditional test) |
its distribution, built once for the p-value and again for each multiset of scores that the interval’s search meets, would take more than 5e9 steps (about 8 s) or one array would hold more than 2^24 values |
phyper / qhyper |
the sum would take more than 2^26 terms / the standard deviation is above 2^21 |
besselJ, besselY, besselI, besselK |
a series or recurrence would take more than 2^25 to 2^31 steps, or a phase is beyond what doubles resolve |
rwilcox / rsignrank |
n is above 2^24 / n is 2^27 or more |
weighted sample, sampleInt |
the size is 2^31 or more (as R) |
tapply, split, aggregate |
the cross-classification has 2^31 cells or more (R’s limit too) |
expm (Taylor and Padé methods) |
order × (n^3 + 64) is above 2^31 |
grad, jacobian, hessian (Richardson) |
r is above 512 |
mcreg (Passing-Bablok) |
n is above 5793 (2^24 pairwise slopes), or all the fits of the call (the data, the jackknife, the resamples and the nested resamples) would take more than 2^27 slopes (about 10 s; 999 resamples of 518 points) |
The details, with timings and R’s behaviour, are in differences from R (Hypothesis tests
and Random numbers under “Deliberate differences”). There is no bound on the length of an output the
caller asks for (seq, approx, spline, sampleInt, the r* functions) or on iteration counts such as optim’s
maxit: their time grows in proportion.
Apple M1 (Mac mini), Node.js 24, the tsc build, 2026-10-26 (version 0.11.0). Distributions and tTest: npm run
bench (scripts/bench.mjs), one call at a time, the minimum of 5 runs; its inputs come from the library’s own
setSeed/runif, so every build sees the same ones. Linear algebra: the minimum of 5 runs of the mean of 20 calls
(3 for n = 300). R 4.5.2 (arm64) in the same machine’s container, the same kind of inputs from R’s own random
numbers. The bundle sizes were measured on 2026-11-05 with the code of 1.0.0 (Node.js 24; the sizes do not
depend on the machine). The mixed models were measured on 2026-10-27 (the code of 0.12.0; the fits of lmer and
lmerSummary again with the code of 0.15.0), with lme4 2.0.6 in the same container.
Background for changing the code; not needed for using the library.
Why it is slower. Where double precision cannot reach 1e-14, the computation is carried in double-double
arithmetic (an unevaluated sum of two doubles, about 32 digits; src/special/dd.ts, src/internal/ddexp.ts): the
exponents of tiny probabilities (Loader’s saddle-point deviance bd0), the pair x, 1 − x and the continued fraction
of the incomplete beta function, the Poisson-mixture sums of the non-central distributions, the sums of squares of
the tests, and the Wilcoxon and studentized-range integrals. Where an error bound shows that double precision is
enough, the double-double steps are left out (the tail of Loader’s series when |x v³| <= 0.12, one logarithm in the
far branch of bd0, dbeta on the linear scale). The continued fraction itself stays in double-double: in plain
double its error reaches 7e-15 for a + b below 10 and 1.5e-13 around 1000, from the rounding of its coefficients.
What keeps it fast.
src/special/incbetTemme.ts) instead of a continued fraction of length
√(a + b); with one parameter large and the other small (pt, pf with large degrees of freedom, pbinom near 0
or n), DiDonato and Morris’s expansion (TOMS 708, eq. 9; src/special/incbetLargeA.ts). Both were checked
against mpmath (scripts/mpmath/largebeta.py, smallb.py) and stay within about 1e-15.[hi, lo] arrays (python3 scripts/ddgen.py ops |
bd0dd | ddlog | scalar | cf), with the same operations in the same order, so the results are bit-identical.ddLog uses a 128-point table and a degree-5 polynomial (maximum relative error 4.4e-32).src/internal/invert.ts): Halley steps in log x with the
closed-form derivative of the log density, one distribution-function call per step giving both scales, and
starting values from Abramowitz and Stegun (qbeta, qf), Cornish–Fisher (qt) and xsf (qgamma); a typical
qbeta takes about 3.5 iterations.Double-precision mode. An internal switch (setDoubleDouble in src/special/dd.ts, not public) replaces every
double-double building block with plain double arithmetic. npm run bench times both modes (in separate processes;
the default mode is now 1.0 to 2.3 times slower than the double mode, and faster for some quantiles whose iterations
converge worse in double precision), and npm run precision-report runs the test suite in both and tabulates the
errors (reports/precision/report.md, not committed). In double precision about half of the fixtures still meet
1e-14; most main d/p functions stay within 1e-10, and some get far worse (beta 3e-5, pt 6e-7, the discrete
quantiles, Pearson’s correlation test near ±1).
Options not taken.
ddLog compiled to Wasm was only 1.2 times faster (1.6 with a relaxed-SIMD fused multiply-add whose
fusing depends on the engine), at the cost of asynchronous loading, no per-function tree-shaking and two code paths.stirlerr for 1 <= n <= 15 (absolute error 8e-16 → 1.8e-16): it moves two fixture cases where R
itself is 1e-14 off outside the tolerance, and hardly changes the speed.lbeta, lgamma, stirlerr): 1 to 6% of the time, no
measurable gain.