| English | 日本語 |
Every exported function is tested against R 4.6.2 (erf and erfc, which R lacks, against
mpmath). The expected values are generated by R scripts
in oracle/ and stored in test/fixtures/ (one JSON line per
case), so the test suite runs without R. Test inputs always include typical values, far tails
(p = 1e-300), boundaries (0, 1, df → ∞), NaN and Infinity, negative values and huge parameters.
substitute in the test file)
and prints the reason with the count.Cases where rStats.js deliberately differs from R are excluded in the test file
(knownDifference) with a written reason and listed in
differences-from-r.md. The test output reports, for each fixture, how
many cases matched R, matched the reference, passed by backward error and were excluded.
The tolerance is a bound, not the aim. Where R is right, rStats.js aims at R’s own bits: a different computation is used only where the result of the usual one is not finite, out of range or has lost its digits, so results in the ordinary range are not changed by a fix for the far range. Where R is wrong, rStats.js returns the right value, only in the range where R is wrong, and the difference is recorded in differences-from-r.md.
Details of the comparison:
test/fixtures/<name>.ref.json (npm run refresh-references). Each oracle registers its
reference function with write_fixture(..., ref = function(row) ...).
Four pnbinom references were computed with mpmath instead (60 digits, betainc), because
the Rmpfr reference is wrong or missing there: pnbinom(1e8, 2.5, mu = 1e6) in both tails
(Rmpfr forms p = size/(size + mu) with 256 bits and then takes 1 − (other tail), which leaves
noise for a result of 8e-106), pnbinom(300, 50, 0.3, lower.tail = FALSE) and
pnbinom(30, 3, 1 - 1e-10, log.p = TRUE). Rerunning npm run refresh-references for
nbinom would bring back the wrong Rmpfr value for the first one.
The incomplete beta function with both parameters large (pbeta and pbinom with a and b from
15 to 1e10, around the region of the uniform asymptotic expansion, both tails and the log scale)
is also checked against mpmath (scripts/mpmath/largebeta.py, audit_pbeta_large and
audit_pbinom_large, tolerance 1e-14; largest errors 9.7e-16 and 3.1e-15). mpmath’s betainc
takes minutes there, so the references are the continued fraction DLMF 8.17.22 at 50 and 80
digits (kept only where the two agree to 1e-28), cross-checked against numerical integration.
With one parameter large and the other small (pbeta with a from 15 to 1e10 and b from 0.1 to 40,
and the mirror image; pt with df from 30 to 1e12, around the region of the expansion for large a)
the references are the positive-term series in the frame of the small parameter at 60 and 90 digits
(scripts/mpmath/smallb.py, audit_pbeta_smallb and audit_pt_smallb, tolerance 1e-14; largest
errors 5.4e-15 and 1.2e-15, from the series of the incomplete beta function; the cases that take
the expansion are within 1.2e-15). The R fixtures of that region (oracle/largea.R) have an
Rmpfr reference for every case.
Five cases where R is off by about 1.0e-14 and the expansion for large a made rStats.js more
accurate (within 4e-16) have mpmath references (100 and 200 digits, betainc, from the double values
of the arguments; scripts/mpmath/mismatch_ref.py), since R was not at hand to compute Rmpfr ones:
pbeta(0.1, 2.5, 1000) (upper tail, and the log of the lower), pnbinom(30, 1000, mu = 2, log.p =
TRUE) and pf(3, 10, 100, ncp = 0) (upper tail, and the log of the lower).
Non-central distributions use Rmpfr sums of the Poisson mixture series or numerical integrals
(oracle/_noncentral.R); regions where no reference can be computed in practical time
(non-central t with df = 1000 or 1e5, non-central F with df2 = 1e9 and x > 30) are left out of
the grid.reference field. Only means, variances and correlation coefficients need high precision, so
this is fast. For the t-test, variances and standard errors are computed with Rmpfr (256 bits)
while the means are first rounded to double, because both R and rStats.js use the difference
of the rounded means reported in estimate; the p-value and quantiles then come from R’s
pt/qt. Fisher’s r × c test enumerates all tables and sums them exactly with gmp. Kendall’s
exact test counts inversions with Rmpfr (1024 bits).choose, lchoose, factorial, lfactorial,
digamma, trigamma, psigamma, the Bessel functions, dmultinom, and gammainc,
erfinv, erfcinv and perm, which R lacks and whose R fixtures are computed from
pgamma × gamma, pracma’s definitions and Rmpfr): besides the R
fixtures, they are checked against mpmath (scripts/mpmath/*.py, 60 to 800 digits,
test/fixtures/<name>.mpmath.json), with the same tolerance of 1e-14 (1e-15 for erfinv,
erfcinv and perm); that comparison uses
rCompatible: false, since the mpmath values are the mathematical function. Where R is
inaccurate or fails, the references in <name>.ref.json were computed with mpmath (100 to 800
digits, from the double values of the arguments) rather than Rmpfr. Where terms cancel
(digamma and psigamma for x < 0, J and Y where they oscillate, I for negative non-integer
orders) the error is measured relative to the size of the cancelling terms, recorded as scale
in the fixture. psigamma for x < 0 is allowed 2e-14: the reflection term is a polynomial of
degree deriv + 1 in cot(pi x), which multiplies the error of cot by up to 101 (also in the R
comparison). In the R comparison of besselJ and besselY, a case whose result is much
smaller than the cancelling terms passes if it is within 1e-14 of the mpmath value relative to
scale (knownDifference); R is closer in relative terms there.log1pmx, lgamma1p, log1mexp, log1pexp, which R has only in its C
API): the fixtures hold correctly rounded values computed with Rmpfr (512 bits, more for
lgamma1p at tiny x so that 1 + x is exact), not R’s output, and are compared with a tolerance
of 1e-14 (the largest observed error is 1.7e-15, for lgamma1p).lmer): lme4’s estimates are where its optimiser stopped, typically 1e-9 to 1e-5 (relative)
from the optimum in θ, so the fixtures carry a reference at the exact optimum: starting from lme4’s θ, the
oracles (oracle/lmer.R, lmermore.R, lmerboundary.R) solve for a zero gradient of the criterion by Newton’s
method in 128-bit arithmetic (Rmpfr; the criterion from the paper’s formulas with dense matrices, the gradient by
central differences), fixing at 0 the elements whose optimum is on the boundary after checking that the
one-sided derivative there is positive, and compute every result field at that θ. The profiled deviance at a
given θ (lmerDeviance) has 200-bit references. The structure (lmerTerms: the model matrices, Z and the layout
of the covariance factor) is compared with lme4’s exactly.solve and determinants have exact references computed with gmp rationals
(every double is exactly representable as a rational). Decompositions are instead checked by
backward error and orthogonality (see below).| Area | Relative tolerance |
|---|---|
| Special functions, densities, distribution functions | 1e-14 (non-central t: 1e-13; ptukey: 2e-14) |
| Distributions outside base R (inverse gamma, Kumaraswamy, Pareto, triangular, Laplace, arcsine) | 1e-14 for d/p/q, against extraDistr (arcsine: R’s beta(1/2, 1/2) moved to (a, b)) and mpmath |
| Quantile functions | 1e-14, or backward error no larger than R’s; discrete quantiles exactly equal or verified by definition |
| Random numbers | discrete values, sample and r2dtable exactly equal, and with L’Ecuyer-CMRG runif and the state after each case (.Random.seed); continuous values within 1e-12; the number of uniform draws consumed is checked exactly |
| Descriptive statistics | 1e-14 (skewness and kurtosis near 0, which come from cancellation, measured absolutely) |
| Hypothesis tests | 1e-12 for every numeric field (cor: 1e-15, pAdjust: exact; roots found by uniroot with tolerance 1e-4, such as the normal-approximation interval of ansariTest: 1e-4); text fields exactly equal |
Linear models (lm, confint, predict) |
1e-12 for every numeric field |
Mixed models (lmer) |
against the reference at the exact optimum: logLik, AIC, BIC, deviance, REMLcrit 1e-12; θ, the fixed effects, their covariance and standard errors, σ, VarCorr, ranef, fitted values and residuals 1e-9 of the largest absolute value of each (largest observed: 1.2e-12, a flat optimum just inside the boundary); elements of θ that are 0 in the reference exactly 0; against lme4 only that the criterion is not larger than lme4’s. The profiled deviance at a given θ 1e-12 |
| Linear algebra | 1e-14 (crossprod, norm, backsolve, chol, qr, the 2-norm from svd, matrix arithmetic, seq, lu) or 1e-13 (solve, determinant, rcond, svd, kappa, ginv, eigen), measured relative to the norm of the result; prcomp 1e-12; pivots, ranks and signs of determinants exactly equal |
| Numerical methods | 1e-15 (approx), 1e-14 (integrate, optimize, uniroot, barylag, hermiteInterp, ode, and grad, jacobian, hessian with an absolute floor of 1e-12) or 1e-13 (spline, splinefun, splinefunH) |
The tolerances were chosen per area from measurement: the largest observed errors against R or the reference are just below them.
ptukey (2e-14): R’s ptukey is accurate to about 8 digits, so nearly every case is
compared with the reference. The reference is a 160-bit trapezoidal rule, used only where
halving the step changed the value by less than 1e-18. Most errors are below 6e-15; the
largest is 1.7e-14 at ptukey(0.5, 1000, 10) ≈ 1.3e-129, where the rounding of log W
(ε |log W|) dominates.The continuous values cannot match exactly because the last bits of Math.log, Math.exp,
Math.pow and Math.tan differ between V8 and glibc’s libm. The largest observed difference is
1.2e-13, for gamma variates with very small shape (e.g. 0.001), where exp(log(u)/shape) magnifies
the error by 1/shape. The number of uniform draws is checked by drawing runif(1) after each
case: a single different accept/reject decision would show up there.
The other RNG kinds (rngRounding, rngLecuyer, rngStreams) are checked the same way, and also
compare the state afterwards (getRandomSeed() against R’s .Random.seed: the kind code, and all
seven integers for L’Ecuyer-CMRG). L’Ecuyer-CMRG’s uniforms are computed with exact integer
arithmetic in doubles, so runif and everything discrete match R exactly; continuous variates
(rnorm, rexp, rgamma) are compared within 1e-12 as above.
method, statistic names and so on) must match exactly. Inputs for which R
throws an error must throw.ansariTest are roots found by uniroot with
tolerance 1e-4, in R as here, so they are compared with a tolerance of 1e-4. Since the same root
finder is evaluated on the same step function, they agree with R to 3e-15. Where R fails with
“object ‘cint’ not found” (24 cases, see differences-from-r.md), the
result is compared with the oracle’s repetition of the case with exact = FALSE.chisqTestSim, fisherTestSim) must match exactly: each case calls
setSeed(seed) before the test, and the runif(1) drawn afterwards must equal R’s, which checks
the number of uniforms consumed.chisqTest(x, y) sorts the levels by code point where R used the locale’s collation; in the one
fixture case where the orders differ, R’s rows and columns are reordered before comparing.crossprod uses the product of the column norms as the scale.args.cond. The allowed error is scale × max(tolerance, 10 ε κ^k), with k from the known
sensitivity: 1/2 for chol; 1 for solve, determinants, rcond, kappa, qr and ginv;
2 for least-squares coefficients (qrCoef, qrSolve), whose error is proportional to κ²
when the residual is non-zero.gap). Vectors for
repeated values (relative gap below 1e-8) and for zero singular values are not unique and are
not compared.qr with dependent columns: Householder vectors for columns after the numerical rank
(singular values above 1e-13 × the largest) are determined only by rounding and are not
compared. For rank-deficient matrices in LAPACK mode, Q, qty, qy and coefficients are not
compared either.test/linalg-invariants.test.ts.qr and chol match R exactly, including R’s limited column pivoting.info block,
because a different OpenBLAS version can change the last bits.None are known in the tested ranges: every case matches R or a high-precision reference within the tolerance.
Outside them, besselJ and besselY with orders in the thousands and x < ν have relative errors
of up to about 1.5e-14 (for example 1.04e-14 for besselY(3000, 3949), 9.9e-15 for
besselJ(1000, 1500) and 1.5e-14 for besselJ(2000, 2999); measured for 0.16.0 against Rmpfr at
25 points with orders 1100 to 9000), from rounding that accumulates over the many steps of the
recurrence and the continued fractions.
R’s qnbinom does not terminate for 13 inputs (e.g. qnbinom(-1, 0.5, 1e-10, log.p = TRUE)).
oracle/nbinom.R limits each case to 2 seconds and drops those rows from the fixture.
Most kernels are ported from SciPy’s xsf (BSD-3) and Boost.Math (BSL-1.0); linear algebra follows the LAPACK reference implementation (modified BSD), LINPACK and EISPACK (public domain). See THIRD_PARTY_NOTICES.md. The main changes made to match or exceed R:
src/special/loader.ts) and the
arguments of exp(−t) in the exponential, Weibull, logistic and log-normal distributions
(src/internal/ddexp.ts) are formed in double-double; at t = 700 the rounding of t alone would
cost 6e-14.src/special/incbet.ts): with large parameters and
x near the mean, the continued fraction is extremely sensitive to the rounding of x (6e-13 at
a = 1e5, x = 0.999). The coefficients a + b and b − 1 are double-double too, since rounding
a + b loses the digits of a small b (5e-11 at a = 1e8, b = 0.01), and from 2^52 on all of them
are, since steps of 1 would be lost.src/special/igam.ts): for
a < 4 and 1.1 < x < 64 the rounding errors of the upper tail’s continued fraction added up to about 120 eps in
double (pchisq(2.3569217160548437, 1, lower.tail = FALSE) was 45 eps off; R 4.6.2 is 25 eps off there and up to
about 60 eps in that region); the recurrence runs in double-double there, within 1 eps, so the upper tails of
pgamma and pchisq (df < 8) are within a few eps (1.2.0). It costs up to about 8 µs instead of 1 µs per call.src/special/igam.ts): in xsf’s region
(|x − a| / a < 0.3 for 20 < a < 200, within 4.5 standard deviations above) the tails were formed from erfc and
e^(−a η²/2) of a rounded η, up to about 35 eps off (pgamma(256.71, 199, lower.tail = FALSE) 22 eps, R 0.2 eps);
the smaller tail now comes from the double-double deviance as farther out, within about 4 eps (1.2.0).src/special/loader.ts): x^a e^−x / Γ(a + 1) for a < 0.01 is
e^−x e^(a log x − lgamma1p(a)); Loader’s form lost up to about 50 eps there through the cancellation of stirlerr(a)
against log(2πa) (pgamma(100, 1e-100, lower.tail = FALSE) was 51 eps off, R 52 eps) (1.2.0).src/special/incbetTemme.ts): with
ab/(a + b) >= 25 and x within the region where the expansion converges, the incomplete beta
function uses Temme’s expansion, whose exponent is the same double-double deviance as Loader’s
binomial density; in the tails a rearranged form avoids the cancellation between the erfc term
and the correction (see performance.md).src/special/incbetLargeA.ts):
with a >= 15, b <= 40 and x in the bulk, DiDonato and Morris’s expansion in incomplete gamma
functions of u = −(a + (b − 1)/2) log x replaces the continued fraction. Its prefactor comes from
the same double-double deviance as Loader’s binomial density, u is formed in double-double, and
the leading incomplete gamma function comes from a continued fraction evaluated backwards, finished
in double-double (the forward evaluation in igam.ts stops too early for small u), so the cases
that take it are within 1.2e-15 (see performance.md).src/special/lgammaDiff.ts) for a ≪ b, used
by lbeta, prefactors and series complements. Near the zeros of lgamma (1 and 2) it uses
lgamma1p (xsf’s formula has relative error 6e-7 at x = 1 + 1e-10).src/internal/invert.ts): Newton/Halley iteration in u = log x, finished
in x itself when |log x| is large (the step in u then limits the precision of x). A step that
does not move because of rounding counts as convergence. qlnorm applies one double-double
Newton correction to the normal quantile, because the rounding of z is magnified by
sdlog |z|. qgeom uses a closed form. Discrete binary searches stop when the answer exceeds
2^53.pnorm tails split q = hi + lo (hi on a 1/16 grid) to compute exp(−q²/2) accurately.
qnorm with logP below −2048 (beyond Cephes’ tail approximation) adds one Newton step using
pnorm.src/internal/poissonMixture.ts): the Poisson mixture is summed
from the mode outwards with positive terms only, recurrences in double-double, and the exponent
kept separately (mantissa × 2^(500k)), so log answers stay accurate below 1e-300. The
non-central t integrates a positive integrand after a one-sided double-exponential change of
variables, so opposite signs of t and δ cause no cancellation; far tails use an integrated-by-
parts form and then the |t|^(−ν) power law. Where the mixture would have too many terms, the
non-central chi-squared density uses an asymptotic expansion of the Bessel function (the
large-argument one, or the uniform one in the order for large df) and its distribution function
integrates that density; the non-central beta and F use limits of X = Y/(Y + W) with Y or W
nearly constant, expanded in its cumulants, or an integral over Y.src/internal/summary.ts): means, variances and correlations are summed in
double-double, and deviations are rescaled by powers of two before multiplying, so data around
1e±200 neither overflow nor underflow.The cost of this in speed, and how to measure what double-double buys, is described in performance.md.