Introduction
When performing statistical analysis in JavaScript, one of the primary concerns is: "Can the calculation results truly be trusted?" In scenarios such as the tails of probability distributions, cases involving extremely large parameters, or instances with infinitesimal probabilities, the choice of numerical computation method can lead to significant discrepancies in the results.
To address this challenge directly, I am developing rStats.js, a statistical computation library for JavaScript. In this project, I use R as a key benchmark to continuously verify the accuracy of these statistical calculations.
In this article, I will explain the methodology used to compare rStats.js against R, as well as the level of precision I aim to achieve.
1. Automated Verification Against R 4.6.2
In rStats.js, the general policy is to compare every public-facing function against R version 4.6.2. This goes beyond simply testing a few representative input values; the test suite incorporates a diverse range of input conditions, including:
- Typical input values
- Extremely small probabilities (e.g., 1e-300)
- The extreme tails of distributions
- Boundary values (e.g., 0, 1, or cases where degrees of freedom approach infinity)
- NaN and Infinity
- Negative values
- Extremely large parameters
The expected values—calculated using R—are generated via dedicated R scripts and stored as JSON fixtures. Consequently, there is no need to have R installed when running the tests; the results from R are pre-saved, allowing the JavaScript test suite to verify against them.
2. Precision Standards Vary by Function Type
When evaluating the accuracy of numerical computations, applying a uniform margin of error across all functions is not appropriate. In rStats.js, the allowable margin of error is tailored to the specific nature of each calculation.
| Calculation Type | Relative Tolerance |
|---|---|
| Special functions, probability densities, and distribution functions | 1e-14 |
| Descriptive statistics | 1e-14 |
| Hypothesis testing | 1e-12 |
| Linear models | 1e-12 |
| Linear algebra | 1e-14 – 1e-13 |
| Numerical methods | 1e-14 – 1e-15 |
| Mixed-effects models | Primarily 1e-12 – 1e-9 |
These values are not merely design targets; they are tolerances established based on the results of running a large number of test cases.
Note that not all results can be evaluated solely by simple relative error. In cases where values are close to zero or where numerical cancellation occurs, evaluation is performed using absolute error or an appropriate scale.
3. R’s Results Are Not Always Correct
This is an aspect to which we pay particular attention during the development of rStats.js.
While R is an excellent environment for statistical computing, its numerical calculations are not necessarily accurate under all conditions. In scenarios involving extremely small probabilities or exceptionally large parameter values, even R can be subject to rounding errors or loss of precision.
Therefore, if our results do not match those of R, we do not automatically assume that the error lies with rStats.js. Instead, we calculate high-precision reference values when necessary to verify which result is mathematically correct.
To calculate these reference values, we primarily use the following methods:
- Arbitrary-precision arithmetic (256–1300 bits) via Rmpfr
- Exact rational number arithmetic via gmp
- High-precision computation via mpmath
- Independent verification using numerical integration or mathematically derived series
For instance, if the value of a probability distribution function does not match R’s output, but it does align with our high-precision reference value, we consider that result to be "more accurate than R." ## 4. Cases where more accurate results are returned than in R
Here are a few specific examples.
4-1. Extremely small probabilities
Consider the following calculation:
pgamma(1e-100, 3, lower.tail = FALSE, log.p = TRUE)
- R result:
-1.6666666666665175e-301 - High-precision reference value:
-1.6666666666666667e-301 - rStats.js:
-1.666666666666667e-301
In this case, the result from rStats.js is much closer to the high-precision reference value.
4-2. Non-central t-distribution calculations
For the non-central t-distribution, there are cases where R fails to correctly calculate extremely small probabilities. For example:
pt(1, 5, ncp = 40)
- R:
1.009e-303 - High-precision reference value:
1.0887e-287 - rStats.js: Matches the high-precision reference value
As shown here, there are instances where simply replicating R's output is insufficient.
Note that the values presented here are taken from the verification examples provided in docs/differences-from-r.md on GitHub.
5. Evaluating quantile functions using a different approach
Quantile functions for probability distributions require a different evaluation method than standard functions. In regions where the cumulative distribution function is very flat, even a tiny difference in probability can lead to a significant discrepancy in the quantile value. Therefore, simply comparing the output against R's results using relative error is not an appropriate approach. rStats.js employs a "backward error" approach, where the result of a quantization function is fed back into the corresponding distribution function for evaluation. In other words, it assesses how closely the value returned by rStats.js—when substituted back into the distribution function—approximates the target probability.
Even if a result differs from that of R, it is considered acceptable provided the backward error is equal to or less than that of R. For discrete distributions, it is also verified that the function returns the smallest integer that satisfies the definition.
6. Verification of Random Number Generation
In addition to statistical functions, rStats.js implements random number generation. It includes a setSeed() function (analogous to R’s set.seed()) and supports various random number generation algorithms, such as Mersenne-Twister and L'Ecuyer-CMRG.
For random number generation, the following aspects are verified:
- Discrete random numbers and sampling results match those of R.
- The state of the L'Ecuyer-CMRG generator matches that of R.
- Continuous random number values fall within an acceptable margin of error.
- The total count of random numbers consumed matches that of R.
For discrete results, an exact match is required, while for continuous results, the standard is generally a deviation of no more than 1e-12.
Note that for continuous random numbers, results may not match down to the least significant bit; this is due to differences between JavaScript’s Math functions and the mathematical libraries used by the R runtime environment. Such discrepancies are intentional and are considered to be within an acceptable range.
7. Documenting Cases Where Results Differ from R
The goal of rStats.js is not necessarily to achieve an exact match with R in every instance. Because the priority is to ensure mathematically and statistically correct results, the library may occasionally produce outputs that differ from those of R.
Rather than silently excluding such cases from the test suite, these instances are documented as knownDifferences—with the underlying reasons explained—and are made publicly available in the docs/differences-from-r.md file on GitHub. In some cases—such as when R returns NaN or Infinity for undefined results—rStats.js is designed to throw an error instead. By explicitly documenting these discrepancies, we ensure that users are aware of how the library’s behavior differs from that of R.
8. Implementation Supporting Numerical Precision
To achieve this level of precision, rStats.js employs not only standard double-precision floating-point arithmetic but also "Double-Double" arithmetic where necessary. Double-Double arithmetic combines two double-precision floating-point numbers to provide a greater number of significant digits than is otherwise possible.
For instance, this technique is used in calculations such as:
- Computing the exponent part of exponential functions.
- Evaluating continued fractions for incomplete beta functions.
- Evaluating continued fractions for incomplete gamma functions.
- Asymptotic expansions for large parameters.
- Operations prone to catastrophic cancellation (loss of precision).
Rather than applying high-precision arithmetic to every calculation, we limit its use to areas where precision is a critical concern, thereby minimizing the impact on the performance of standard calculations.
9. Transparency Regarding Validation Results and Limitations
rStats.js provides full transparency regarding its validation methods and its differences from R via GitHub:
Within the scope of our validation, no known precision issues have been identified. However, this does not constitute a guarantee of mathematical perfection for all possible input values. Certain precision limitations remain outside the validated range—for example, when dealing with Bessel functions of extremely high order.
We believe that what matters most is not merely claiming high precision, but being transparent about the validation methods used and the inherent limitations of the implementation. ## Conclusion
I believe that JavaScript is not merely a language for developing web applications; it can also serve as a platform for performing reliable statistical computations.
In rStats.js, we verify the accuracy of numerical calculations by using R as a key benchmark, alongside high-precision reference values. Moving forward, I intend to not only expand the range of available statistical methods but also further enhance the reliability and verifiability of these numerical computations.
rStats.js is an open-source project licensed under the MIT License.
GitHub: https://github.com/numero-uno-jpn/rStats.js/
I welcome feedback from developers interested in statistical computing, numerical analysis, and data analysis using JavaScript.
Top comments (0)