Changelog

All notable changes to this project will be documented in this file.

The format is based on Keep a Changelog, and this project adheres to Semantic Versioning.

[Unreleased]

Added

scripts/precision-refs-test.py (generates test/precision-test.js, checked via scripts/eval-test.js): a reference-value precision gate for the 9 src/test/ hypothesis tests (bartlett, levene, brownForsythe, mannWhitney, welch, hsic, andersonDarling, cramerVonMises, kolmogorovSmirnov), closing the gap noted in #1270 — every other numeric module already had one, src/test/ had none. Unlike the mpmath-sourced gates, references come from R 4.3.3 (base R + goftest 1.2.3), since these are canonical statistical procedures rather than bare mathematical functions; hsic has no R equivalent (R's dHSIC defaults to a permutation test, not the Gamma-approximation variant ranjs implements) and is instead an independent from-first-principles Python reimplementation of Gretton et al. (2008). levene/brownForsythe reference the identical textbook one-way-ANOVA-on-absolute-deviations formula car::leveneTest computes internally, via base R only (car's lme4/quantreg/RcppEigen dependency chain failed to link natively in the generating environment). Three documented convention gaps between ranjs and R's defaults are reconciled rather than silently worked around: mannWhitney has no tie-variance or continuity correction (reference datasets are drawn tie-free, R forced to correct=FALSE, exact=FALSE); kolmogorovSmirnov is always asymptotic (R forced to exact=FALSE); cramerVonMises omits R's finite-sample correction (referenced via goftest::pCvM(stat)'s default n=Inf, not the top-level cvm.test() wrapper — the same public-wrapper-bundles-an-extra-correction failure mode previously hit with scipy's cramervonmises()). Two further numeric gaps were found and deliberately left unfixed per this issue's explicit scope (fixing disagreements is out of scope; flagged for separate follow-up): andersonDarling's p-value, which already includes ranjs's own Marsaglia & Marsaglia (2004) finite-sample correction (_errfix), differs from R's independent implementation of the same correction by up to ~7e-5 relative; ran.dist.Kolmogorov's survival function differs from R's own Kolmogorov-distribution CDF by up to ~1e-4 relative at some arguments (isolated to the distribution itself, not the kolmogorovSmirnov test, by evaluating both at an identical statistic).
scripts/difftest-special.py (npm run difftest:special): a differential-testing harness for src/special/ that sweeps randomized, seeded input grids (10,000 points per function by default) far denser than the committed precision gate, evaluating both mpmath (mp.dps=50, live) and ranjs at every point via the existing scripts/eval-special.js bridge, and reporting per-function max/median/p99 ULP error plus the worst-case reproducer as JSON. Unlike test/precision-special.js, this harness commits no reference literals and runs entirely out-of-band from npm test — it is a non-blocking diagnostic/audit layer, not a merge-blocking gate; see ADR-0052. Initial coverage is the six functions already bridged by eval-special.js (besselI, besselISpherical, besselInu, besselK, besselKnu, digamma); broader coverage is a follow-up (#1271). Confirms the known besselK/besselKnu series/asymptotic crossover degradation near x=6 (_X_K_SERIES, already accepted and tolerance-documented in scripts/precision-refs-special.py) as elevated max-ULP, and additionally surfaces the same crossover at larger besselKnu order (nu up to ~5) than the fixed grid probes (which stops at nu=3.3) — a defect the fixed grid cannot reach, left unfixed per this issue's scope (#1264). #1271's first follow-up extends this same bridge/gate/harness trio to the gamma/beta cluster — gamma, logGamma, gammaLowerIncomplete, gammaUpperIncomplete, gammaLowerIncompleteInv, beta, logBeta, betaIncomplete, regularizedBetaIncomplete, logBinomial — with the same threshold-focused grid philosophy (each new precision-refs-special.py grid builder straddles a documented internal dispatch threshold, e.g. gammaLowerIncomplete/gammaUpperIncomplete's x=s+1 series/continued-fraction crossover and the shared _deviance.js stirlerr/bd0 thresholds) and calibrated (never blind) ulp_ceiling values in difftest-special.py; the zeta/polylog cluster (#1414) and the remainder cluster (#1415) are further follow-ups. Surfaced two pre-existing defects in the newly-swept functions, both fixed in this same change (see Fixed below): gamma()'s premature overflow and beta()'s sign loss for negative non-integer arguments.
scripts/difftest-dist.py (npm run difftest:dist): extends the differential-testing harness to distribution pdf/cdf sweeps, piloted on the gamma/beta family (Gamma, Beta, Chi2, F, StudentT, InverseGamma), which share the gammaLowerIncomplete/gammaUpperIncomplete and regularizedBetaIncomplete composition chains #1264 already sweeps directly. A new scripts/eval-dist.js bridge constructs a distribution from a name and parameter tuple and evaluates pdf/cdf at a point derived from the distribution's own quantile method (q(p) for a randomly-drawn p ~ Uniform(0.001, 0.999)), guaranteeing every probe lands in the distribution's interior without the cost of an mpmath bisection per point. 10,000 randomly-drawn parameter tuples per distribution (60,000 total) surface two real defects the curated 3×5-point precision gate misses: Gamma.pdf returns NaN for large alpha (e.g. alpha≈96.5) because Math.pow(x, alpha-1) overflows to Infinity while Math.exp(...) underflows to 0, giving 0 * Infinity = NaN; InverseGamma.pdf silently returns 0 instead of a tiny representable density for extreme quantile draws because x*x overflows to Infinity in super._pdf(1/x) / (x*x). Both are left unfixed per this issue's explicit scope (#1265) and flagged for follow-up. A _formula_self_check(), run unconditionally from main() alongside the existing ulp_diff self-check, validates the six hand-derived pdf/cdf reference formulas themselves against exact closed-form identities (e.g. Gamma(1,1) ≡ Exponential(1), StudentT(nu=1) ≡ standard Cauchy, the DLMF 8.17 F/beta CDF relation) so a future transcription error in a reference formula is caught immediately rather than silently corrupting every sweep result. scripts/eval-dist.js restricts distribution construction to an explicit whitelist matching scripts/difftest-dist.py's DIST_SPEC keys (Gamma, Beta, Chi2, F, StudentT, InverseGamma), the same pattern scripts/eval-special.js's FN map already uses, so a JSON input's dist field can no longer instantiate an arbitrary exported distribution by name — only the six the harness actually drives.
scripts/difftest-quantile.py (npm run difftest:quantile): extends the differential-testing harness to quantile accuracy, its own dimension per #1269 with a different error mechanism (numerical CDF inversion) than pdf/cdf. Two metrics, both live-mpmath and out-of-band from npm test (ADR-0052): a round-trip sweep (|cdf(q(p)) - p|) over every one of the ~146 distributions in test/dist-cases-*.js, needing no external reference since cdf and q are both ranjs's own methods, for p drawn log-uniformly toward both tails ([1e-6, 1e-1] ∪ [1-1e-1, 1-1e-6] — not 1e-12: several discrete _cdf are O(k) recurrence sums whose double-precision summation saturates before reaching that close to 1, and probing past the saturation point sends the base class's bracket-expansion search hunting for a k that doesn't exist, hanging indefinitely); and an absolute ULP accuracy sweep against an independently mpmath-derived inverse-CDF reference (bracket expansion + bisection in log-space/logit-space for domains bounded at 0 or to (0,1), since an additive step is meaningless once the true root is hundreds of orders of magnitude from the seed) for the #1265 pilot family. A new scripts/eval-quantile.js bridge adds a catalog mode reading each distribution's canonical parameter tuple, type, and closed-form-vs-numerical quantile status (typeof instance._q === 'function', only readable from JS) straight from test/dist-cases-*.js, so the swept population and the closed-form/numerical split can never drift from what's actually tested. Non-convergence (NaN), out-of-support returns, and non-monotonicity are reported as hard failures distinct from the error statistics. The sweep surfaces several real quantile defects, left unfixed per this issue's explicit scope and flagged for follow-up: LogCauchy's closed-form _q underflows to 0 / overflows to Infinity well before the true quantile is unrepresentable; StudentT's closed-form _q is non-monotonic and returns a wrong-signed value in the extreme lower tail for small nu; Beta's numerically-inverted _q underflows to its lower boundary for extreme shape parameters at small p; Gamma/InverseGamma's closed-form _q returns NaN (non-convergence) for extreme shape parameters near p=0 or p=1.
docs/accuracy.md: a committed, documented-accuracy-bounds table generated by the new scripts/generate-accuracy-docs.js (npm run accuracy, chaining accuracy:special → accuracy:dist → accuracy:docs) from the #1264/#1265 differential-testing harness JSON reports, closing the "Documented accuracy bounds" gap todo.md tracked under Publication-Grade Gaps. Every special function src/special/index.js exports and every distribution src/dist/index.js exports is listed — swept ones with their measured domain (read straight from the harness report's own domain field, so the table can never drift from what was actually sampled), max/median ULP, and sample count; unswept ones as an explicit "not yet measured" row rather than a silent omission. Gamma.pdf/InverseGamma.pdf (#1265's NaN/overflow defects, now #1363/#1364) and besselK/besselKnu's x=6 series/asymptotic crossover (#1140) render with their actual measured (bad) values — ∞ for the NaN-mismatch divergences, ULP counts in the billions for the Bessel crossover — each linked to its tracking issue, never rounded away. docs/accuracy.md is committed rather than generated at docs-build time, so it is readable on GitHub without a Python + mpmath environment and its diffs are reviewable per-PR; see ADR-0053. Both harness scripts' build_report() gained a domain/mp_dps field (read from the same SWEEP_SPEC/DIST_SPEC dict the sweep itself draws from) to support this without risking the "for |x| <= Y" claim drifting from what was actually measured. generate-accuracy-docs.js's statusFor() now composes every applicable flag instead of stopping at the first match: a divergence count no longer disappears once an entry also carries a KNOWN_ISSUES link (Gamma.pdf's row now reads "31 divergence(s) ... known accuracy gap" instead of swallowing the count), and a report's errors field (points where the harness's Node eval bridge threw rather than returning a value — a distinct, more severe failure mode than a returned-but-wrong value) is now rendered too, previously tracked in the JSON report but never surfaced in the table at all. Each of the four flags (thrown error, divergence, ceiling breach, known-issue link) now gets an emoji chosen for what it means rather than a shared ⚠️ or an arbitrary severity color — 💥 the eval crashed, ❌ a value came back but is NaN/nonsensical against mpmath, ⚠️ a real finite value worse than its calibrated ceiling, 🔗 a pointer to an already-tracked non-new problem, ✅ OK — so the failure mode reads at a glance without parsing the status text.
.github/workflows/difftest.yml: runs the #1264/#1265 differential-testing harness on a weekly schedule (plus workflow_dispatch for manual runs), separately from ci.yml — the harness needs a Python + mpmath environment and sweeps far denser grids than test/precision-*.js, so it stays out-of-band from the fast, merge-blocking unit-test gate (ADR-0052). The job runs accuracy:special/accuracy:dist, uploads both JSON reports as a workflow artifact, then runs the new scripts/difftest-ci-gate.js (npm run difftest:ci-gate), which fails the job when any function/distribution exceeds its declared ulp_ceiling and writes the exact reproducer — function/distribution, parameter tuple, evaluation point, ranjs value, mpmath value, ULP distance — to the job summary. Deliberately does not auto-file or update a tracking issue on failure: a red run plus the uploaded report is a sufficient signal, and a per-run auto-filed issue would duplicate weekly on top of what stale.yml already manages. The gate is a separate script rather than added to the harness scripts themselves, since npm run difftest:* also doubles as a plain local diagnostic run that should not start failing the shell over an already-tracked, already-calibrated defect (e.g. besselK's #1140 crossover). difftest-ci-gate.js now also fails the job on divergences > 0 or errors > 0, not only ceiling_exceeded — an inf ULP distance (the harness's encoding for a NaN/Infinity mismatch against a finite mpmath reference) is deliberately excluded from ceiling_exceeded's comparison, so without this a regression to NaN/Infinity on a previously-clean function/distribution would have stayed invisible to the gate (#1369). A KNOWN_ISSUES allowlist (mirroring generate-accuracy-docs.js's own map) keeps the two already-tracked divergence sources, Gamma.pdf (#1363) and InverseGamma.pdf (#1364), from turning the job permanently red on ship; allowlisted entries still appear in the job summary table, tagged with their tracking issue, rather than silently disappearing. The summary table gained a Reason column so a divergence or eval-error failure reads distinctly from a plain ceiling breach instead of leaving a reviewer to guess why a row with a small max_ulp still failed.

Changed

Code Health of test/process.js improved from 9.09 to 10.0 by splitting it into a test/process/ directory, one file per process (mirroring the existing test/mc/ per-sampler layout), plus a shared _helpers.js and a reference-values.js consuming test/process-cases.js. No behavior change — same test cases, same assertions.
Code Health of test/special.js improved from 9.09 to 10.0 by splitting it into a test/special/ directory, one file per src/special/ module (mirroring the existing test/mc//test/process/ per-module layout; hurwitz-zeta.js/riemann-zeta.js share a single zeta.js test file since one existing test asserts an identity across both), and extracting a checkReferenceValues helper in error.js to eliminate duplication between the erf/erfc reference-value tests. No behavior change — same test cases, same assertions.
Code Health of test/test-utils.js improved from 9.17 to 10.0 by extracting the discrete/continuous branches of runX and Tests.cdf2pdf into named helpers (runXDiscrete/runXContinuous, cdf2pdfDiscrete/cdf2pdfContinuous/cdf2pdfContinuousAt), and extracting chiTest's frequency-map construction and chi-square binning into frequencyMap/binChiSquare, removing a dead no-op if block along the way. No behavior change — same test cases, same assertions (verified by mutation-testing a deliberately injected Normal._pdf bug against the refactored helpers before reverting it).
Code Health of src/special/owen-t.js improved from 9.09 to 9.38: findSectorRow/findSectorColumn's duplicated linear-search loops merged into a single findSector(value, ranges) helper; runAlgorithm's argument count reduced from 5 to 4 by deriving order from code internally instead of passing it in; _t2/_t3's duplicated setup (hh/vi/ph/y) and final-result computation extracted into shared _t2t3Setup/_t2t3Result helpers, leaving each function's distinct series recurrence untouched. No behavior change — same algorithm, same values.
Code Health of src/dist/doubly-noncentral-beta.js improved from 9.38 to 10.0: _pdfRelocated/_cdfRelocated and their per-r inner sums _pdfTermSumOverS/_cdfTermSumOverS each had two forward/backward convergent-series walk loops flagged as a "Bumpy Road" (nested conditional logic repeated twice per function). Both loops of each pair now go through one shared _walkSum(z, count, term) accumulator, preserving the exact same accumulation order — the backward half's convergence check still compares against the combined forward+backward running total, not a fresh partial sum. No behavior change — same series, same values.
Code Health of src/dist/categorical.js improved from 9.38 to 10.0: skewness()/kurtosis() had a duplicated raw-moment accumulation loop, now extracted into a shared _rawMoments(order) helper that computes E[X], ..., E[X^order] in one pass over the pmf table. No behavior change — same formulas, same values.
Code Health of src/dist/irwin-hall.js improved from 9.38 to 10.0: _pdf/_cdf had a duplicated alternating-sign, log-domain term array construction and sort-then-Neumaier-sum, now extracted into a shared _alternatingLogSum(y, power) helper (power = n - 1 for the pdf series, power = n for the cdf series). No behavior change — same series, same values.
Code Health of src/dist/doubly-noncentral-t.js improved from 9.43 to 10.0: _findStartIndex's two nested loops (Fibonacci bracket search, then bisection) split into _bracketMaximum/_narrowBracket helpers, each taking the bracket as a single { j1, j2, f1, f2 } object rather than four separate arguments; _pdf's x*mu >= 0 forward/backward series computation extracted into _pdfSameSignSeries(x). No behavior change — same algorithm, same values.
Code Health of src/dist/davis.js improved from 9.38 to 10.0: mean()/variance()/skewness()/kurtosis() each recomputed the same raw-moment expression (b^k * Γ(n-k) * ζ(n-k) / (Γ(n)·ζ(n))), now extracted into a shared _rawMoments(maxOrder) helper. No behavior change — same formulas, same values.
Code Health of src/dist/noncentral-t.js improved from 9.58 to 10.0: fnm's AS243 series-constant setup (75-line method, over the 70-line "Large Method" threshold) extracted into a new _fnmSeriesInit(nu, delta, x) helper, leaving the forward/backward summations untouched. No behavior change — same series, same values.

Fixed

ran.special.besselKnu(nu, x): silently returned up to ~77% relative error for orders nu whose magnitude was comparable to x, just past the x=6 series/asymptotic crossover (e.g. nu=4.82, x=7.18 returned ≈0.000358 vs. the correct ≈0.00154), with no error, warning, or NaN to signal the defect. The unconditional dispatch to _KAsymptotic(nu, x) (the DLMF §10.40.2 large-x asymptotic expansion) for x > 6 ignored how nu compared to x; its "optimal truncation" only bounds error correctly when the expansion's first correction term (4ν²−1)/(8x) is already small, which fails once nu is comparable to x. Fixed by reducing the order to mu = |nu| - round(|nu|) ∈ [−0.5, 0.5] (where the existing connection formula and _KAsymptotic both stay accurate at any x) and reaching the target order via the same upward recurrence (DLMF §10.29.1) besselK already uses for integer order — no new algorithm was needed. Verified against mpmath (mp.dps=50) across nu ∈ [3,10], x ∈ [6,15] (#1361).
src/dist/gamma.js: Gamma.pdf(x) returned NaN for large alpha (e.g. alpha≈96.5), since the density was computed as two separate factors — Math.exp(logNorm - beta*x) * Math.pow(x, alpha-1) — and for large alpha/x the Math.pow term overflowed to Infinity while the Math.exp term underflowed towards 0, giving 0 * Infinity = NaN. _pdf(x) now accumulates the full exponent (logNorm - beta*x + (alpha-1)*Math.log(x)) before a single Math.exp() call, with an explicit x === 0 branch to avoid a new 0 * -Infinity NaN the log-space rewrite would otherwise introduce when alpha === 1 at the closed boundary (#1363).
src/dist/inverse-gamma.js: InverseGamma.pdf(x) silently returned 0 instead of a tiny representable density for extreme x (e.g. InverseGamma(0.01017360968553757, 0.22993683529824133).pdf(7.584718518060176e+162) returned 0 instead of mpmath's ≈2.927e-167), since _pdf(x) was computed as super._pdf(1 / x) / (x * x) and the x * x intermediate overflowed to Infinity for large x, making finite / Infinity collapse to 0. _pdf(x) now computes the density directly in log-space (logNorm - (alpha+1)*Math.log(x) - beta/x, exponentiated once), matching the class's own documented closed form and avoiding the overflow intermediate entirely (#1364).
.github/workflows/docs-deploy.yml: the release-channel cleanup step (find "$SITE" ... ! -name 'v*' -exec rm -rf {} +) deleted the gh-pages worktree's own .git metadata file on every tagged release, since .git starts with a dot and matches none of the exclusion patterns — every subsequent git command then failed with fatal: not a git repository, silently leaving the published docs site (and its / root, which is supposed to mirror the latest tagged release) stuck on the previous version. Added ! -name '.git' to the exclusion list, and added a workflow_dispatch version input as a recovery path to manually redeploy an already-tagged release without moving the (immutable) release tag.
scripts/generate-accuracy-docs.js: extractDistributionNames()/extractSpecialFunctionNames() had no sanity check on the export-line regex used to build docs/accuracy.md's coverage registry, so a future change to src/dist/index.js/src/special/index.js's export syntax could silently break the match and render the table missing rows for real distributions/functions with no error anywhere in CI. Both functions now throw if a known canary name (Gamma, digamma) is absent from the extracted list, failing loud instead of silently under-reporting coverage.
scripts/difftest-ci-gate.js: KNOWN_ISSUES suppressed a gate failure for the entry's name as a whole rather than for the specific failure reason that was originally allowlisted, so Gamma.pdf/InverseGamma.pdf (allowlisted for their tracked #1363/#1364 divergences) would have stayed silently green even if they developed a brand-new, unrelated failure reason (e.g. eval errors) on the same key. KNOWN_ISSUES now maps each entry to { issue, reasons }, and isGateFailure() fails the gate on any reason not in that list while the originally-allowlisted reason(s) stay suppressed; the job summary table annotates each reason individually as allowlisted/tracked or failing, rather than annotating the whole row (#1372).
scripts/difftest-ci-gate.js and scripts/generate-accuracy-docs.js: the KNOWN_ISSUES allowlist still tracked InverseGamma.pdf (#1364), Gamma.pdf (#1363), and besselK/besselKnu (#1140) as known accuracy defects even though all underlying issues were resolved by earlier PRs (#1364's own fix, #1375, #1216) that never triggered a recalibration. A fresh differential-testing sweep confirmed all four are now clean, so their entries were removed — both KNOWN_ISSUES maps are now empty; docs/accuracy.md was regenerated and the weekly CI accuracy gate (.github/workflows/difftest.yml) no longer silently allowlists a regression on any of them (#1377).
src/dist/_gamma.js: the shared small-shape (alpha < 1) gamma sampler's boost branch (gamma(r, a+1, b) * Math.pow(r.next(), 1/a)) underflowed to an exact 0.0 for a measurable fraction of draws at very small alpha (e.g. ~0.05% at alpha≈0.0102), a value outside Gamma's open (0, Infinity) support and, more visibly, one whose reciprocal in InverseGamma.sample() returned Infinity. The boost branch now rejects and redraws whenever the result underflows to exactly 0, protecting Gamma, InverseGamma, Beta, BetaPrime, StudentT, and every other distribution sharing this sampler. InverseGamma._generator() additionally guards against the reciprocal of a subnormal (nonzero but near-underflow) gamma draw overflowing to Infinity, resampling until the reciprocal is representable. BetaPrime._generator() and StudentT._generator() had the same subnormal-denominator overflow risk — BetaPrime's x / y ratio and StudentT's gamma(r, 0.5) / gamma(r, nu/2) ratio inside the sqrt both route through the boost branch when beta < 1 / nu < 2 respectively — and now resample the same way until the result is finite (#1379). That redraw-on-underflow guard itself became an unconditional infinite loop for alpha below ≈3.13e-13: at that scale, every one of the xoshiro128+ PRNG's 2^32 possible outputs underflows the boost factor, so the loop's acceptance probability is exactly zero and it never terminates (e.g. new Gamma(1e-15, 1).sample() hung forever). The boost factor is now computed in log-space (exp(ln(X) + ln(U)/a) instead of X * U^(1/a)) so it never underflows before combining with X, and an analytically-derived threshold short-circuits the provably-hopeless regime to the correctly-rounded 0 directly, with a generous iteration cap retained as a backstop above the threshold; InverseGamma and StudentT's own reciprocal/ratio rejection loops are capped the same way, returning Infinity (sign-adjusted for StudentT) — the IEEE-754-correct rounding of a value astronomically beyond Number.MAX_VALUE, per the return-value convention's "answer diverges" channel — instead of also looping forever. BetaPrime and Beta compose two independent gamma draws, so when both shape parameters are below the threshold the ratio is 0/0 (NaN) rather than a provable one-directional overflow; both now resolve the direction analytically (toward whichever boundary the smaller shape parameter's draw is pulled to, split randomly when the shapes are equal) before ever retrying, rather than falling through to a fixed sentinel. Two further consumers of Beta's sampler broke in a new way once it stopped hanging and started returning exact boundary values: BetaGeometric.sample() divided by Math.log(1 - 0) = +0 and returned -Infinity (a sign flip from the true +Infinity, outside its {1, 2, 3, ...} support), and BetaNegativeBinomial.sample() drove _poisson.js's lambda to Infinity, which fell through that function's large-lambda loop with no return and yielded undefined — an explicitly forbidden sentinel per this project's return-value conventions. Both are now guarded at their own call sites and return the correctly-rounded Infinity (#1384, see ADR-0054).
src/special/gamma.js: gamma(z) returned Infinity prematurely for z roughly in [143, 171] (e.g. gamma(143) returned Infinity where the true value is ≈2.695e245, well within double range), since the Lanczos tail computed Math.pow(t, z+0.5) and Math.exp(-t) as separate factors and the Math.pow term alone overflowed to Infinity before the Math.exp term could bring the product back down to its true finite value. The two factors are now combined into a single Math.exp((z + 0.5) * Math.log(t) - t), the same log-space technique log-gamma.js's Lanczos tail already uses, which never forms the oversized intermediate; gamma(171) (≈7.257e306) is now the largest finite integer argument, matching the true float64 overflow boundary at gamma(172).
src/special/beta.js: beta(x, y) could never return a negative value for negative non-integer arguments (e.g. beta(-0.5, -0.4), whose true value is ≈-1.249), since its fallback path composed three logGamma() calls — which intentionally return ln|Γ(z)|, discarding sign — with a single Math.exp(), and Math.exp() is always non-negative. Fixed by multiplying the result by the correct sign of each Γ factor, derived from Euler's reflection formula (Γ(z)Γ(1-z) = π/sin(πz)): sign(Γ(z)) = sign(sin(πz)) for non-integer z < 0, and +1 for z > 0.

[1.32.0] - 2026-08-09

Added

ran.dist.Distribution.prototype.copy(): returns a fully independent copy of a distribution instance, including its current PRNG state — a thin named wrapper around the existing this.constructor.load(this.save()) round-trip, added so cloning a Distribution instance doesn't require knowing that trick. Used internally by params()'s new Distribution-instance-valued-field cloning (see ### Fixed), and useful standalone — e.g. running two MCMC chains seeded from the same fitted distribution without them sharing PRNG state. See ADR-0051.
scripts/precision-refs-process.py and the generated test/precision-process.js: a stochastic-process precision gate, giving src/process/ the same arbitrary-precision verification standard src/dist/ already has from scripts/precision-refs-continuous.py/-discrete.py. Process densities were previously checked only against scipy doubles at a uniform 1e-10 over a handful of hand-picked points; the new gate covers all nine processes that expose a closed-form time-t marginal — AR1, BrownianBridge, BrownianMotion, CompoundPoisson, CoxIngersollRoss, GeometricBrownianMotion, OrnsteinUhlenbeck, Poisson, and RandomWalk — over a systematic 3-parameter-sets × 3-times × 5-interior-points grid, with the probe x-values obtained by inverting the high-precision marginal CDF at p ∈ {0.1, 0.3, 0.53, 0.72, 0.9} (integer lattice points for the discrete Poisson and RandomWalk). Each reference gates three independent code paths — pdf(x, t), marginal(t).pdf(x), and marginal(t).cdf(x), the last of which previously had no external reference at any tolerance; marginal() derives its law's parameters separately from pdf(), so checking the two only against each other (as test/process.js does at 1e-10) would let a shared parameterization slip cancel out. Every marginal law in the generator is re-derived from the process's own SDE or update rule rather than read off the JavaScript: CompoundPoisson's reference in particular is summed directly as a Poisson-weighted mixture of Gammas, never through the compound-Poisson → Tweedie parameter mapping that marginal() applies, so it gates that mapping as well as Tweedie's own Dunn & Smyth series. The generator self-checks 25 of those re-derivations against the values already vetted in test/process.js and verifies that all nine compound Poisson-gamma mixtures normalize to 1 and reproduce Wald's mean, aborting before emitting a single literal on any mismatch. Seven of the nine processes hold at 1e-14 with no exception; RandomWalk at p = 0.3 (3e-14 pdf / 2e-14 cdf, log-gamma ULP amplification at t = 30) and CompoundPoisson (6e-14 pdf, Tweedie series — its cdf stays gated at 1e-14, since the two floors genuinely diverge) carry documented, mechanism-named bounds pinned just above their measured worst case. No process behavior changed: this is a pure regression guard, not a bug fix (#1223).
scripts/check-subpath-runtime.js (npm run check-subpath-runtime): dynamically imports one representative built ESM subpath module from each of the three minified categories (dist/beta.esm.js and dist/poisson.esm.js for distributions, dist/process/brownian-motion.esm.js for processes, dist/mc/rwm.esm.js for MCMC samplers) and asserts instantiation succeeds, constructor.name survives minification, and a known computed value matches — a pdf/cdf value against an mpmath/scipy-sourced reference for Beta/Poisson/BrownianMotion, and a seeded, pinned sample() array (in addition to its shape) for RWM. This is a direct regression guard for the keep_classnames: true fix in #1220, wired into CI's build job (.github/workflows/ci.yml), since npm test only ever exercises src/ and never imports from dist/ (#1227).
ran.process.Process.fit(path, dt) (static, per subclass): estimates process parameters from an observed discrete-time path, added as a throw-by-default hook on Process (mirroring marginal(t)'s rollout, #1132) and implemented for BrownianMotion, GeometricBrownianMotion, and OrnsteinUhlenbeck via their exact closed-form MLE — increments (or log-returns, or the AR(1) transition already coded into OrnsteinUhlenbeck._next()) are i.i.d./exactly linear-Gaussian, so sample mean/variance (or OLS regression of X_{n+1} on X_n) recovers the parameters to machine precision as the path grows. CoxIngersollRoss.fit() instead uses two-stage Conditional Least Squares (Overbeck & Rydén, 1997), since CIR's true one-step conditional transition is a scaled noncentral chi-squared with generally non-integer degrees of freedom — a different object from the Gamma marginal already implemented as its pdf(x,t)/marginal(t) (valid only because the class hardcodes x0 = 0) — and ran.dist.NoncentralChi2 rounds its k to the nearest integer, so it cannot represent CIR's non-integer degrees of freedom for a true conditional MLE. CLS is consistent but not maximally efficient, and its accuracy degrades near the Feller boundary and at large dt. See ADR-0044 (#1133). Extended to AR1.fit(path) (OLS regression of X_{n+1} on X_n, reusing the shared ols() helper — the true transition has no intercept, but fitting through the intercept-plus-slope form still recovers phi consistently since the true intercept is exactly 0); RandomWalk.fit(path) (the exact MLE: the fraction of +1 steps among all observed increments, algebraically identical to recovering p from the sample mean of increments since every step is exactly ±1); and BrownianBridge.fit(path, T, dt) (the exact MLE for sigma, since each step's conditional variance is fully determined by the known, fixed T/dt — unlike the other four processes, T is a required argument here rather than something to estimate, since the bridge's defining feature is a fixed, given endpoint). AR1 and RandomWalk have no dt parameter in their own model, so their fit() drops it entirely rather than taking an unused argument (#1212). Extended to the counting-process family: Poisson.fit(path, dt) recovers the exact MLE lambda = totalCount / (n*dt) from the path's net increase, since increments are i.i.d. Poisson(lambda*dt). CompoundPoisson.fit(path, dt, jumpDistConstructor) estimates lambda the same way, treating every non-zero increment as exactly one jump — individual arrival counts within a single dt interval are not observable from the cumulative path alone, so this is an approximation valid when lambda*dt is small enough that multi-jump intervals are rare — and fits the jump-size distribution's own parameters by handing the recovered non-zero increments to the caller-supplied jumpDistConstructor's static fit() (#1213).
ran.process.Process.prototype.lnL(path): transition log-likelihood of an observed discrete-time path, added as a throw-by-default hook on Process (mirroring marginal(t)'s and fit(path, dt)'s partial rollout) and implemented for BrownianMotion, OrnsteinUhlenbeck, and GeometricBrownianMotion via a new protected _transitionLnPdf(xPrev, xNext) hook that each overrides with its closed-form one-step Gaussian (BM, OU) or log-Gaussian-with-Jacobian (GBM) transition density — the same law already encoded in each class's _next() and reused by fit()'s sufficient statistics, so no new numerical machinery was needed. Computes transition, not marginal, likelihood: a realized path's points are a dependent, Markov-correlated sequence, not independent draws from the marginal distribution, and conflating the two is a known trap in this codebase (CoxIngersollRoss's marginal/conditional Gamma mismatch, #1133) — see ADR-0046. GeometricBrownianMotion.lnL() returns -Infinity (not a thrown error) for a path that visits a non-positive state, mirroring pdf(x,t)'s existing x <= 0 => 0 convention (#1153).
ran.dist.Tweedie(mu, phi, p): the Tweedie exponential dispersion model for the compound Poisson-Gamma power range 1 < p < 2 — a point mass at zero (P(Y=0) = exp(-lambda)) plus a continuous positive tail, used for insurance claims, rainfall accumulation, and zero-inflated continuous GLM responses. Neither the PDF nor the CDF has a closed form: _pdf evaluates the Dunn & Smyth (2005) infinite series for the compound Poisson-Gamma density in log-space (all terms are positive for 1 < p < 2, so no cancellation), locating the peak term via a closed-form Stirling estimate before summing; _cdf sums a Poisson-weighted series of gammaLowerIncomplete evaluations with a purely relative convergence check (no absolute floor, avoiding the false-early-convergence failure mode documented for DoublyNoncentralBeta, #1108). Both series are capped a number of terms past their peak that scales with sqrt(peak) rather than by a constant, since the peak's own width grows the same way — a constant slack silently truncates both sums mid-peak once the peak clears MAX_SERIES_ITER (at Tweedie(50, 0.02, 1.5), lambda = 707, it left pdf 0.5% low, cdf plateauing at 0.970 instead of reaching 1, and q(p) returning NaN above that plateau). _generator() samples via the exact compound Poisson-Gamma simulation (N ~ Poisson(lambda), then the N events' total drawn as a single Gamma(N * shape, rate), which is an identity rather than an approximation and keeps a sample at O(1) instead of O(lambda)); _q(p) returns 0 for any p <= P(Y=0) (the base class's root-finder cannot find a sign change in that region, since cdf(x) - p >= 0 everywhere) and root-finds otherwise; mean()/variance()/skewness()/kurtosis() are closed-form via EDM cumulant theory; _fitInit() seeds p at the literature-typical 1.5 (no closed-form estimator exists) with method-of-moments for mu/phi (#1136).
ran.dist.ExponentiallyModifiedGaussian(mu, sigma, lambda): the exponentially modified Gaussian (EMG) distribution, the convolution of a Normal(mu, sigma^2) and an Exponential(lambda) random variable — used for right-skewed data with exponential tails (chromatography peak modeling, reaction-time analysis, neuroscience). PDF/CDF use the closed-form erfc-based formula (Wikipedia: Exponentially modified Gaussian distribution), rewritten via the scaled complementary error function erfcx to avoid the exp(large)·erfc(large→0) cancellation the naive formula hits for large lambda·sigma — the same technique already used for InverseGaussian's CDF. _generator() samples as the sum of independent Normal and Exponential draws; mean()/variance()/skewness()/kurtosis() are closed-form; _fitInit() uses method-of-moments (#1131).
ran.process.Process.prototype.marginal(t): returns the process's marginal distribution at time t as a fully-functional ran.dist.Distribution instance, unlocking quantile(), hazard(), survival(), likelihood(), aic(), bic(), and test() on process marginals without any new numerical machinery. Implemented by composing each process's already-existing mean()/variance()/pdf() formulas: BrownianMotion, OrnsteinUhlenbeck, and BrownianBridge return Normal; GeometricBrownianMotion returns LogNormal; CoxIngersollRoss returns Gamma, reusing the shape/scale already derived for its own pdf() — valid since the process always starts at x0 = 0, which collapses the general noncentral-chi-squared transition density to a plain Gamma. Throws for t outside the domain where the marginal is genuinely a continuous distribution (t <= 0 for all five; additionally t >= T for BrownianBridge, where the process is pinned to a point mass) (#1132). Extended to Poisson and AR1, which return ran.dist.Poisson/Normal instances the same way and likewise throw for t <= 0 (the target class's own parameter validation can't express the degenerate zero-mean/zero-variance case at t = 0); and to RandomWalk, which returns an instance of a new private ShiftedBinomial distribution (src/dist/_shifted-binomial.js, not part of the public ran.dist API — see ADR-0045) representing the pushforward of Binomial(t, p) under x = 2k - t. Unlike Poisson/AR1, RandomWalk.marginal(0) does not throw, since a point mass at 0 is directly representable as ShiftedBinomial(0, p) (#1156). CompoundPoisson (and its deprecated alias CompoundPoissonProcess) now overrides marginal(t): for a ran.dist.Gamma jumpDist, X_t is by definition the compound Poisson-gamma total that ran.dist.Tweedie already represents, so marginal(t) returns a Tweedie instance via a closed-form parameter mapping derived from matching each representation's Poisson rate and gamma shape/rate — no new special function or Distribution subclass was needed, since Tweedie already shipped in #1136. Every other jumpDist throws a specific, documented error instead of inheriting the generic base-class message: an arbitrary caller-supplied distribution makes X_t a Poisson mixture over sums of an unknown distribution, with no general closed form reducible to a single existing ran.dist class (#1157).
"engines": { "node": ">=20" } added to package.json, documenting the Node.js version constraint that CI's test matrix and nyc@18 (engines.node: "20 || >=22", see #960) already impose in practice, so npm/Yarn warn contributors and downstream consumers installing on Node 18 or earlier instead of failing later with a confusing nyc internal error (#1137).
.github/dependabot.yml: weekly automated npm devDependency updates, restricted to patch-level bumps (minor/major are ignored, since devDependency major bumps like ESLint 7→9 or a Rollup major often carry breaking config/plugin-API changes that warrant manual review), grouped into babel, lint, test, build, and docs buckets to keep PR volume low while staying atomic and reviewable; each PR runs through the existing CI gates (lint, jsdoclint, test+coverage, typecheck, docs-build, build) before merge (#1142).
Versioned API docs: the published site now mirrors the latest tagged release at /, keeps every past release permanently at /vX.Y.Z/, and publishes tip-of-main at /unreleased/ with an "unreleased" banner — instead of redeploying the entire site from whatever was on main on every push (which had let unreleased distributions such as Tweedie leak into the live docs ahead of their release). A version dropdown and an "outdated release" banner are populated client-side from a versions.json manifest. See decisions/0043-versioned-docs-deployment.md.
ran.test.cramerVonMises(values, cdf, alpha): the Cramér-von Mises single-sample goodness-of-fit test, testing the null hypothesis that values is drawn from the distribution cdf represents. The statistic T = n·ω² = 1/(12n) + Σᵢ[(2i-1)/(2n) − F(xᵢ)]² is computed over sorted, CDF-transformed order statistics (the same EDF-comparison family as the private andersonDarling helper in src/dist/_tests.js, but with squared-deviation rather than log-weighted terms); the asymptotic p-value sums the Csörgő & Faraway (1996, JRSS-B 58(1), eq. 1.2) convergent series for the n → ∞ limiting distribution's CDF, built entirely from besselKnu/logGamma already in src/special/ — no new special function or algorithm was needed. Returns {stat, passed, pValue}, the shape adopted by ADR-0042 for single-sample GoF tests newly exported from ran.test (extending, rather than replacing, the plain {stat, passed} shape the module's existing multi-sample comparison tests use) (#1134).
ran.test.kolmogorovSmirnov(x, y, alpha): the two-sample Kolmogorov-Smirnov test, testing the null hypothesis that x and y are drawn from the same distribution. The statistic D = sup|F1(x) - F2(x)| is computed over the pooled empirical CDFs of the two samples by binary-searching each sample's sorted values at every pooled point; the asymptotic p-value is obtained from the existing ran.dist.Kolmogorov distribution's survival(), evaluated at sqrt(n1*n2/(n1+n2))·D. Returns {stat, passed, pValue} per ADR-0042 (#1138).
ran.test.andersonDarling(values, cdf, alpha): the Anderson-Darling single-sample goodness-of-fit test, testing the null hypothesis that values is drawn from the distribution cdf represents. The statistic A² = -n - (1/n)·Σᵢ(2i-1)[ln F(xᵢ) + ln(1-F(x_{n+1-i}))] is computed over sorted, CDF-transformed order statistics; the asymptotic p-value uses the Marsaglia & Marsaglia (2004, JSS 9(2)) rational-function approximation to the limiting distribution, with their finite-sample correction. This is a thin public wrapper around the private andersonDarling helper already implemented and tested in src/dist/_tests.js (which continues to back Distribution.prototype.test() unchanged, with its own hardcoded α=0.01) — no new math was needed. Returns {stat, pValue, passed} per ADR-0042, which explicitly named this function as the next to adopt that shape (#1144).
ran.process.Process.prototype.params(): returns the process's parameters (this.p), mirroring ran.dist.Distribution.prototype.params() so that .fit() results and other downstream consumers can inspect a process's parameters through a stable public accessor instead of reaching into the internal this.p storage convention (#1251).
ran.dist.guess(data, options): fits a set of candidate distributions to a dataset and ranks them by BIC weight — Δᵢ = BICᵢ − BIC_min, wᵢ = exp(−0.5·Δᵢ) / Σⱼ exp(−0.5·Δⱼ) — the estimated probability that each candidate is the best-fitting model in the set, given the data. "Guess" is intentional: this is a heuristic exploratory tool, not a verdict. Candidates are pre-filtered before the expensive fit() call: hard filters exclude type (continuous/discrete) and support mismatches, and soft, statistically-principled filters exclude symmetric-only or positive-skew-only families against sample skewness, Exponential-like families against an out-of-range coefficient of variation, and Poisson-like/NegativeBinomial-like families against an incompatible dispersion index. Throws if data.length is below 20 * max_k (BIC's asymptotic approximation needs roughly 20 observations per parameter, evaluated against the largest parameter count among surviving candidates), and skips (rather than propagates) any candidate whose fit() throws. Returns a sorted array of {name, params, bicWeight, pValue}, carrying a warning string property when every surviving candidate fails goodness-of-fit at α=0.05. The default candidate pool covers all distributions, including VonMises, Rice, NoncentralChi2, NoncentralChi, and Skellam — an initial exclusion for their per-point Bessel-function evaluation cost was lifted after benchmarking showed their fit() cost is comparable to already-included distributions of the same parameter count (#813, #1051). The soft filters' false-exclusion rates are empirically measured by Monte Carlo simulation (scripts/guess-filter-validation.js): the skewness filter's ~5% analytical target was initially confirmed for Normal (4.2%-4.7% measured) but found badly miscalibrated for Laplace (34.7%-51.1% measured, 7-10× the target) under a single normal-only threshold (2·√(6/n)) shared across every SYMMETRIC family; the threshold is now computed per family as 2·√(c/n), where c is each family's own asymptotic skewness-estimator variance (Normal → 6, Uniform → 72/35, Laplace → 63, derived from Var(g1)·n ≈ μ6/μ2³ − 6·μ4/μ2² + 9), bringing measured false exclusion to 4.2%-4.7% for Normal, 4.4%-5.6% for Uniform, and 1.4%-4.2% for Laplace (#1054, #1064); the coefficient-of-variation and dispersion-index filters measured ~0% false exclusion for their representative distributions, well within safe bounds.
ran.shape.max/ran.shape.min are now exported from src/shape/index.js. Both files existed with public-style JSDoc (@memberof ran.shape) but were only reachable via direct relative imports (e.g. from src/dispersion/range.js), not through the public ran.shape namespace — missing wiring, not a missing implementation (#1233).
ran.dist.WrappedCauchy(mu, rho): the wrapped Cauchy circular distribution, the standard heavy-tailed alternative to VonMises, parameterized by mean direction mu and concentration rho in (0, 1). Unlike VonMises, whose CDF requires an infinite Bessel-function series, wrapped Cauchy's PDF, CDF, and quantile are all elementary closed forms built from sin/cos/tan/atan2 — no new special functions were needed. Support is the mu-centred window [mu-pi, mu+pi] (matching scipy's vonmises(loc=mu) convention) rather than a fixed [-pi, pi], since a circular distribution has no canonical cut point independent of its own location parameter; _cdf uses atan2 (rather than a plain atan ratio) to avoid the tan((x-mu)/2) singularity at the support boundary. mean()/variance()/skewness()/kurtosis() are left to the base class's numerical quadrature fallback (matching VonMises's precedent) since these are the arithmetic, not circular, moments and are always finite on the bounded support. _fitInit() uses the trigonometric moment estimator (mean resultant length/angle), since no closed-form MLE exists in general (Kent & Tyler, 1988) (#1135).

Changed

ran.dist.VonMises gains a location parameter mu (the mean direction), matching the parameterization on Wikipedia: pdf(x) = exp(kappa*cos(x-mu)) / (2*pi*I0(kappa)), with support [mu-pi, mu+pi] instead of the previously fixed [-pi, pi]. The constructor signature changes to new VonMises(mu, kappa), matching every other location-shape distribution in the library (e.g. Cauchy(x0, gamma)); .k (the parameter count .aic()/.bic() penalize against) is now 2, was 1. _fitInit() recovers both parameters from the sample's circular resultant vector: mu as its angle, re-anchored to the 2*pi "sheet" nearest the sample's own extremes and clamped into [xmax-pi, xmin+pi] so the fixed-width support is guaranteed to contain every sample (mirroring how Uniform/Triangular derive their own support-defining parameters directly from the data extremes — needed so ran.dist.guess()'s pre-fit probe never excludes VonMises over an estimation-noise-driven support miss); kappa is unchanged, still from the resultant length.
ran.dist._tests.chi2(values, pmf, c) and ran.dist._tests.andersonDarling(values, cdf) (and therefore Distribution.test() for both discrete and continuous distributions) now return a pValue field alongside the existing statistics/passed fields. chi2PValue() and andersonDarlingPValue() — sibling helpers that briefly exposed this without changing the parent functions' return shape — are removed now that both parents carry the field directly; ran.dist.guess()'s per-candidate p-value now reads chi2(...).pValue/andersonDarling(...).pValue instead (#1052, #1053).
Hot-path _pdf/_cdf/_generator/_q methods on 14 distributions now read parameter-only constants (log-gamma normalizers, log-binomial/log-beta terms, Bessel/Poisson-mixing terms) from this.c instead of recomputing them on every call: Gamma (and its subclasses Chi2, Erlang, which now share the parent's cached log-normalizer instead of each calling logGamma again), InverseChi2, Poisson, NegativeBinomial, NoncentralChi2, NoncentralBeta (also speeding up NoncentralF, which delegates to it), DoublyNoncentralBeta, BetaBinomial, NegativeHypergeometric, Hypergeometric, Muth, and VonMises (which also caches the ratio-of-uniforms sampling constant used by _generator()). BrownianMotion, OrnsteinUhlenbeck, and GeometricBrownianMotion's _transitionLnPdf hot path (called once per step from Process.prototype.lnL(path), potentially many times in an MLE-calibration/MCMC loop) likewise now reads its precomputed log-scale constant (this.c.logSigmaDt/this.c.logNoise) instead of calling Math.log() on every transition. No behavior or return-value change.
mean()/variance()/skewness()/kurtosis() on 24 distributions now share cached parameter-only raw/central moment terms (gamma/beta/Hurwitz-zeta/Riemann-zeta evaluations, series sums) instead of each method recomputing them independently: Frechet, GeneralizedExtremeValue, InvertedWeibull, Weibull (and DoubleWeibull, which now reuses Weibull's cached terms instead of calling gamma() again), Burr, Kumaraswamy, GeneralizedLogistic, FisherZ, Zeta, BenktanderII, HeadsMinusTails, Hyperexponential, ShiftedLogLogistic, Soliton (caches the harmonic number instead of re-summing an O(N) loop per method), UniformProduct, JohnsonSU, LogLogistic, LogSeries, LogGamma, LogLaplace, and ExponentiatedWeibull. GeneralizedNormal and HalfGeneralizedNormal now read GeneralizedGamma's already-cached log-gamma terms instead of bypassing the cache with their own logGamma() calls. No behavior or return-value change.
The per-distribution, per-process, and per-MCMC-sampler subpath builds (dist/<name>.esm.js, dist/process/<name>.esm.js, dist/mc/<name>.esm.js) are now minified with @rollup/plugin-terser (module: true, preserving ESM-safe mangling for downstream tree-shaking), the same way dist/ranjs.min.js already was — these were previously emitted with full variable names, JSDoc, and whitespace intact. keep_classnames: true is set (at a negligible size cost) since Distribution.load()/Distribution.fit() (src/dist/_distribution.js) and HMC/NUTS's resumed-state validation (src/mc/_mcmc.js) interpolate this.name/this.constructor.name into thrown error messages — without it, minification would silently replace e.g. Beta.fit() requires a _fitInit()... with a mangled single-letter class name in every subpath-imported distribution's error output. A representative subpath build (dist/beta.esm.js) shrinks from 93106 to ~22973 bytes raw (-75%) and from 30113 to ~10144 bytes gzipped (-66%), matching the single-distribution-import path README.md recommends as the low-footprint usage pattern (#1220).
ran.process.AR1.marginal(t) no longer performs its own variance(t) <= 0 pre-check, matching the pattern every other process's marginal() already used (BrownianMotion, BrownianBridge, OrnsteinUhlenbeck, CoxIngersollRoss, GeometricBrownianMotion, Poisson, PoissonProcess, CompoundPoisson, RandomWalk all construct their target law straight from mean(t)/variance(t) and let its constructor validate the scale). The guard's only real-world trigger was the variance() cancellation bug fixed earlier in this same release, which returned exactly 0 for near-unit-root phi with small fractional t — so it was converting a silent precision defect in its own dependency into a confusing AR1.marginal(): variance is not positive at t domain error rather than protecting against a genuinely non-positive variance. A 29700-combination sweep of variance(t) (dense phi grid straddling the 1e-14 reformulation boundary, sigma and t spanning underflow through overflow) found no strictly negative result for any t > 0; the explosive |phi| >= 1 branch diverges to +Infinity but never flips sign, since its numerator and denominator change sign together. v <= 0 remains reachable only by floating-point underflow (t below ~1e-322, or sigma below ~1.6e-161 so sigma² underflows), and those inputs are still rejected with an Error — now Invalid parameters. ... sigma > 0 from Normal's own validation, so only the message changes (#1244). pdf(x, t)'s parallel v <= 0 => NaN guard is deliberately left in place: it predates the guard under discussion and uses a different return channel.

Deprecated

ran.dist.VonMises's single-argument constructor form new VonMises(kappa) (implicitly mu = 0, the library's previous fixed behavior) is deprecated in favor of new VonMises(mu, kappa). The old form still constructs and behaves identically but emits a one-time console.warn on first use; it will be removed in v1.33.0.

Fixed

ran.dist.Skellam(mu1, mu2).cdf(k) lost 3-4 orders of magnitude of precision (5e-10 to 6e-9 relative error, vs. the ~1e-12 to 1e-14 floor elsewhere) for k close to mu1 in highly asymmetric configurations (e.g. Skellam(5000, 1).cdf(k) for k in [4988, 4997]). Contrary to the issue's initial suspicion, src/special/marcum-q.js's _transitionBand is not implicated — for this call shape (marcumQ(k+1, mu2, mu1) with mu2 < 30), the dispatcher always routes through _series, whose only non-recurrence value is a single gammaUpperIncomplete(mu, mu1) call. The bug is entirely in src/special/gamma-incomplete.js's _gui (the upper-incomplete-gamma continued fraction): (1) its loop was capped at the shared MAX_ITER=100 with no regime-aware extension, unlike its sibling _gli, silently truncating before the ~150-160 iterations the near-diagonal s≈mu1≈x regime needs (the same failure class #1286 fixed in _fc); (2) its shared prefactor with _gli, f * Math.exp(-x + s*Math.log(x) - logGamma(s)), cancels three O(mu1)-magnitude terms down to an O(1) result, an unavoidable ~1e-11 to 1e-12 floor no compensated summation of those specific terms can beat. Both are fixed via a new src/special/_deviance.js module (log1pmx, relocated verbatim from marcum-q.js's private _log1pmx; stirlerr, the Stirling series remainder; bd0, the Loader (2000) binomial-deviance term) that lets _gli/_gui compute f * Math.sqrt(s/(2*Math.PI)) * Math.exp(-bd0(s,x) - stirlerr(s)) with every intermediate term O(1) or O(log s) instead of O(s), plus a regime-aware iteration budget and a throw-on-non-convergence guard (_assertGuiConverged, mirroring _fc's _assertFcConverged, ADR-0049) for _gui. bd0 routes only x/s near 1 through the cancellation-safe log1pmx path; far from 1 it uses the direct x - s - s*Math.log(x/s) (no cancellation there, and routing extreme ratios through log1pmx(x/s - 1) would itself lose accuracy, since x/s - 1 rounds to exactly -1 once x is ~16 orders of magnitude below s). Deriving _gui's iteration budget also surfaced a second, unrelated latent bug: for s near zero (not just large s), the continued fraction needs up to ~99 iterations at the x=s+1 boundary regardless of how small s is — previously silently wrong (caught live by Tweedie.test()'s Anderson-Darling sweep once the new throw guard was in place); _gui's floor is raised from MAX_ITER=100 to 200, empirically confirmed ≥2x the worst-case measured need across s from 1e-20 to 20000. Skellam(5000,1).cdf(k) for k in [4988,4997] now matches mpmath (mp.dps=50) to ~1e-14 to 2e-15 relative error (previously up to 6e-9); the [5000,1]/[1000,1]/[2000,1] precision-gate groups' tolerances are unchanged since their floor is now set by Skellam._pdf's own, separate log-space cancellation (#1321, ~9e-12 worst case), not by this fix. scripts/precision-refs-discrete.py's [5000,1] k-grid gains two points inside the previously-withheld band (k=4990, 4995) (#1348).
ran.dist.Skellam(mu1, mu2).pdf(x) returned NaN for highly asymmetric mu1/mu2 (e.g. Skellam(1000, 1).pdf(999)) with x near the mean, distinct from and un-fixed by #1309's earlier symmetric-large-mu overflow fix. _pdf multiplied three independently-scaled factors -- expNegScaled (exp(-(√mu1-√mu2)²), which underflows to exactly 0 once the asymmetry between mu1 and mu2 grows large, contrary to a doc comment inherited from #1309's fix, which only holds for the symmetric case), Math.pow(sqrtRatio, x) (overflows to Infinity), and besselIExpScaled(|x|, twoSqrtProd) (also underflows to exactly 0, since the true scaled Bessel value at this Bessel order/argument combination -- e.g. order 999 against argument ~63.25 -- is genuinely non-representable as a double, ~1e-1092) -- a three-way 0 * Infinity * 0 collision even though the true pmf is a normal, representable number (~0.01-0.2). ran.special.bessel.js gains logBesselIExpScaled(n, x), the log-domain analogue of besselIExpScaled: it delegates to besselIExpScaled and takes its log whenever that stays representable, falling back to a convergence-checked Taylor-series evaluation in log-space (leading term via the already-exported logGamma) only when besselIExpScaled underflows to exactly 0 -- purely additive, with zero change to besselIExpScaled's own behavior or precision-gated callers. Skellam._pdf now combines all three log-space terms into a single exponent and calls Math.exp exactly once, matching the codebase's established convention for this shape of computation (Poisson, Borel, Delaporte, etc.), rather than multiplying three separately-materialized factors. Combining terms whose individual magnitude grows with mu1 while their sum stays O(1) near the mean does cost some precision at very large mu1 (measured worst case ~5.7e-13 relative error in pdf at mu1=1000, up to ~6e-12 at mu1=5000) -- an inherent, honestly-documented trade-off (_LOG_CANCEL tolerance override in scripts/precision-refs-discrete.py), and a dramatic improvement over the prior NaN. Closes #1321.
ran.special.besselISpherical(n, x) threw a confusing "_hi: continued fraction failed to converge for n=..., x=... after NaN iterations" for n > 1 and negative x with |x| >= 1 (the branch that delegates to the Wronskian-based continued-fraction helper _hi). _hi's iteration budget computes Math.ceil(7 * Math.sqrt(x)), which is NaN for negative x, so its for loop's condition was always false and the loop never ran even once — the throw fired on an un-iterated ratio, not on genuine non-convergence. besselISpherical(n, x) is entire with only x^(n+2k) terms in its Taylor series, so it has definite parity i_n(-x) = (-1)^n i_n(x); the default branch now maps negative x to (n % 2 === 0 ? 1 : -1) * besselISpherical(n, -x) before reaching _hi, returning the mathematically correct value instead of throwing. Not reachable through any production call path — NoncentralChi/NoncentralChi2, the only internal callers, always pass a non-negative argument — only reachable via a direct call to the exported besselISpherical(n, x) (#1324).
ran.dist.NoncentralT's internal CDF helper (fnm, an AS243-series implementation) rounded to exactly 1.0/0 whenever the true survival probability was closer to the boundary than a double can represent — not a fixable precision bug in fnm itself (no double "1 minus something" can distinguish a gap smaller than ~1.11e-16), but a caller-visible information loss whenever two such saturated values were differenced (or summed and then subtracted from 1). This broke ran.dist.DoublyNoncentralT.pdf(x) in the x*mu < 0 branch at extreme parameters: DoublyNoncentralT(5, 5, 120).pdf(-0.7) remained ~1.7x off from its true value even after #1235's cancellation fix (the fix's own documented residual limitation). NoncentralT gains a direct survival sibling, snm(nu, mu, x) (computed via tanh-sinh quadrature over the noncentral-t's mixture representation, never as 1 - fnm(...)), which DoublyNoncentralT._pdfPoissonMixture now falls back to for any Poisson-mixture term whose fnm difference cannot be trusted (gated on nu magnitude, where fnm's own regularizedBetaIncomplete-derived series genuinely loses precision, and on the raw difference's magnitude) — DoublyNoncentralT(5, 5, 120).pdf(-0.7) now matches the mpmath reference (8.08e-15) to ~1e-14 relative error. The same root cause independently broke .cdf(x) for x < 0 at the same extreme parameters — _cdf sums Poisson-weighted fnm terms directly and subtracts from 1, so high-weight terms saturating to exactly 1.0 silently overcounted (DoublyNoncentralT(5, 5, 120).cdf(-0.7) returned 6.66e-16 against an mpmath reference of 2.62e-16, ~154% relative error) — found while validating the .pdf() fix above; _cdf now accumulates the x < 0 complement termwise (sum(weight_i * (1 - fnm_i)), falling back to snm under the same gating) and matches the mpmath reference to ~1e-10 relative error, with no regression to .fit()/quantile-root-finding performance (#1250). The identical saturation was confirmed directly on ran.dist.NoncentralT.pdf(x) itself (not just DoublyNoncentralT's use of it): NoncentralT(30, 5).pdf(40) returned exactly 0 while the mpmath reference is ~1.54e-18, since _pdf's own nu * (fnm(nu+2, mu, x*nuScale) - fnm(nu, mu, x)) / x differences two fnm calls that both saturate to exactly 1. _pdf now routes through the same nu-magnitude/diff-magnitude-gated snm fallback (reusing DoublyNoncentralT's thresholds verbatim), matching the mpmath reference to ~3e-15 relative error with no change to any ordinary (non-saturating) NoncentralT evaluation (#1302). Separately, that same nu-magnitude/diff-magnitude gate (as originally shipped by #1250, before the fix described next) had two further blind spots in DoublyNoncentralT._fnmDiff/_cdfTerm, both closed under #1298: (1) _fnmDiff missed a single "knife-edge" nu0 per x, where one of the two fnm calls being differenced had separated from fnm's phi = 0.5*(1+erf(-mu/sqrt2)) plateau and the other hadn't — their raw difference was then dominated by the still-plateaued operand's own error, which is wrong but not small (~1e-7, evading a < 1e-9 magnitude check) — solely responsible for DoublyNoncentralT(5, 5, 120).pdf(-0.2)'s remaining ~2e-3 relative error; (2) _cdfTerm missed an entire low-nu0 plateaued range whose raw complement is pinned at exactly 1 - phi (~2.87e-7 for mu=5, also not < 1e-9) — solely responsible for DoublyNoncentralT(5, 5, 120).cdf(-0.1) being ~14.5x wrong, a case #1298 itself did not anticipate (its own acceptance criteria assumed cdf was unaffected, having only measured cdf(-0.2)). Both helpers now check two independent conditions, since a raw fnm value can be untrustworthy either way and neither implies the other: whether it is still stuck at phi (no nu-magnitude pre-filter needed — this only fires when the nu-dependent correction is genuinely unresolved), or — the original #1250 mechanism, still needed since a value that has resolved away from phi can independently saturate toward the opposite 0/1 boundary as nu grows — the pre-existing nu0 >= 30 && |raw value| < 1e-9 magnitude check. pdf(-0.2) and cdf(-0.1) — the two points issue #1298 itself reported broken — now match their mpmath references to ~1.9e-14 and ~4.9e-14 relative error respectively (worst case across all three reported points: pdf 8.75e-14, cdf 3.80e-8); test/guess.js's fit-all-distributions tests rose from a post-#1250 baseline of ~23-24s to ~46-50s in isolation, since the added phi-check fires more often during .fit()'s optimizer exploration than the magnitude check alone did — combined with #1302's own new, independent NoncentralT._fnmDiff cost (above), this pushed both tests past their previous 60000ms mocha timeout under full-suite --parallel CPU contention (isolated runs stayed under 60s; the full suite did not), so both timeouts were raised to 120000ms (test/guess.js, matching .mocharc.yml's own global default). NoncentralT._fnmDiff (added by #1302, above) reused the original magnitude-only gate and was NOT updated with this phi-check, so NoncentralT.pdf(x) still silently returned 0 (or, in a nearby regime, a badly wrong nonzero value) whenever both fnm calls stayed stuck at phi without ever separating — confirmed at NoncentralT(5, 6).pdf(-0.5) (returned 0, mpmath reference ~3.34e-10) and NoncentralT(1, 8).pdf(-0.3) (also 0, reference ~4.78e-16), both at nu far below the 30 floor the magnitude gate needs to even evaluate; NoncentralT(10, 6).pdf(-1.0) returned a nonzero but ~480x wrong value, showing the blind spot isn't only an exactly-zero case. Porting DoublyNoncentralT's corrected phi-equality gate verbatim was not sufficient on its own: NoncentralT.snm (its designated fallback) is only accurate for nu >= 30, per its own documented limitation, and NoncentralT._pdf's call site — unlike DoublyNoncentralT's, which never invokes snm below that floor — needs it down to nu = 1. NoncentralT._fnmDiff is removed; _pdf now inlines the corrected gate (phi = 0.5*(1+erf(-mu/sqrt2)), computed unconditionally — the sign-flip fnm's own internal x<0?-mu:mu uses is fully internal to that function's x>=0 ? z : 1-z return-value flip and does not propagate to callers) and, when it fires, falls back to a new NoncentralT._pdfDirect(nu, mu, x): a direct tanh-sinh quadrature of the density's own defining formula (already documented in the class JSDoc) rather than a CDF difference, so there is no cancellation to lose precision to at any nu. All three reported cases now match their mpmath references to ~1e-14-1e-15 relative error, with test/guess.js's .fit()-exploration timing unaffected (#1318). Separately, _pdf's other saturation gate, nearOppositeBoundary (nu >= 30 && |a - b| < 1e-9, unchanged by #1318), missed a large-nu regime the flat 1e-9 threshold was never tuned for: fnm's own absolute noise floor grows roughly linearly with nu, and the fast path's nu * (a - b) / x identity amplifies that noise by the same nu/x factor, so NoncentralT(10000, 0).pdf(0.5) returned 0.3520526413036684 against a true 0.35205267468981716 (~9.5e-8 relative error, nine orders of magnitude worse than _pdfDirect's own ~1e-13) while |a - b| = 1.76e-5 sailed straight past the flat threshold. nearOppositeBoundary's threshold is now scaled by nu (nu * Number.EPSILON * 1e10, empirically validated across nu from 30 to 100000), correctly routing large-nu evaluations to _pdfDirect while leaving the already-accurate nu in [30, 300] regime #1318 validated untouched. This made NoncentralT.fit() pay _pdfDirect's ~80x per-call cost whenever Powell's optimizer explores large nu — harmless for genuinely noncentral-t-shaped data (small interior optimum, few such evaluations), but data with no good t fit (e.g. bounded/circular samples) has no interior optimum in nu and drove the unbounded search into the tens of thousands, multiplying that cost across hundreds of thousands of likelihood evaluations (test/guess.js's VonMises-in-default-pool test regressed from ~34s to over 150s). NoncentralT gains a static _powellOptions() bounded search budget ({ tol: 1e-3, maxIter: 15 }), mirroring the identical DoublyNoncentralBeta/DoublyNoncentralF fix for the same class of problem (#1063) — cuts the pathological case back to ~9s while reproducing genuinely-matched-data fits' converged (nu, mu) to within floating-point noise (#1325). DoublyNoncentralT._pdfPoissonMixture's own _fnmDiff helper had the structurally identical flat-1e-9 nearOppositeBoundary gate, never updated by #1325 (whose scope was restricted to NoncentralT._pdf) — each Poisson-mixture term multiplies its fnm-difference by nu0 (the term's own degrees of freedom), the same amplification shape as NoncentralT._pdf's nu*(a-b)/x, so DoublyNoncentralT.pdf(x) accumulated the identical nu-scaled precision loss at large nu. Porting #1325's nu * Number.EPSILON * 1e10 threshold into _fnmDiff's gate (keeping its existing nu0 >= 30 guard) closes a real, already-reachable gap: DoublyNoncentralT(5, 2, 120).pdf(-0.7) tightens from ~1.7e-9 to ~7.3e-15 relative error, a ~235,000x improvement at parameters this library's own precision-gate suite already exercises. At genuinely extreme nu (>= 10000, unreachable via .fit() or any realistic dataset) the fix only partially helps — _fnmDiff's fallback (NoncentralT.snm(lo) - NoncentralT.snm(hi)) is itself a difference of two ~1e-11-accurate quadratures, unlike NoncentralT._pdf's cancellation-free _pdfDirect fallback, so it re-encounters a smaller-scale version of the same cancellation problem one level down once the true difference itself shrinks to a comparable magnitude — DoublyNoncentralT(50000, 0.01, 0.1).pdf(-0.5) improves from ~3.15e-6 to ~6.0e-7 relative error but is not made fully precise; this residual is a documented, accepted limitation, not a regression. _cdfTerm (used by .cdf()) does not share this amplification (its Poisson-mixture sum does not multiply by nu0) and is left unchanged. The wider-firing gate has a second effect, caught only by the full test suite (not by any test targeting the fix itself): .fit()'s Powell optimizer, on data with no interior optimum in nu/theta (e.g. the same VonMises(0,2)-sampled data #1325 used), now pays the added NoncentralT.snm-fallback cost across hundreds of thousands of likelihood evaluations — DoublyNoncentralT.fit() on that data went from ~6s to ~68s with an unbounded search, exactly the class of regression #1325's own solution doc warned a hot-path-to-expensive-fallback fix must be checked for separately. DoublyNoncentralT gains a static _powellOptions() ({ tol: 1e-2, maxIter: 15 }, matching DoublyNoncentralBeta's values), bounding the pathological case back to ~18s alone / ~34s inside guess()'s full default-pool sweep, with no intolerable quality loss on well-matched data. See solutions/correctness/2026-08-04-0823-doubly-noncentral-t-nu-scaled-fnmdiff-gate-fix.md (#1332).
src/algorithms/powell.js's fractional convergence test (2*|fStart-fret| <= tol*(|fStart|+|fret|)) tolerates an absolute log-likelihood gap that grows with sample size n, since Distribution.fit()'s objective is -lnL(data): issue #1338 measured this across every _powellOptions()-bounded distribution and found DoublyNoncentralT(5,1,2)'s bounded-vs-unbounded gap growing roughly with n, from ~0.12 at n=100 to ~3.08 at n=3000, and DoublyNoncentralF(3,8,1,1)'s ranging non-monotonically from ~0.74 at n=100 to ~2.48 at n=3000 (peaking at ~3.51 at n=1000) — both non-trivial and not shrinking with more data. powell() gains an optional capAbs field (default Infinity, so every existing caller not passing it is unaffected) that bounds the threshold via Math.min(tol*(|fStart|+|fret|), capAbs), and Distribution.fit() now merges in a calibrated capAbs=2 default — chosen via Wilks'/LRT theory (the lnL gap at a confidence-region edge is ~chi2_p/2, an O(1) quantity independent of n) and confirmed against every affected distribution's own worst-case pathological-data wall-clock/call-count ceiling — unless a subclass's own _powellOptions() already sets capAbs itself. Closes DoublyNoncentralT's gap from ~1.41/~3.08 to ~0.0003/~0.018 at n=1000/3000, and DoublyNoncentralF's from ~3.51/~2.48 to ~0.002/~0.045 at the same sample sizes; a no-op for NoncentralT (its 2-parameter (nu, mu) gap is already ~1e-11 to 1e-13 at every n) and only a partial improvement for DoublyNoncentralBeta, consistent with part of its gap being a genuine shape/noncentrality ridge (#1063) rather than purely a convergence-tolerance artifact. DoublyNoncentralF.fit()'s own custom ridge-penalized objective calls powell() directly rather than through Distribution.fit(), so it does not receive the injected default. See solutions/testing/2026-08-05-1736-powell-fractional-convergence-n-scaling.md (#1342).
ran.dist.Distribution.prototype.params() and ran.process.Process.prototype.params() returned this.p by live reference, letting a caller silently corrupt a distribution's or process's internal state (e.g. const p = dist.params(); p.mu = 999). Both now return a shallow copy ({ ...this.p }); nothing in the codebase relied on the previous mutable-reference behavior. See ADR-0047 (#1257). The same live-reference issue was found in ran.dist.Distribution.prototype.support(), which fed the mutable boundary objects directly into pdf/cdf/quantile/sample's internal _belowSupport/_aboveSupport/_atClosedBoundary checks; it now returns this.s.map(b => ({ ...b })), copying the nested {closed, value} boundary objects as well as the array, since a shallow array spread alone would still leave them shared. The shallow { ...this.p } copy itself left one gap: array-valued parameter fields (Hyperexponential's weights/rates, Categorical's weights) were still shared by reference, so dist.params().weights[0] = 0 still reached this.p.weights through the copied top-level key. Both params() implementations now additionally copy every array-valued field (Array.isArray(p[key]) ? [...p[key]] : p[key]), a targeted per-field copy rather than a generic recursive/structured clone — the latter would also try to clone non-array object fields such as CompoundPoisson's jumpDist (a live Distribution instance), which ADR-0047 scoped out of this accessor's copy guarantee on the reasoning that "a caller mutating a nested distribution's own state goes through that distribution's own params()/setters, not through the outer process's accessor." See ADR-0050 (#1299). That reasoning turned out to be wrong: CompoundPoisson._next() samples directly from this.p.jumpDist on every step with no per-step reseed (only CompoundPoisson.prototype.seed() reseeds it, once, at seed time), so jumpDist's PRNG stream is live process state, not an isolated implementation detail — confirmed empirically, seeding two identically-constructed processes the same way but calling .seed() on one's params().jumpDist in between produced different path() output from the other, with neither process's own .seed() called again. params() now also clones any Distribution-instance-valued field (via the new copy() method, above), superseding ADR-0047's carve-out for this field shape specifically; its shallow-copy decision for plain primitive/array fields is unaffected. See ADR-0051.
ran.dist.ReciprocalInverseGaussian.cdf(x) returned a value quantized to multiples of 2^-53 (essentially garbage) for small x, where the internal 1 - InverseGaussian.cdf(1/x) subtraction catastrophically cancelled because InverseGaussian.cdf(1/x) rounds to within 1 ULP of 1 in that regime. InverseGaussian gains a numerically stable _survival(x) (mirroring its own _cdf's erfc/erfcx cancellation fix, applied symmetrically), which ReciprocalInverseGaussian.cdf(x) now calls instead of subtracting from 1.
test/dist-cases-continuous.js's Normal[0,2] far-tail (x = ±14) refVals were stale — 1 ULP off for pdf, ~2.3e-6 relative error for cdf — predating the cancellation-safe far-tail fix already shipped for test/precision-continuous.js under #808, which was never back-ported to this file. scripts/precision-refs-continuous.py's self_check() (only made to actually run under #1110) caught the discrepancy; the correct values were independently re-derived and confirmed via three agreeing mp.dps=50 formulations (erf, erfc, mpmath's built-in ncdf) (#1193).
ran.special.marcumQ/ran.special.marcumP returned NaN in the quadrature branch (large x, deep lower tail) whenever the scaled argument y/mu was far below 1 — _zetaxy()'s saddle-point formula catastrophically cancelled once sqrt(1 + 4*x*y/mu²) rounded to exactly 1.0, collapsing a denominator to 0. This broke ran.dist.Rice.cdf(x)/.q(p), ran.dist.NoncentralChi.cdf(x)/.q(p), and ran.dist.NoncentralChi2.cdf(x)/.q(p) near x = 0 and, for .q(p), at any probability p — the base class's quantile root-finder always probes cdf(Number.EPSILON) first, and the resulting NaN silently defeated the root-finder's own bracket-validity guard (NaN comparisons are always false in JS). _zetaxy now uses the exact identity d1 - eps = d2 to fold the two near-cancelling terms into one well-conditioned expression whenever 4*x*y/mu² < 0.5, leaving the existing near-transition-line formula (y/mu close to x/mu + 1) unchanged (#1179).
scripts/precision-refs-continuous.py --emit --allow-prune --only Name1,Name2 (dev-only tooling) silently ignored --only and recomputed every distribution instead of scoping to the named ones, because --only's parsing checked a fixed argv position while --allow-prune's was already position-independent. --only is now detected by argv.index('--only') in both the --emit and self-check branches, so it works regardless of where it appears relative to --allow-prune.
scripts/precision-refs-continuous.py's existing_groups() (dev-only tooling), the guard render() relies on to preserve hand-maintained precision-gate groups it can never reproduce (e.g. TruncatedExponential), silently dropped any group whose raw text didn't match its expected name: '...', params: ..., tol: ... shape instead of preserving or flagging it — a future hand-edited group with different field order or an inserted field would then be neither reproduced nor preserved, reintroducing the exact silent-loss failure mode this mechanism exists to prevent. It now raises RuntimeError naming the unparseable span so a maintainer can fix it before --emit runs.
scripts/precision-refs-continuous.py's bare/--check self-check (dev-only tooling) hung for 100+ minutes once it reached DoublyNoncentralBeta's LARGE_LAMBDA_ANCHORS regression case ((2,2,1200,1200)), never completing and never reaching the remaining ~90 distributions — dncbeta_cdf() called mpmath's expensive regularized-incomplete-beta function (betainc) fresh for every one of the ~800k-1M (r, si) pairs its nested double-Poisson-mixture summation visits at this lambda scale (dncbeta_pdf(), which needed no such call, was never the bottleneck). dncbeta_cdf() now tracks the incomplete-beta value itself via an exact recurrence (a standard DLMF 8.17.20-style contiguous relation, independently re-derived and numerically verified against direct betainc() calls at both toy and production scale before use) instead of recomputing it from scratch at every step, cutting DoublyNoncentralBeta(2,2,1200,1200).cdf(0.3) from ~1235s to ~66s and .cdf(0.5) from ~2659s to ~67s with no change to the walk's structure, floor, or convergence semantics (the #1108/#1086 anti-regression fix), and no change to any already-vetted reference value. self_check() --only DoublyNoncentralBeta now completes in ~4 minutes with 0 mismatches (#1194).
npm run standard/npm run lint silently skipped every file sitting directly in src/ or test/ (e.g. src/index.js, test/ad.js, test/core.js, test/algorithms.js) because the lint/standard scripts passed an unquoted src/**/*.js test/**/*.js glob to the shell — under a POSIX /bin/sh/dash shell (how npm actually invokes scripts on Linux, absent bash's non-default globstar option), ** behaves like a single *, so only files exactly two path segments deep were ever linted. Both scripts now quote the globs ('src/**/*.js' 'test/**/*.js') so standard's own bundled glob engine expands ** correctly instead of the shell. Fixing the scope surfaced several previously-hidden, genuinely-live violations, now fixed: an over-precision numeric literal in test/ad.js (no-loss-of-precision) shortened to the value that round-trips exactly as a double, two similarly over-precision refVals/moments reference literals in test/dist-cases-continuous.js corrected the same way, and two new SomeClass(...) calls used only for their deprecation-warning side effect in test/process.js (no-new) now capture the instance into a variable and assert instanceof on it.
ran.process.CoxIngersollRoss.pdf(0, t) returned +Infinity when the Feller condition is violated (alpha < 1), disagreeing with the Gamma(alpha, 1/scale) instance marginal(t) returns for the same process, whose own pdf(0) is 0 there — Gamma's support (like Beta's and Weibull's) is open at 0 whenever the shape parameter is below 1, so the boundary point is excluded rather than evaluated. pdf(0, t) now returns 0 for alpha < 1, matching marginal(t).pdf(0); the already-correct alpha === 1 (1/scale) and alpha > 1 (0) cases are unaffected.
ran.dist.NoncentralBeta.pdf(1) returned 0 for beta < 1 instead of the correct +Infinity. The density carries a (1 - x)^(beta - 1) factor that diverges as x → 1 when beta < 1 (dominated by the k = 0 Poisson term regardless of alpha/lambda), but the Poisson-mixture series evaluated at exactly x = 1 produced Infinity - Infinity = NaN, which the base pdf() silently collapsed to 0 via its NaN→closed-boundary guard. _pdf now short-circuits x === 1, beta < 1 to Infinity; beta >= 1 is unaffected ((1 - x)^(beta - 1) is 0 for beta > 1, or 1 for beta === 1, giving the finite Poisson mean alpha + lambda/2, both handled correctly by the existing series). The mpmath reference generator (scripts/precision-refs-continuous.py, dev-only) had the mirror-image bug — a blanket x >= 1 → 0 early return that never inspected beta — and now returns +inf/alpha + lambda/2/0 for beta < 1/beta == 1/beta > 1 respectively (#1121).
ran.core.Xoshiro128p.next() is uniform on [0, 1) and can legitimately return exactly 0 (~1-in-2³² per call). Six generators fed that raw draw straight into Math.log(...), which sends Math.log(0) = -Infinity through the rest of the formula and can leak a literal Infinity (or, for UniformProduct, a silent 0 that violates its open lower bound) as a returned sample: the shared _exponential() helper (and therefore Exponential and HyperExponential), YuleSimon, UniformProduct, LogSeries, FlorySchulz, and PolyaAeppli. All six now take 1 - r.next() instead of r.next() into the log, which is uniform on (0, 1] and can never hit the singularity at 0. LogSeries's default GoF test seed sweep changes from [0, 42, 12345] to [1, 42, 12345] since the fix necessarily changes the deterministic sample sequence for a given seed, and seed 0's new sequence happened to land in the chi-squared test's ~1% rejection region by chance (empirically confirmed: sample mean matches theory, and a 200-seed sweep shows a ~2% failure rate consistent with the test's own false-positive rate, not a systematic bias).
ran.dist.BetaRectangular's parameter count (.k) was inherited from Beta's constructor (2) despite BetaRectangular having 5 free parameters (alpha, beta, theta, a, b), causing .aic()/.bic() to under-penalize its complexity. .k now correctly reports 5. A follow-up audit of every reparametrizing Distribution subclass found the same defect in 11 more distributions and fixed all of them: PERT (3, was 2 from Beta), Bates (3, was 1 from IrwinHall), BetaBinomial (3, was 2 from Categorical), SkewNormal (3, was 2 from Normal), BirnbaumSaunders (3, was 2 from Normal), JohnsonSB (4, was 2 from Normal), and JohnsonSU (4, was 2 from Normal) all under-counted their true free-parameter count; Gilbrat (0, was 2 from LogNormal/Normal), PowerLaw (1, was 2 from Kumaraswamy), QExponential (2, was 3 from GeneralizedPareto), and R (1, was 2 from Beta) went the other way — each fixes one or more of its parent's parameters to a constant, so the inherited .k over-counted and over-penalized complexity (#1049). A further audit of every remaining Distribution subclass extending a concrete distribution class found the same under-counting defect in 2 more Categorical subclasses: Hypergeometric (3, was 2 from Categorical) and NegativeHypergeometric (3, was 2 from Categorical); every other such subclass was confirmed to already report the correct .k (#1094).
ran.dist.PowerLaw, R, Gilbrat, JohnsonSU, JohnsonSB, SkewNormal, BirnbaumSaunders, and PERT — reparametrizing Distribution subclasses that call super(...) with transformed or dummy values — leaked the parent constructor's internal parameter keys (and, for PowerLaw/R/Gilbrat, values the caller never supplied) into the public .params() method instead of exposing only the constructor's own declared natural parameters; BirnbaumSaunders additionally stored its location parameter under the wrong key mu2 instead of its declared mu, so .params().mu always returned the leaked 0 rather than the constructor's actual value. .params() now returns exactly the natural parameters named in each constructor's JSDoc, matching the fix already applied to Chi2/Erlang/MaxwellBoltzmann/Rayleigh/DoubleWeibull/HalfNormal/Slash/LogCauchy/StudentZ under ADR-0018 (#1057). ran.dist.QExponential — the one distribution deliberately left out of that fix, since it previously relied on GeneralizedPareto's canonical {mu, sigma, xi} for its moment methods and IEEE-754 divergence-boundary tests — now follows the same convention: .params() returns {q, lambda}, and the relocated GP-space values live in this.c, with no change to pdf/cdf/quantile results. Bringing skewness()/kurtosis() in line with GeneralizedPareto's own three-tier formula/Infinity/NaN split surfaced a latent discrepancy between the two: for xi >= 1/2 (variance itself infinite, e.g. q = 1.8), QExponential returned Infinity where GeneralizedPareto, given the identical xi, correctly returns NaN for the same indeterminate ∞/∞ ratio (decisions/0015-return-value-and-error-conventions.md); QExponential.skewness()/.kurtosis() now return NaN in that range, matching GeneralizedPareto (#1058). The same leak is fixed for the remaining 9 reparametrizing subclasses: ran.dist.F, BaldingNichols, Weibull, NoncentralF, DoublyNoncentralF, GeneralizedGamma, GeneralizedNormal, DoublyNoncentralChi2, and ExponentiatedWeibull. Weibull and GeneralizedNormal had the same wrong-key-collision pattern as BirnbaumSaunders: Weibull.params().lambda returned the leaked dummy 1 passed to the internal Exponential(1) transform while the constructor's real scale was hidden under a synthetic lambda2; GeneralizedNormal.params().alpha/.beta were similarly shadowed by leaked Gamma-space values, hidden under alpha2/beta2 (ExponentiatedWeibull, which reparametrizes Weibull, inherited the same lambda/lambda2 split and is fixed alongside it). DoublyNoncentralChi2.params() no longer exposes the internal collapsed k/lambda it computes internally (DoublyNoncentralChi2(k1,k2,λ1,λ2) ≡ NoncentralChi2(k1+k2,λ1+λ2)) alongside its own k1/k2/lambda1/lambda2. NoncentralF, DoublyNoncentralF, and DoublyNoncentralChi2 — whose immediate parent's pdf/cdf/sampling are non-trivial series algorithms rather than a one-line special-function call — now cache a correctly-parameterized instance of that parent and delegate to it (ADR-0039), instead of duplicating its internals or modifying the parent class (NoncentralBeta, DoublyNoncentralBeta, NoncentralChi2 are themselves independent public distributions, unaffected). ran.dist.HalfGeneralizedNormal, which extends GeneralizedNormal, is updated alongside it since it read the same leaked keys directly (#1070). HalfGeneralizedNormal itself was inadvertently left out of both that effort's and #1057/ADR-0018's scoped file lists: its own constructor never reassigned this.p after super(0, alpha, beta), so .params() returned the inherited { mu: 0, alpha, beta } instead of its own two natural parameters. .params() now returns exactly { alpha, beta }; since GeneralizedNormal.prototype._generator/_pdf/_cdf read this.p.mu directly, HalfGeneralizedNormal's own overrides of those three methods are now inlined against the mu = 0-folded formulas (mirroring the Weibull/Exponential pattern) instead of delegating to super, with no change to sampled values, pdf/cdf results, or .seed()-determined output (#1087).
The generated API docs (npm run docs) now render the individual fields of every options-object constructor (e.g. ran.mc.RWM's options.logDensity, options.config, options.initialState; ran.mc.HMC's additional options.gradLogDensity) as indented rows in the Parameters table, instead of silently dropping them behind a single opaque options: Object row. documentation.js nests dotted @param tags (e.g. @param {Object} options.config) under the parent param's properties array rather than returning them as flat top-level params; docs/src/param-parser.js never read that array, so every JSDoc'd nested field for every options-object constructor in the codebase (RWM, AdaptiveMetropolis, HMC, NUTS, MALA, Gibbs, Slice, ParallelTempering) was invisible in the rendered docs even though it was correctly documented in the source JSDoc. docs/index.js's call-signature renderer is updated alongside to use only the top-level (depth-0) params, so signatures still read e.g. RWM(options) rather than incorrectly listing the newly-surfaced nested fields as separate positional arguments.
ran.dist.DoublyNoncentralBeta.fit() (and DoublyNoncentralF.fit(), which delegates its _pdf/_cdf to DoublyNoncentralBeta) could take 13-30+ seconds on ordinary data, driven by two compounding issues in doubly-noncentral-beta.js: (1) _pdfRBackward/_cdfRBackward's Poisson-mixing outer loop had no iteration cap, unlike its MAX_ITER-bounded forward counterpart, so it could run arbitrarily long as Powell's optimizer explored large trial non-centrality parameters — now capped at MAX_ITER to match; (2) more significantly, on data that does not genuinely belong to this family, the log-likelihood surface carries a long, near-flat ridge between the shape and non-centrality parameters that a full-precision Powell search (the base class's default tol=1e-8, maxIter=200) chases almost indefinitely — worse, each step further along the ridge is itself more expensive to evaluate, since larger non-centrality parameters require more series terms. DoublyNoncentralBeta now overrides static fit() with a bounded Powell search budget (tol=1e-2, maxIter=15), empirically verified to still recover parameters within this class's existing fit tolerances on well-matched data (matching the default optimizer's result within ordinary finite-sample noise across multiple seeds) while bounding worst-case cost to roughly 1-2s (#1063).
ran.dist.DoublyNoncentralBeta.pdf()/.cdf() (and DoublyNoncentralF, which delegates to it) returned NaN instead of a finite probability once both non-centrality parameters were large (e.g. lambda1 = lambda2 = 2000), driven by two compounding overflow/underflow bugs in the double-Poisson-mixture summation: (1) the Poisson-weight speed-up constants pr0/ps0 were computed as the unnormalized lambda^k/k! with the compensating e^{-lambda} deferred to a later multiplication, overflowing Number.MAX_VALUE once lambda1/lambda2 exceeded ~1418 — before the compensator was ever applied; (2) independently, Beta(alpha+r0, beta+s0) underflows to exact 0 in double precision once both r0 = round(lambda1/2) and s0 = round(lambda2/2) are large (e.g. Beta(1002,1002) ≈ 1e-604, far below Number.MIN_VALUE), while power-of-x/power-of-y terms in the same term are comparably extreme in the opposite direction — even though the combined term is an ordinary, representable double, 0 * Infinity (or equivalent) produced NaN once these isolated linear-space factors combined. The Poisson-weight normalization now defers to a single outer-scale multiplication (bit-identical to the prior, more precise formulation) whenever it is safe to do so, only folding it in directly when the unnormalized magnitude would itself overflow; the Beta-function constant and power-of-x/power-of-y terms are now tracked as logarithms, updated additively through the existing forward/backward recurrence, and combined via a single exp() per term rather than ever being materialized in isolation. pdf/cdf are now finite for lambda1 = lambda2 up to at least 50000, matching the issue's acceptance criteria, with existing small-lambda precision-gate values unchanged (#1075). Series-truncation precision at very large lambda was left out of scope for that fix and separately tracked as #1063/#1086: pdf/cdf could return a finite-looking but silently wrong value — off by up to ~10 orders of magnitude — once lambda1 + lambda2 ≳ 400-600 and x moved away from 0.5 (e.g. DoublyNoncentralBeta(2,2,1200,1200).pdf(0.3) previously returned 9.5e-31 against an mpmath (dps=50) reference of 3.03e-21). Two compounding truncation bugs are now fixed: (1) the outer Poisson-mixing loops (_pdfRForward/_pdfRBackward/_cdfRForward/_cdfRBackward) were capped at MAX_ITER (100) steps from the x-independent Poisson mean (r0, s0), but the true summand peak shifts away from (r0, s0) as x moves from 0.5 (e.g. a shift of ~146 steps for lambda1=lambda2=1200, x=0.3) — now capped at the wider MAX_SERIES_ITER (500), matching the cap already used elsewhere for this class of series; (2) more fundamentally, the inner per-r sum over s (_pdfSumOverS/_cdfSumOverS) relied on the shared recursiveSum helper's convergence check, which floors its relative-error tolerance at EPS * max(|sum|, 1) — an absolute floor that falsely declares convergence after only 1-2 terms whenever a sum's true converged value is itself far below 1 in magnitude (routine here, since these densities can be astronomically small). A new file-local _seriesSum helper drops that floor (safe here specifically because every summed term is a non-negative probability-weighted value, never subject to cancellation), fixing the truncation at its root rather than merely widening the outer loop's window. DoublyNoncentralF.fit()'s bounded Powell search budget (#1063) is unaffected — empirically re-verified at ~114000 _pdf calls and ~8s on the original #1063 reproduction, matching the pre-fix baseline. A residual gap remained even after the MAX_SERIES_ITER widening: once lambda1 + lambda2 grows large enough (empirically >= ~8000) combined with x far enough from 0.5, the true peak shifts beyond even that wider window, and pdf()/cdf() silently returned exactly 0 — not merely imprecise, flatly and incorrectly zero for parameter combinations already within this class's own tested range (#1102). _pdf/_cdf now detect this case directly (the standard walk's own last term stays non-negligible relative to its total after exhausting its window, rather than estimating in advance whether relocation will be needed — an estimated-shift heuristic was tried first and found to misroute cases the standard window already handles correctly) and fall back to a walk centered on a closed-form peak-index estimate instead of (r0, s0), bounded by a separate, smaller iteration cap (RELOCATE_MAX_ITER) chosen specifically to keep this fallback's inherently costlier per-term evaluation from reintroducing the #1063 fit()-search-cost regression. This fallback trades some precision for that bound — large-lambda values a few x away from 0.5 are now correct to within an order of magnitude rather than exact-0, not to full machine precision; already-correct small/moderate-lambda behavior is unchanged.
ran.dist.Distribution.load(state) restored this.p/this.c directly from a serialized state with no shape validation, so loading a malformed or version-skewed snapshot (e.g. one saved before a distribution migrated its this.p/this.c split under ADR-0018) silently read missing keys as undefined and propagated to NaN from pdf()/cdf()/sample() instead of throwing. load() now constructs a throwaway probe instance from the restored params (padded to the constructor's declared arity, so distributions like Categorical whose this.p intentionally holds fewer keys than constructor arguments are still validated correctly) and compares its this.p/this.c key sets against the restored state's, throwing a clear Error on any mismatch before the state is otherwise used unchanged. Because the probe runs the real constructor, load() can also throw on a snapshot whose this.p/this.c shape is unchanged but whose saved values now violate a constructor constraint that has since been tightened (e.g. a parameter that used to allow >= 0 now requires > 0) — an intentional, accepted trade-off, not a regression (#1074, decisions/0038-distribution-load-probe-validation.md's "Consequences → Harder" section).
ran.dist.NoncentralChi(1, lambda).pdf(0) hardcoded a return of 0 for every k, but the true limit at k=1 is finite and nonzero (sqrt(2/pi)*exp(-lambda^2/2), since only the underlying non-central chi-squared pdf's j=0 Poisson term diverges as v^(-1/2) near v=0 for df=1) — matching the fix already applied to ran.dist.Chi(1).pdf(0). k >= 2 is unaffected, since the true limit there is genuinely 0 (#1122).
ran.dist.DoublyNoncentralF's constructor built its internal DoublyNoncentralBeta delegate (the one pdf()/cdf()/sample() actually compute against) from raw, un-rounded d1/d2, while .params() reported the rounded integers its own JSDoc promises — a silent internal/public mismatch that also broke save()+load() round-trips for non-integer inputs, since _afterLoad() rebuilt the delegate from the rounded post-restore params instead of the original raw ones. d1/d2 are now rounded once, before any internal use, matching the pattern already used by NoncentralF/DoublyNoncentralChi2, so .params(), pdf()/cdf()/sample(), and a save()+load() round trip are now always internally consistent. Rounding early on its own discretizes the log-likelihood surface fit()'s Powell search explores, re-triggering the #1063 bounded-search regression at roughly double the _pdf call count; DoublyNoncentralF now overrides static fit() to search DoublyNoncentralBeta's continuous space directly (only rounding the final returned instance) with a smooth squared-hinge penalty — zero within a plausible region of the moment-matched initial guess, growing quadratically only beyond it — that keeps Powell off #1063's near-flat ridge without ever excluding a finite parameter value outright, so genuinely large-parameter fits are never silently underfit the way a hard cutoff would (#1084).
ran.dist.VonMises.cdf(x) (and therefore .q(), whose root-finder samples cdf() at arbitrary internal points) could return values far outside [0, 1] for concentrated distributions (kappa gtrsim 6-9) whenever x was at or near a multiple of pi/4 — e.g. VonMises(9).cdf(-Math.PI / 4) returned -0.0074 instead of 0.0119, and VonMises(9).q(VonMises(9).cdf(-1)) returned -pi/4 instead of -1. The underlying Fourier-series summation checked convergence on each raw term, which happens to collapse to machine-epsilon at x = k*pi/4 (sin(4x) ≈ 0 there) well before the series had actually converged for concentrated kappa; convergence is now checked on the term's non-oscillating envelope instead, which cannot be fooled by an incidental zero of sin(i*x).
ran.special.besselI(0, x) (and therefore ran.dist.Rice, VonMises, Skellam (at k=0), and NoncentralChi/NoncentralChi2 (at k=2) wherever the effective Bessel argument fell in the same range) was off by up to ~1.2e-9 relative error for x in roughly (10, 14], well outside the library's usual ~1e-14 precision — a "cold start" gap immediately after _besselIBackward's Miller backward-recurrence takes over from the |x| <= 10 Taylor series, recovering smoothly by x ~ 15-16. The recurrence's run-up-margin formula scales its extra headroom as sqrt(40 * n), which degenerates to exactly 0 for n = 0 (the order besselI(0, x) dispatches to) while every n >= 1 order already receives adequate margin from the same term; n is now clamped to Math.max(n, 1) inside that formula, so n = 0 inherits n = 1's already-validated margin with zero behavioral change for any n >= 1. Also corrects a pre-existing self-referential reference literal in test/special.js's |x|=10 routing-boundary test (it asserted a value computed from the pre-fix buggy code path instead of mpmath), and adds the Rice[3.16,1]/NoncentralChi[2,3.5]/NoncentralChi2[2,8]/Skellam[6,5] precision-gate parameter sets that issue #1143's boundary-grid work deliberately withheld because they surfaced this exact gap (#1185).
ran.dist.DoublyNoncentralT.pdf(x)'s general (mu != 0) branch returned significantly wrong densities (up to ~13% relative error observed) once mu was non-zero and large relative to nu, combined with large theta — e.g. DoublyNoncentralT(5, 5, 120).pdf(1.3) returned 0.8149681936132279 against an mpmath (mp.dps=50) reference of 0.71818185584468099.... The series walk advanced Kummer's ₁F₁(a,b,z) across the series index via a three-term contiguous recurrence in a (_f11Forward/_f11Backward), which is numerically unstable in both directions once the series' peak index pushes a large relative to b — confirmed by direct measurement (forward: 330% error from the very first recurrence step; backward: growing to 9 orders of magnitude of error near the series' start). Both private methods are removed; every series term now calls the already-correct f11() special function directly, matching the mpmath reference to ~1e-11 to ~1e-15 relative precision. See solutions/correctness/2026-07-30-1600-doubly-noncentral-t-pdf-f11-recurrence-instability.md (#1207).
ran.dist.DoublyNoncentralT.cdf(x) returned badly wrong, non-monotonic values once theta was large enough that exp(-theta/2) underflowed below Number.EPSILON (e.g. DoublyNoncentralT(5, 5, 120).cdf(-1) returned 1 while .cdf(0) returned ~1.5e-31) — the Poisson-mixture summation's leading term satisfied recursiveSum's default absolute-floor convergence check after a single iteration, the same failure mode previously fixed for DoublyNoncentralBeta (#1086/#1103). Fixed by passing { useFloor: false }, the opt-out recursiveSum gained for that earlier fix. Discovered, and the boundary-adjacent DoublyNoncentralT[5, 0, 120] precision-gate parameter set added, while extending #1143's boundary-grid methodology to f11's |z|=50 dispatch threshold (issue #1189).
ran.special.besselInu(nu, x) returned Infinity for very negative fractional order (e.g. nu = -1.5, -2.5, -3.3) at x near the ~710 series-overflow boundary, even though the true value is a large but finite number (e.g. besselInu(-1.5, 709) returned Infinity against an mpmath (dps=50) reference of ~1.23e+306) — the internal recursiveSum accumulator representing the series sum before the (x/2)^nu prefactor is applied overflowed past Number.MAX_VALUE, since for very negative nu that prefactor is tiny and the unnormalized sum must be proportionally larger to compensate. besselInu now uses a hand-written loop that rescales the running sum and current term in lockstep (mirroring _besselIBackward's existing overflow-guard pattern) whenever the sum approaches double overflow, tracking a log-scale offset combined into the final result only when a rescale actually occurred — preserving the original direct-multiplication precision for every case that never needs it, including besselKnu's connection-formula cancellation path (#1215).
ran.test.hsic() and ran.test.mannWhitney() silently mis-calibrated their Type-I error rate, discovered via new Monte Carlo calibration tests added while extending the hypothesis-test suite's rigor bar (#1229). hsic() fit a Gamma null approximation following Gretton et al.'s hsicTestGamma.m reference, whose b parameter is computed in MATLAB's shape/scale convention (Gamma mean = a*b), but passed it directly as ran.dist.Gamma's rate parameter (mean = a/rate) without inverting it, and additionally queried the lower alpha-quantile instead of the upper (1-alpha)-quantile appropriate for HSIC's right-tailed test (large statistic implies dependence) — together these suppressed the empirical Type-I error to ~0% instead of the nominal 5% (200-trial simulation: 0/200 rejections under H0 before the fix). Now uses new Gamma(a, 1 / b).q(1 - alpha); re-simulation gives 12/200 (6%, consistent with alpha=0.05). mannWhitney() compared its already-folded U = min(U1, U2) statistic against Normal(0,1).q(1 - 2*alpha), but a folded two-sided statistic's correct critical value is the alpha/2-tail (P(U1<=c or U2<=c) = 2*Phi((c-m)/s) = alpha implies z = q(1-alpha/2)) — the original formula inflated empirical Type-I error to ~17.5% (35/200 rejections under H0 before the fix). Now uses Normal(0,1).q(1 - alpha / 2); re-simulation gives 11/200 (5.5%). Both fixes are verified against the pre-existing seeded regression tests (hsic's dependent-data rejection, mannWhitney's same/different-distribution pass/reject cases), which are unaffected.
ran.dist.DoublyNoncentralT.pdf(x) had large relative error (up to ~130x observed) whenever x*mu < 0, even after #1207 replaced the unstable ₁F₁ recurrence in the same branch with direct f11() calls. The branch summed a series that alternates sign when x*mu < 0, accelerated via wynnEpsilon; series acceleration cannot recover precision already lost to cancellation between individual terms many orders of magnitude larger than the converged sum. _pdf's x*mu < 0 branch now uses a new private _pdfPoissonMixture(x), the term-by-term derivative of the cancellation-free Poisson(theta/2)-mixture-of-noncentral-t formula _cdf already uses — every term is a Poisson weight times a difference of two NoncentralT.fnm CDF values, never an alternating-sign term. The x*mu >= 0 branch is unchanged. See solutions/correctness/2026-07-31-1300-doubly-noncentral-t-pdf-cancellation-x-mu-negative.md (#1235).
test/precision-continuous.js's NoncentralChi2([268, 64]) quantile round-trip gate (qtol) was too tight at 1e-13, consistently failing (measured ~1.015e-13-1.05e-13) under full-parallel-suite npm test runs while passing in isolation — the same JIT-order-dependent floating-point summation-order sensitivity already documented for sibling marcumQ-adjacent groups. qtol is now 5e-13, matching the established tolerance already used for NoncentralChi2([5, 58]), NoncentralChi2([5, 62]), NoncentralChi2([270, 64]), and NoncentralChi([5, 7.5]); no reference value or pdf/cdf tolerance changed.
ran.process.AR1.variance(t) lost all significance for near-unit-root phi (phi² just outside the existing 1e-14 special-case band) combined with small fractional t (< 0.1): Math.pow(phi2, t) rounds to exactly 1.0 in double precision there, so 1 - Math.pow(phi2, t) evaluated to exactly 0 instead of the true small positive variance — e.g. variance(1e-6) returned 0 instead of ~1e-6 for phi2 = 1 - 2e-14. A numerical sweep (phi2 deltas 1e-14–1e-1, t up to 1e300) found this was the only real failure mode — the originally-suspected large-t scenario (negative/NaN variance) never occurred. Fixed by replacing 1 - Math.pow(phi2, t) with the cancellation-safe -Math.expm1(t * Math.log(phi2)), matching the existing expm1/log1p idiom used elsewhere in the codebase (e.g. ran.dist.Pareto, ran.dist.Weibull); the 1e-14 special case is unchanged (still required at phi2 === 1 to avoid 0/0) (#1243). ran.process.AR1.covariogram(s, t) carried a second, independent copy of the same 1 - Math.pow(phi2, min(s, t)) expression and was left unfixed by that pass; it failed identically, and was caught by the sweep run for #1244. Since Cov(X_t, X_t) = Var(X_t) by definition, the two methods openly disagreed: for phi2 = 1 - 2e-14, covariogram(1e-6, 1e-6) returned exactly 0 against variance(1e-6)'s correct ~1e-6 (100% error), and covariogram(0.01, 0.01) was off by 11%. The same -Math.expm1(...) reformulation is now applied there, and covariogram() gains the min(s, t) === 0 fast path variance() already had at t === 0 — without it the reformulation would have turned 0 * Math.log(phi2) into NaN whenever phi2 underflows to 0 or overflows to Infinity, which the old Math.pow(phi2, 0) === 1 identity had made safe (covariogram(0, 3) for phi = 1e200 returned NaN even before this change, since Infinity * -0 is already NaN). Trade-off, stated plainly: -expm1(n·log(x)) amplifies log's rounding error by n, so for a strongly explosive process at large min(s, t) the new form is less accurate than Math.pow was — e.g. phi = 1.5, s = t = 200 moves from 3.4e-17 to 9.6e-15 relative error against an mpmath mp.dps=60 reference. That is a deliberate exchange of ~2 digits in a regime whose value is already ~1e70 and diverging, for the elimination of a 100% error near the unit root; variance() has made the identical trade since #1243, and keeping both methods on one formulation is what makes the Cov(t,t) = Var(t) identity hold exactly.
ran.special.marcumQ/marcumP's _fc(nu, z) (the modified-Lentz continued fraction for I_nu(z)/I_{nu-1}(z), seeding the mu < 135 transition-band recurrence) silently returned an unconverged value once z grew past roughly 250-300, because its loop was capped at the shared MAX_ITER = 100 with no convergence check on exit — the required depth scales as ~6.2*sqrt(z), not a constant, so e.g. NoncentralChi2(200, 2000).cdf(2080) (z ≈ 2038, needing 189 iterations) was off by 5.1e-08 relative instead of the library's usual ~1e-14 floor, and the mu = 134/mu = 135 transition-band boundary carried a six-orders-of-magnitude accuracy discontinuity (_largeMu, used for mu >= 135, never calls _fc and was unaffected). _fc now computes a regime-aware local iteration budget (Math.max(MAX_ITER, Math.ceil(7 * Math.sqrt(z)) + 20), stress-tested across nu in (0, 135) and z up to 1e5 with zero non-convergent cases) and throws if that budget is ever exceeded, rather than returning the unconverged value — matching the existing "throw on exceeded iteration budget" convention in src/algorithms/rejection.js. NoncentralChi2(200, 2000).cdf(2080) now matches the mpmath (mp.dps=50) reference to 1.5e-12 relative, the same value an effectively-uncapped _fc produces, confirming the residual is _recurrence's own pre-existing seed/amplification floor rather than further _fc truncation. Adds the large-x recurrence-regime precision-gate set (NoncentralChi2[76, 692]) that #1190/#1143 deliberately withheld until this fix landed. See solutions/special-functions/2026-08-02-1200-marcum-fc-slow-convergence.md (#1286). scripts/precision-refs-continuous.py's existing_groups() (dev-only tooling) separately gains a fix for a different pre-existing parsing failure surfaced while regenerating this gate: its brace-depth scan tracked every {/} character in the file including ones inside // comments (e.g. a comment referencing the JS snippet { useFloor: false }), so a balanced brace pair inside a comment was misread as a whole REFS group, corrupting every span parsed after it. Comment-only lines are now blanked out before the structural scan runs.
ran.special.besselISpherical's _hi(n, x) continued-fraction helper (whose iteration budget #1292 widened to Math.max(MAX_ITER, Math.ceil(7 * Math.sqrt(x)) + 20), described below) had no convergence check on loop exit, unlike marcum-q.js's sibling _fc, which already throws via _assertFcConverged (#1286) instead of returning an unconverged value silently. _hi now gains an equivalent _assertHiConverged check, thrown when |del/h| > EPS after the loop exits, matching the "throw on exceeded iteration budget" convention in src/algorithms/rejection.js. No valid distribution parameterization in this codebase (NoncentralChi, NoncentralChi2, the only two distributions that call _hi, always with a non-negative sqrt(lambda*x)-derived argument) is known to reach non-convergence within the existing regime-aware budget; the check is a defensive availability guard against an extreme, currently-unreached caller-supplied argument (e.g. NoncentralChi2(...).pdf(1e12)) rather than a fix for an observed wrong value. _hi and _fc's shared "budget grows with sqrt(argument), no fixed upper ceiling" design is now documented as an explicit accepted trade-off in both files, made safe by each throwing on non-convergence instead of truncating silently or running unbounded (#1311). besselISpherical/besselISphericalExpScaled (src/special/bessel.js) and Distribution.prototype.pdf() (src/dist/_distribution.js) now carry @throws JSDoc documenting this exception where it is actually reachable by a caller — pdf() carries a single precise tag naming its narrow scope (NoncentralChi/NoncentralChi2, odd k only) rather than repeating it across hazard()/lnPdf()/lnL()/aic()/bic(), which all call pdf() and inherit the same documented exception. See ADR-0049, which reconciles throw (over NaN) against decisions/0015-return-value-and-error-conventions.md: the continued fraction's true value is finite and well-defined for any valid argument, so a non-convergence is an algorithmic budget failure, not the mathematically-indeterminate case NaN is reserved for (#1326).
ran.dist.NoncentralChi2.pdf(x)/ran.dist.NoncentralChi.pdf(x) returned NaN once lambda * x (or lambda^2 * x^2 for NoncentralChi) grew past roughly 5e5 — e.g. NoncentralChi2(100, 900).pdf(1000) (an ordinary parameterization evaluated near its own mean) — because _pdf combined a log-space prefactor (exp(-0.5*(x+lambda)), which underflows to exactly 0 in this regime) with a linear-space Bessel factor (besselI/besselISpherical evaluated at sqrt(lambda*x), which overflows past Number.MAX_VALUE once its argument exceeds ~710-720): 0 * Infinity is NaN even though the true density is an ordinary, representable double. The same class of bug as #1075's DoublyNoncentralBeta overflow. src/special/bessel.js gains two exponentially-scaled accessors — besselIExpScaled(n, x) = exp(-|x|) * I_n(x) (reusing _besselIBackward's existing internal ratio, which is already this exact quantity before its final * exp(x) step) and besselISphericalExpScaled(n, x) = exp(-x) * i_n(x) for x >= 0 (a Wronskian rebuilt from _knRaw's un-normalized upward-recurrence values instead of _kn's exp(-x)-scaled ones, so the exponent never has to be materialized and immediately inverted back out) — and both _pdf methods now fold the Bessel argument's exponent into the existing log-space prefactor before exponentiating, relying on the identity -0.5*(x+lambda) + sqrt(lambda*x) = -0.5*(sqrt(x)-sqrt(lambda))^2 <= 0 (AM-GM) to keep the combined exponent always finite. Reaching this newly-representable regime also exposed a second, previously-unreachable defect: _hi's continued fraction (used by besselISpherical's Wronskian branch) shares the same MAX_ITER = 100 cap _fc was fixed for above (#1286) and silently under-converged past x ~ 250 — it now uses the identical regime-aware budget, Math.max(MAX_ITER, Math.ceil(7 * Math.sqrt(x)) + 20). NoncentralChi2(100, 900).pdf(1000), NoncentralChi2(200, 2000).pdf(2080), and NoncentralChi(200, 44.7).pdf(45.6) now return finite values matching an independent Poisson-mixture cross-check and an mpmath (mp.dps=50) reference; existing small-lambda precision-gate values are unchanged (#1292). The same defect shape hit ran.dist.Skellam.pdf(x): _pdf combined a separate expNeg = exp(-mu1-mu2) prefactor (underflowing to exactly 0 once mu1+mu2 > ~745) with besselI(|x|, twoSqrtProd) (overflowing to Infinity once twoSqrtProd = 2*sqrt(mu1*mu2) > ~709-720, a lower threshold than mu1+mu2 itself) — e.g. Skellam(360, 360).pdf(0) returned Infinity (only the Bessel factor had overflowed) and Skellam(400, 400).pdf(0) returned NaN (0 * Infinity, both factors past their threshold). The constructor's speed-up constants now precompute expNegScaled = exp(-mu1-mu2+twoSqrtProd), which stays in (0, 1] since -mu1-mu2+twoSqrtProd = -(sqrt(mu1)-sqrt(mu2))^2 <= 0 always, and _pdf combines it with besselIExpScaled(|x|, twoSqrtProd) (added by #1292) instead of the unscaled besselI. Skellam(360, 360).pdf(0), Skellam(400, 400).pdf(0), and Skellam(2000, 2000).pdf(0) now return finite values matching mpmath (mp.dps=50) to the project's 1e-14 precision-gate tolerance; existing small-mu precision-gate values are unchanged (#1309). The same defect shape also hit ran.dist.VonMises(mu, kappa).pdf(x)/.cdf(x), NaN for kappa past roughly 710-720 — e.g. VonMises(0, 720).pdf(0), VonMises(0, 800).pdf(0.001), VonMises(0, 800).cdf(0.5) — because exp(kappa*cos(x-mu)) and besselI(0,kappa) both independently overflow to Infinity there, and _cdf's Fourier series hit the identical Infinity/Infinity in every term. _pdf is rewritten as exp(kappa*(cos(x-mu)-1)) / (2*pi*besselIExpScaled(0,kappa)), whose numerator exponent is bounded <= 0 by cos(x-mu) <= 1; _cdf's series envelope substitutes besselIExpScaled(i,kappa)/(besselIExpScaled(0,kappa)*i) for the old besselI(i,kappa)/(besselI0Kappa*i), an algebraically exact substitution since the shared exp(-kappa) factor cancels, preserving the existing oscillating-term-safe convergence check (solutions/correctness/2026-07-26-1339-vonmises-cdf-oscillating-term-premature-convergence.md) unchanged. _cdf's return is now also clamped to [0, 1] (Math.max(0, Math.min(1, ...)), the same guard already used in noncentral-beta.js), since 0.5*(1+dx/pi) + sum/pi cancels two O(1) terms and can round a few ULPs outside [0, 1] for x far from mu — a pre-existing characteristic of that formula, only reachable now that large kappa no longer immediately overflows to NaN (#1308). That cancellation is now fixed by #1320: _cdf no longer computes 0.5*(1+dx/pi) + sum/pi at all, replacing the Fourier series entirely with direct tanhSinh quadrature of the already cancellation-free _pdf over the tail interval (using the pdf(mu+t) = pdf(mu-t) symmetry to always integrate on the side away from the density's peak at mu), which is monotonic by construction and accurate arbitrarily deep into the tail instead of merely clamped to [0, 1] — e.g. VonMises(0, 730).cdf(-0.357) now returns ~4.29e-22 (matching an mpmath mp.dps=50 reference) instead of the previous 2.78e-16 cancellation noise that made .cdf(-0.355) come out below .cdf(-0.357) despite -0.355 > -0.357. Existing pdf/cdf precision-gate values for kappa in {0.5, 1, 2, 9, 11, 1000, 1500, 2000} are unchanged within their existing tolerances.

[1.31.0] - 2026-07-20

Added

ran.mc.MCMC.state() now round-trips the sampler's complete PRNG-and-adaptation state, making a resumed sampler's subsequent draws bit-for-bit identical to an uninterrupted run instead of merely statistically equivalent (#1033, ADR-0035). state() gains a top-level prng key (the Xoshiro128+ stream position, restored by the constructor via Xoshiro128p.save()/.load(), mirroring ran.dist.Distribution.save()/.load()'s existing prngState precedent), and every subclass's _internal() now also serializes its own proposal/momentum generator (ran.mc.RWM, ran.mc.AdaptiveMetropolis, ran.mc.HMC, ran.mc.NUTS, ran.mc.MALA) and adaptation-batch accumulators — Robbins-Monro counters (RWM, MALA, ran.mc.Slice), the running covariance accumulator (AdaptiveMetropolis), and the dual-averaging and mass-matrix accumulators (HMC, NUTS) — superseding the prior "serialize effective state only" convention (ADR-0020 §2, ADR-0021, ADR-0029) for these specific fields. ran.mc.Gibbs needed no changes: restoring the base class's PRNG alone is sufficient since Gibbs has no subclass-owned generator or adaptation state. Old snapshots (missing the new fields) continue to construct valid instances with the pre-#1033 behavior — additive, non-breaking. Two scope boundaries are documented rather than silently unmet: samplingRate (thinning) is not guaranteed to reproduce exactly across a mid-warm-up resume for any subclass, since it depends on the base-class autocorrelation accumulator, which stays intentionally unserialized (ADR-0023); and RWM's per-dimension proposal scale (_base) is not guaranteed bit-for-bit reproducible across a resume that lands strictly mid-batch-window, because _refreshBase() depends on the same never-serialized base-class Welford accumulator — every other subclass's adaptation state is unaffected by this gap.
ran.mc.NUTS now reports sampler-health diagnostics, matching the per-iteration divergent/maxTreeDepthReached signals Stan/PyMC/NumPyro expose. Every iterate() result carries a divergent boolean (a leapfrog leaf whose Hamiltonian drifted past the energy-divergence threshold — step size too large or target geometry too extreme) and a maxDepthHit boolean (the doubling tree saturated MAX_TREE_DEPTH without a U-turn — step size too small), and two aggregate accessors, divergenceCount() and maxDepthCount(), report the per-sampling-phase totals. The counts ride the same accumulator lifecycle as ar() (reset at construction and at each sample() start, so a read afterwards reflects the sampling phase only); a well-behaved run reports both as 0. Diagnostic-only — no change to sampling behavior (#1037, ADR-0048).
ran.mc.NUTS now supports Euclidean metric (mass matrix) adaptation via config.metric, matching ran.mc.HMC: 'diag' (default) adapts a per-dimension variance and 'dense' adapts the full covariance matrix (factored via Matrix.ldl()) during warm-up. Momentum is resampled from N(0, M) instead of a standard Normal, the leapfrog integrator and kinetic energy apply the metric, and the no-U-turn criterion is evaluated on the velocity M⁻¹r; the adapted metric round-trips through state()/_internal(). This removes the previous capability regression where poorly-scaled or correlated targets mixed better under HMC than NUTS (#1035, ADR-0034).
All 11 ran.mc samplers (AdaptiveMetropolis, ARS, gelmanRubin, Gibbs, HMC, MALA, NUTS, ParallelTempering, runChains, RWM, Slice) are now available as tree-shakeable subpath imports under a dedicated mc namespace (import RWM from 'ranjs/mc/rwm', import gelmanRubin from 'ranjs/mc/gelman-rubin'), matching the per-distribution (ranjs/dist/<name>) and per-process (ranjs/process/<name>) subpath export patterns. Each resolves to a self-contained ESM bundle at dist/mc/<name>.esm.js (#1036).
ran.mc.AdaptiveMetropolis(logDensity, config, initialState): full-covariance adaptive Metropolis sampler (Haario, Saksman & Tamminen, 2001). Adapts the joint proposal covariance Sigma_proposal = (2.38^2 / dim) * Cov(x) + epsilon * I from the chain's own history during warm-up via an online covariance accumulator and Matrix.ldl(), then freezes the covariance for the sampling phase. Mixes substantially better than RWM's diagonal-only adaptation for correlated multi-dimensional targets (#823).
ran.mc namespace (RWM, gelmanRubin) is now exported from the library's entry point, wiring it up to ran.mc after it was inadvertently left unexported during PR #615's cleanup (#617).
seed(value) method on ran.mc.MCMC (and ran.mc.RWM, which additionally reseeds its internal proposal distribution) for deterministic, reproducible sampling. Internally, both classes now use a per-instance Xoshiro128p PRNG instead of the shared module-level generator, so seeding a sampler no longer affects unrelated code sharing that singleton. If the initial position was not explicitly supplied, seed() also redraws it from the newly seeded generator so that .seed(s).sample(n) is fully reproducible (#912).
ran.mc.runChains(logDensity, config, options): runs multiple independently-seeded RWM chains and computes the gelmanRubin() diagnostic across them in one call — the recommended workflow (ADR-0024) for gating MCMC convergence, since no signal computable from a single chain can distinguish "converged" from "stuck". Defaults to 2 chains seeded [1, 2]; options.chains, options.warmUpBatches, options.sampleSize, options.seeds, and options.maxLength are all configurable. Returns { samples, rhat } (#935).
ran.mc.Gibbs(conditionals, config, initialState): component-wise (systematic-scan) Gibbs sampler, implemented as an MCMC subclass. Cycles through each dimension in order, replacing it with a draw from the caller-supplied full conditional given the current state. Every draw comes directly from the exact conditional, so there is no accept/reject step and ar() is always 1.0 (#821).
ran.mc.HMC(logDensity, gradLogDensity, config, initialState): Hamiltonian Monte Carlo sampler, implemented as an MCMC subclass. Uses a leapfrog integrator over config.pathLength steps of size config.stepSize to propose distant moves along Hamiltonian trajectories, with momenta resampled from N(0, I) each iteration and Metropolis accept/reject on the augmented (position, momentum) system. Step size is adapted during warm-up via Robbins-Monro dual averaging (Hoffman & Gelman 2014) toward a target acceptance probability, and jittered multiplicatively (ε ~ Uniform(0.9ε, 1.1ε)) each iteration to avoid periodicity artifacts (#824). Now also supports Euclidean metric (mass matrix) adaptation via config.metric: 'diag' (default) estimates a per-dimension variance online during warm-up so the sampler mixes efficiently on targets whose parameters span very different scales; 'dense' (opt-in) estimates the full covariance matrix via an online accumulator, regularized and factored through Matrix.ldl(), so the sampler also compensates for correlated parameters. The adapted metric round-trips through state()/_internal() alongside stepSize/pathLength (#826).
ran.mc.MALA({ logDensity, gradLogDensity, config, initialState }): Metropolis-Adjusted Langevin Algorithm sampler, implemented as an MCMC subclass. Proposes a single gradient-informed Langevin step per iteration (x' = x + (stepSize² / 2) * ∇log p(x) + stepSize * z, z ~ N(0, I)) and applies a Metropolis-Hastings correction for the proposal's asymmetry, making the chain exact. Step size is adapted during warm-up via batch Robbins-Monro (Roberts & Rosenthal 2009), the same scheme RWM uses, toward the MALA-optimal 0.574 acceptance rate (Roberts & Rosenthal 1998) (#828). Ships with an options-object-only constructor from its first release — there is no positional form and no deprecation warning, unlike RWM/Slice/AdaptiveMetropolis (#970).
ran.mc.NUTS({ logDensity, gradLogDensity, config, initialState }): No-U-Turn Sampler, implemented as an MCMC subclass using the identity-mass leapfrog integrator extracted to src/mc/_leapfrog.js, combined with Hoffman & Gelman's (2014) doubling-tree algorithm. Automatically tunes the trajectory length each iteration by recursively extending a leapfrog trajectory forward or backward in a random direction until the trajectory's outer endpoints start turning back toward each other (the U-turn criterion) or a maximum tree depth is reached, selecting the transition via slice sampling over the tree's valid states — eliminating the need to hand-tune pathLength. Step size is adapted during warm-up via the same Robbins-Monro dual averaging as HMC, driven by the tree-averaged acceptance statistic, toward a target acceptance probability of 0.8 (#825). Ships with the options-object constructor form from its first release — no positional form, no deprecation cycle (#972).
ran.mc.ARS(logDensity, support, derivative): Gilks-Wild (1992) Adaptive Rejection Sampling for univariate log-concave densities on a finite support bracket. Builds a piecewise-exponential upper envelope (and a secant lower "squeeze" hull) from tangent lines to the log-density, adaptively tightening on every rejection so acceptance probability increases monotonically; throws Error for non-log-concave targets. Unlike the rest of ran.mc, it is not an MCMC subclass — it produces exact i.i.d. draws directly, with no warm-up or accept/reject Markov-chain machinery (#820).
ran.mc.Slice(logDensity, config, initialState): coordinate-wise slice sampler (Neal 2003), implemented as an MCMC subclass. Requires only logDensity — no proposal tuning, no gradient. Each dimension is updated via stepping-out and shrinkage; the interval width w (default 1.0) is the only tunable parameter and is adapted per dimension during warm-up. Every sweep produces an accepted draw, so ar() is always 1.0. A prior, non-functional slice.js (100% commented out, never wired into the base class) was removed as dead code in PR #615; this is a fresh implementation (#822).
ran.mc.ParallelTempering(logDensity, options): Parallel Tempering / Replica Exchange MCMC (Geyer 1991) for multimodal targets. Runs N independent replica samplers (default RWM, or a caller-supplied sampler factory) at inverse temperatures beta_1 = 1 > beta_2 > ... > beta_n — an explicit options.temperatures array, or an auto-generated geometric ladder from options.nReplicas/options.tempMax. warmUp() tunes every replica independently; sample() runs all replicas in lockstep and, after each thinned step, proposes a swap between one alternating-parity set of adjacent replica pairs, accepted with probability min(1, exp((beta_i - beta_j)(log p(x_j) - log p(x_i)))) per detailed balance on the joint replica distribution, returning the cold (beta = 1) replica's samples; swapRate() reports the accepted/attempted fraction per adjacent pair. Unlike the rest of ran.mc, it is not an MCMC subclass — it coordinates an array of replicas rather than driving a single chain, and does not support state()/resumption (ADR-0028, #830).
ess() method on ran.mc.MCMC: computes the Effective Sample Size per dimension using Geyer's initial positive monotone sequence estimator (IPSM), N / (-1 + 2 * sum_m Gamma_m), where Gamma_m = rho[2m] + rho[2m+1] pairs consecutive lags starting at lag 0 (from the existing ac() accumulators, so the first pair always includes rho[0] = 1) and is clamped to be non-increasing, summed until the first pair whose clamped value is not positive (falling back to ess = N if even the first pair is non-positive). A fully stuck (zero-variance) chain, where ac() returns NaN at every lag, reports ess = 1 rather than saturating to N. Reads directly from the online accumulators already backing ac() and statistics() — no new accumulator state (#827, #975).

Changed

ran.mc.runChains() is generalized to drive any ran.mc.MCMC subclass instead of hardcoding RWM: the new signature is runChains(Sampler, samplerOptions, runOptions), where samplerOptions is forwarded verbatim to new Sampler(samplerOptions) for every chain — the same options-object shape that sampler's own constructor accepts ({logDensity, config, initialState} for RWM/AdaptiveMetropolis/Slice, {logDensity, gradLogDensity, config, initialState} for HMC/MALA/NUTS, {conditionals, config, initialState} for Gibbs). runOptions keeps the previous {chains, warmUpBatches, sampleSize, seeds, maxLength} shape. ran.mc.gelmanRubin() is unaffected, since it only ever consumed the returned per-chain sample arrays, never the sampler that produced them. See ADR-0033 (#967).
ran.mc.RWM now uses a consistent joint diagonal adaptive-Metropolis proposal in both warm-up and sampling instead of tuning per-component (Metropolis-within-Gibbs, 0.44 target) during warm-up and switching to joint proposals for sampling. Warm-up adapts a single global step scale via batch Robbins-Monro toward the optimal acceptance rate (0.44 for dim = 1, 0.234 for dim > 1) and tracks per-component scales from the running marginal standard deviations, so the proposal that is tuned is the proposal that samples. Behavior for dim = 1 is unchanged (the two schemes coincide); multi-dimensional targets are now correctly tuned. See ADR-0022.
ran.mc.MCMC.ar() now reports the acceptance rate over a sliding window of the most recent config.arWindow iterations (default 1000) instead of the cumulative rate since the last reset, so mid-warmUp() reads aren't dragged down by early untuned batches. During the partial-fill phase (fewer than arWindow iterations since reset) the value is unchanged from before. See ADR-0021 (#920, #926).
Code Health improved across three files by extracting shared/named helpers: test/special.js (8.28 → 9.09, shared check/checkBesselIdentity/checkF11Recurrence helpers), src/special/marcum-q.js (8.67 → 10.0, _expansionSum/_transitionBand/_initPhi helpers eliminating three Complex Method smells), and test/dist.js (8.76 → 9.09, assertFitSpec/assertParamRecovery helpers eliminating a Complex Method and Excess Arguments smell).
ran.mc.HMC's class-level documentation and pathLength parameter docs now disclose that a fixed pathLength can still produce genuine resonance-driven negative lag-1 autocorrelation at certain target correlations, even with the existing per-iteration stepSize jitter — confirmed empirically (an investigation swept both target correlation and pathLength, finding resonance bands as narrow as 2-3 integer pathLength steps that a ±10%-scale jitter cannot reliably escape) — and point affected users to ran.mc.NUTS, which adapts trajectory length automatically. No behavior change (#1005).

Deprecated

ran.mc.ParallelTempering's positional constructor form new ParallelTempering(logDensity, options) is deprecated in favor of the options-object form new ParallelTempering({ logDensity, ...options }), bringing it in line with every other ran.mc sampler and coordinator (RWM, AdaptiveMetropolis, Slice, HMC, MALA, NUTS, Gibbs per ADR-0030; ARS per ADR-0031) and removing the last positional-constructor wart in ran.mc. The positional form still constructs and samples correctly but emits a one-time console.warn on first use; it will be removed in v1.32.0 (#1034).
ran.process.PoissonProcess and ran.process.CompoundPoissonProcess are renamed to ran.process.Poisson and ran.process.CompoundPoisson: no ran.process.Process subclass name should redundantly repeat "Process" (every other subclass — BrownianMotion, OrnsteinUhlenbeck, AR1, RandomWalk, etc. — already follows this). The old names still construct and behave identically (PoissonProcess extends Poisson, CompoundPoissonProcess extends CompoundPoisson) but emit a console.warn on construction and will be removed in v1.33.0 (ADR-0041).

Removed

ran.mc.RWM's, ran.mc.AdaptiveMetropolis's, ran.mc.Slice's, ran.mc.HMC's, and ran.mc.Gibbs's deprecated positional constructor arguments (e.g. new RWM(logDensity, config, initialState)) are removed; only the options-object form (new RWM({ logDensity, config, initialState })) remains. ran.mc.runChains()'s deprecated legacy call form, runChains(logDensity, config, options), is likewise removed; only the generalized runChains(Sampler, samplerOptions, runOptions) form remains. This closes out the deprecation cycle introduced in #962–#967 (ADR-0030, ADR-0031, ADR-0033). Note on the deprecation cycle: CLAUDE.md's normal deprecation-cycle rule requires a released minor version containing the warning to ship and hold for a full release before the removal lands; here the #962–#967 ### Deprecated entry never left the [Unreleased] section of this changelog (v1.30.0 predates it), so the removal is landing without that hold, at explicit maintainer request overriding the standard process (#968). No published version of ranjs ever carried the positional forms as deprecated-but-working, so no downstream user is exposed to a behavior change without warning — the positional forms and their console.warn deprecation notices are simply gone, as if they had never been introduced.

Fixed

ran.mc.RWM, ran.mc.AdaptiveMetropolis, ran.mc.Slice, ran.mc.HMC, and ran.mc.Gibbs constructors now throw a clear, class-specific Error (e.g. "RWM: constructor requires an options object: new RWM({ logDensity, config, initialState })") when called with null, any other non-plain-object argument, or no argument at all, instead of either a generic, engine-dependent TypeError: Cannot destructure property '...' of 'null' as it is not an object. or, for a zero-argument call, silently constructing an unusable instance that only fails later inside _iter(). These five constructors became options-object-only in #968 but never received the guard ran.mc.MALA/ran.mc.NUTS already had (#970, #972), so they regressed to the confusing destructuring error MALA and NUTS were already fixed to avoid (#1029).
ran.mc.MCMC (and all subclasses, e.g. ran.mc.RWM) now reject config.dim above 10000, config.maxLag above 10000, config.arWindow above 10000, and any individually-valid dim/maxLag combination whose combined accumulator footprint (dim * maxLag * 16 bytes) exceeds 100MB, throwing a clear Error instead of allocating oversized arrays and crashing the process with an out-of-memory error (#916, #922, #928). ran.mc.HMC now rejects config.pathLength above 1024 (2^10, matching the NUTS sampler's own literature-derived MAX_TREE_DEPTH ceiling — the Stan/PyMC/NumPyro default), throwing instead of letting warmUp()/sample() hang indefinitely on the per-iteration leapfrog cost of an unreasonably large path length (#947, #989).
ran.mc.MCMC warm-up thinning no longer inverts for slow-mixing chains: when a dimension's autocorrelation never decays to ≤ 0.05 within maxLag, _thinningLag() now falls back to the largest measured lag instead of reporting 0. Previously a chain that mixed slower than maxLag could resolve was treated as already-decorrelated, driving samplingRate down toward 1 and under-thinning sample() — the opposite of the intended "slowest-mixing dimension wins" rule (ADR-0020 §3).
ran.mc.MCMC.warmUp(progress, maxBatches) now runs exactly maxBatches batches (was maxBatches + 1 due to a batch <= maxBatches loop bound) and reports 100 at completion instead of firing a redundant 0% callback at the start.
ran.mc.MCMC.sample(progress, size) now reports each integer percentage of progress exactly once. Previously the i % (iMax/100) check used a fractional modulus whenever the total iteration count was not a multiple of 100, silently skipping most progress callbacks.
ran.mc.Gibbs's conditionals now receive the sampler's own PRNG as a second argument (conditionals[d](x, rng)), so seed() can produce reproducible chains for conditionals that draw their randomness from rng.next() instead of an independently-seeded generator. Previously Gibbs._iter() never read this.r, so gibbs.seed(42).sample(null, N) silently failed to reproduce, violating the contract documented on MCMC.seed() (ADR-0026, #938).
ran.mc.RWM, ran.mc.AdaptiveMetropolis, and ran.mc.Gibbs constructors now have a dedicated JSDoc @param block directly on constructor(), so tsc's generated .d.ts resolves the true parameter types (Function/Function[]) instead of any; Gibbs previously had no constructor signature in the generated declaration at all (#944).
ran.mc.ARS's hull segment-mass (_build), envelope inverse-CDF (_sampleEnvelope), and tangent-intersection breakpoint (_tangentIntersection) formulas now treat a tangent slope, or a difference between two tangent slopes, as zero once it falls below a Math.cbrt(EPS)-scaled tolerance (matching the noise floor already used elsewhere in the same file for finite-difference-derived slopes), instead of only below raw Number.EPSILON. A small-but-nonzero slope (or near-indistinguishable slope pair) previously fell through to a general-case formula — differencing two nearly-equal exponentials in _build/_sampleEnvelope, or dividing by a near-cancelled denominator in _tangentIntersection — a catastrophic-cancellation pattern that could place a hull breakpoint far outside its valid bracket (#941, #957).
The docs build's assembleLinks() (docs/src/desc-parser.js) no longer skips every second {@link} construct in a JSDoc paragraph: an off-by-one advanced the loop index by 3 instead of 2 after each converted (text, link) pair, silently dropping the next pair instead of converting it to a hyperlink. src/mc/adaptive-metropolis.js's previously-masked bare {@link ran.mc.RWM} reference is now written in the codebase's standard bracketed [RWM]{@link ran.mc.RWM} form so it renders as a working link instead of tripping the #980 bare-link guard (#997).
ran.mc.RWM, ran.mc.AdaptiveMetropolis, ran.mc.HMC, ran.mc.MALA, and ran.mc.NUTS no longer alias their proposal/momentum generator with their accept/reject generator after seed(). MCMC._reseedCachedLogDensity() seeded the subclass-owned _q generator with the same raw value passed to this.r; since Xoshiro128p.seed() is a pure function of its argument with no per-instance salt, both generators produced byte-identical streams, so the Metropolis acceptance uniform was a value already consumed to build an earlier proposal — violating the independence the MH ratio assumes. _q is now seeded from a derived value (`${value}-q`), mirroring ParallelTempering's per-replica seeding. Reproducibility is preserved (deterministic derivation).
ran.mc.HMC and ran.mc.NUTS no longer permanently freeze when the caller's gradient returns NaN at a visited state (e.g. a hand-written gradient that yields NaN near a support boundary instead of a rigorous -Infinity). A non-finite acceptance statistic previously flowed unchecked through the Robbins-Monro dual-averaging recursion into _daHbar → _daLogEpsBar → stepSize, becoming a sticky NaN that silently stopped the sampler from ever moving again. _adjust() now treats a non-finite acceptance statistic as a fully-rejected (divergent) step, driving the step size down so warm-up recovers, matching Stan's divergent-proposal handling.
ran.mc.Slice now throws for a w (in initialState.internal.w) that is neither a number nor an array (e.g. a string, boolean, object, or null). Such values were silently coerced to the 1.0 default before validation ran, so a documented-parameter type error passed unchecked instead of failing fast per the library's return-value conventions.
ran.mc.Slice now throws (rather than hanging indefinitely) when logDensity returns NaN at the current point: logY then becomes NaN, lnp(candidate) > logY is always false, and the shrinkage loop narrows forever without accepting. _shrink() is now bounded by a MAX_SHRINK cap — the shrink analogue of the existing w: Infinity stepping-out guard.
ran.mc.AdaptiveMetropolis's proposal-covariance regularization now scales the epsilon term by s_d, matching Haario, Saksman & Tamminen (2001)'s C_n = s_d * (Cov(x) + epsilon * I). Previously the fixed epsilon = 1e-6 sat outside the s_d = 2.38^2/dim factor (s_d * Cov(x) + epsilon * I), so the regularization floor grew relatively more influential as dimension increased — the opposite of the reference formula's proportional shrinkage.
ran.mc.gelmanRubin() now throws when the supplied chains do not all have the same length. _gri used chains[0].length as the sample-variance divisor for every chain, so an unequal-length chain (reachable via direct calls, though never via runChains) silently read past its end (undefined → NaN) and mismatched its divisor, producing a wrong or NaN R-hat instead of an error.

Security

Remediated npm audit findings (#960): @babel/core patched to a version above the arbitrary-file-read range (GHSA-4x5r-pxfx-6jf8) via npm audit fix; nyc bumped ^15.1.0 → ^18.0.0, which pulls a fixed istanbul-lib-processinfo that no longer depends on the vulnerable uuid (GHSA-w5hq-g745-h8pq) — verified against the full test suite, including its coverage-threshold gate; serialize-javascript pinned to ^7.0.7 via a new overrides entry to close mocha's transitive RCE/DoS vulnerabilities (GHSA-5c6j-r48x-rmvq, GHSA-qj8w-gfj5-8c6v), since mocha's own package.json range (^6.0.2) predates the fix even on its latest release. All three changes are devDependency-only; none affect src/ or the published package. Accepted risk, documented and left unresolved because no upstream fix exists: documentation@14.0.3 (latest release) bundles vue-template-compiler@2.7.16 (latest ever published, Vue 2 tooling is EOL) which has an XSS advisory (GHSA-g3ch-rx76-35fx) exploitable only via untrusted template compilation — not applicable here, since npm run docs only compiles maintainer-authored templates locally; the mathjax-node-page toolchain (mathjax, mathjax-node, jsdom, request, request-promise-core/-native, form-data, qs, tough-cookie, nested uuid, yargs/yargs-parser) is abandoned upstream (last publish 2022, itself depending on the long-deprecated request library), so npm audit fix --force's suggested resolution is an older mathjax-node-page release that carries the identical vulnerable subtree — it doesn't fix anything. Both chains are used exclusively by docs/index.js for local, maintainer-invoked API doc generation (npm run docs); docs/ is excluded from the package's files field, and neither dependency runs during npm test, npm run build, or at library runtime. Since nyc@18 declares engines.node: "20 || >=22", the CI test-job matrix (.github/workflows/ci.yml) drops Node 18 (now [20, 22]); this only affects the project's own contributor/CI tooling — the published dist/ bundle carries no Node version requirement.

[1.30.0] - 2026-07-09

Added

ran.process.CompoundPoissonProcess(jumpDist, lambda, dt): compound Poisson process accumulating random-magnitude jumps at rate λ. At each step draws K ~ Poisson(λ·Δt) arrivals and sums K independent samples from the supplied ran.dist Distribution instance. Exposes mean(t) (λ·t·E[J]), variance(t) (λ·t·E[J²]), and covariogram(s,t) (λ·E[J²]·min(s,t)) using the jump distribution's analytical moments; pdf is not implemented as it has no closed form for arbitrary jump distributions (#860).
ran.process.CoxIngersollRoss(kappa, theta, sigma, dt): Cox-Ingersoll-Ross mean-reverting process for positive-valued dynamics (interest rates, variance). Uses an Euler-Maruyama step with reflection max(X_n, 0) inside the noise term to prevent negative states. Warns (but does not throw) when the Feller condition 2κθ > σ² is not met. Exposes mean(t), variance(t), pdf(x,t) (Gamma marginal for x0=0), and covariogram(s,t) with closed-form analytical values (#858).
ran.process.RandomWalk(p): discrete-time random walk on the integers; at each step the state moves by +1 with probability p (default 0.5, symmetric) or −1 with probability 1 − p. Exposes mean(t) (t·(2p−1)), variance(t) (4p(1−p)·t), pdf(x,t) (exact binomial PMF), and covariogram(s,t) (4p(1−p)·min(s,t)) (#859).
ran.process.AR1(phi, sigma): first-order autoregressive process with update rule X_{n+1} = φ·X_n + σ·Z where Z ~ N(0,1). For |φ| < 1 the process is stationary with marginal distribution N(0, σ²/(1−φ²)); for |φ| ≥ 1 the process is non-stationary and a console.warn is emitted (no error thrown). Exposes mean(t) (always 0), variance(t), pdf(x,t), and covariogram(s,t) with closed-form analytical values (#857).
pdf(x, t) method on ran.process.BrownianMotion, ran.process.OrnsteinUhlenbeck, ran.process.GeometricBrownianMotion, ran.process.BrownianBridge, and ran.process.PoissonProcess: returns the marginal density/mass of the process at state x and time t. BM and OU use the Normal closed-form; GBM uses the log-normal closed-form; BrownianBridge uses Normal(0, σ²t(T−t)/T) returning Infinity/0 at the pinned endpoints (t = 0 or t ≥ T); PoissonProcess uses the Poisson PMF formula. All return NaN for out-of-domain inputs as documented (#879).
Process.ensemble(m, n) method: generates m independent paths of n steps each, returning an Array of m arrays each of length n+1; validates m ≥ 1 and n ≥ 1 and throws Error otherwise (#878).
covariogram(s, t) is now a required method on Process: the base class throws Error('Process.covariogram() is not implemented') when called, mirroring _next(). BrownianBridge now implements it with the closed-form formula σ²·min(s,t)·(T−max(s,t))/T for 0 ≤ s,t ≤ T, returning 0 for s > T or t > T and NaN for negative arguments. All five ran.process subclasses now provide covariogram.
covariogram(s, t) method on all four ran.process subclasses (BrownianMotion, OrnsteinUhlenbeck, GeometricBrownianMotion, PoissonProcess): returns the theoretical covariance C(s, t) = Cov(X(s), X(t)) between process values at times s and t. Returns NaN when either argument is negative. Satisfies covariogram(t, t) === variance(t) and covariogram(s, t) === covariogram(t, s) (#876).
ran.process.BrownianBridge(sigma, T, dt): Brownian bridge process conditioned to return to 0 at time T, driven by the discrete-time rule X_{n+1} = X_n − X_n·dt/(T−n·dt) + σ·√dt·N(0,1) and pinned to 0 at step N = T/dt. Exposes mean(t) (always 0) and variance(t) (σ²·t·(T−t)/T for 0 ≤ t ≤ T, 0 for t > T). reset() restores both state and internal time index; path(n) always starts from the initial condition without corrupting the caller's time index (#855).
ran.process.GeometricBrownianMotion(mu, sigma, dt): Geometric Brownian Motion with drift, using an exact discrete-time sampler (X_{n+1} = X_n · exp((μ − σ²/2)·dt + σ·√dt·N(0,1))). Starts at 1; paths stay positive by construction; log-returns are Normal((μ−σ²/2)·dt, σ²·dt) (#854).
ran.process.PoissonProcess(lambda, dt): Poisson counting process where arrivals in each interval Δt follow Poisson(λ·Δt); state is cumulative event count, guaranteed non-decreasing and integer-valued (#853).
ran.process.BrownianMotion and ran.process.OrnsteinUhlenbeck are now available as tree-shakeable subpath imports (import BrownianMotion from 'ranjs/process/brownian-motion', import OrnsteinUhlenbeck from 'ranjs/process/ornstein-uhlenbeck'), matching the per-distribution subpath export pattern.
ran.process.OrnsteinUhlenbeck(theta, mu, sigma, dt): mean-reverting stochastic process with exact discrete-time sampler (x = x·exp(−θ·dt) + μ·(1−exp(−θ·dt)) + σ·√((1−exp(−2θ·dt))/(2θ))·N(0,1)). Exposes mean(t) and variance(t) with closed-form analytical values; converges to stationary Normal(μ, σ²/(2θ)) (#846).
Demo page (demo.html) now includes an interactive Stochastic Processes section below the Distributions section. Select BrownianMotion, adjust parameters (μ, σ, dt) and path length, and see 7 semi-transparent sample paths with a theoretical mean E[X(t)] line and ±1σ envelope overlaid (#862).
ran.process.BrownianMotion(mu, sigma, dt): Brownian motion (Wiener process) with drift, using an exact O(1) discrete-time sampler (x += μ·dt + σ·√dt·N(0,1)). Exposes mean(t) and variance(t) with closed-form analytical values (#848).
Process.seed(s) method: seeds the internal PRNG for reproducible paths; delegates to this.r.seed(s) and returns this for chaining, mirroring Distribution.seed() (#861).
ran.process: new Process abstract base class (src/process/_process.js) with next(), path(n), reset(), and state() public interface; prerequisite for BrownianMotion and OrnsteinUhlenbeck (#847).
ran.test.welch(x, y, alpha) — Welch's two-sample t-test for equality of means using the Welch–Satterthwaite degrees of freedom. Returns { stat, passed } consistent with all other ran.test functions (#815).
DiscreteLaplace(p, mu) distribution: the bilateral geometric distribution supported on all integers ℤ, parameterized by decay p ∈ (0, 1) and integer location μ. PMF (1−p)/(1+p)·p^|k−μ|. Implements closed-form _pdf, _cdf, _q (inverse CDF), mean, variance, skewness (= 0), kurtosis (= (p²+4p+1)/(2p)), O(1) difference-of-geometrics sampler, and method-of-moments _fitInit (#812).
besselK(n, x) and besselKnu(nu, x) — Modified Bessel function of the second kind K_ν(x) for integer orders (combined series seeded upward recurrence + asymptotic expansion) and real orders (connection formula via besselInu for x ≤ 6, asymptotic expansion for x > 6), exported from ran.special (#809).
TruncatedExponential(lambda, a, b) distribution: exponential distribution restricted to a finite interval [a, b] (λ > 0, a ≥ 0, b > a). Subclass of Exponential. Implements closed-form _pdf, _cdf, _q (inverse CDF), mean, variance, inverse-CDF sampling, and method-of-moments _fitInit (#806).
AsymmetricLaplace(mu, sigma, kappa) distribution: a three-parameter two-sided exponential family parameterized by location μ ∈ ℝ, scale σ > 0, and asymmetry κ > 0. Reduces to Laplace(μ, σ/√2) at κ = 1. Implements closed-form _pdf, _cdf, _q (inverse CDF), mean, variance, skewness, and kurtosis, with exact inverse-CDF sampling and method-of-moments _fitInit (#805).

Changed

Code Health of src/dist/doubly-noncentral-beta.js improved from 7.89 → 10.0 by extracting the r-direction and s-direction summation loops into named private helpers (_pdfRForward, _pdfRBackward, _pdfSumOverS, _cdfRForward, _cdfRBackward, _cdfSumOverS), eliminating deep nesting, bumpy road, large method, and complex method smells.
Process.path(n) now advances the PRNG by n steps on each call (matching Distribution.sample() behaviour) instead of restoring the PRNG stream afterward. Consecutive calls return independent realizations; seeding before a call still guarantees reproducibility. Code that called path() twice without re-seeding and expected identical results will now receive two distinct paths. No deprecation cycle was applied: ran.process was introduced in the same release cycle and repeated idempotent path() calls without re-seeding have no legitimate use case (#869).
Distribution.test() now uses the Anderson-Darling test (Marsaglia & Marsaglia 2004 asymptotic series + finite-n correction, α = 0.01) instead of Kolmogorov-Smirnov for continuous distributions. The passed field is unaffected. The statistics field now carries the A² statistic (typical scale 0.2–5) rather than the KS D-statistic (scale 0–1); code that reads the raw value will silently see a different number (#816).

Fixed

CompoundPoissonProcess.seed(v) now also seeds the jump distribution's internal PRNG with a derived seed, making consecutive .path() calls with the same seed before each produce identical arrays. Previously only the arrival PRNG (this.r) was seeded; jump magnitudes drawn from jumpDist.r were not reset, silently breaking the reproducibility contract advertised by Process.seed() (#893).
gammaLowerIncomplete(s, x) no longer silently returns a truncated (wrong) value for large shape parameters (e.g. s = x = 1000). The series loop in _gli now uses an adaptive per-call iteration limit ceil(sqrt(2·(s+1)·log(1/ε))) instead of the fixed MAX_ITER = 100; for s ≈ 1000 the old cap truncated at 100 iterations while convergence requires ~265. Downstream distributions (Chi2, Gamma, Erlang, Poisson, Nakagami) that delegate CDF computation to this function are also corrected (#837).
docs/index.js now uses sass.compile() instead of the deprecated sass.renderSync(), eliminating deprecation warnings on every npm run docs invocation (#817).

[1.29.0] - 2026-07-04

Added

Documentation pages now display a "source" link next to each API entry, pointing to the exact file and line in the versioned GitHub tree (v{version}/{file}#L{line}). Links are constructed in docs/index.js from the context.file / context.loc fields exposed by documentation.js, and rendered as a small .source-link anchor in docs/templates/index.pug (#739).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for Laplace (μ, 2b², 0, 3), IrwinHall (n/2, n/12, 0, −6/(5n)), Bates ((a+b)/2, (b−a)²/(12n), 0, −6/(5n)), JohnsonSU (closed-form via E[e^{tU}] MGF with ω=exp(1/δ²)), and UniformProduct (raw moments E[X^k]=(1/(k+1))^n assembled into central moments), overriding the numerical fallback from #403 (#584).
Distribution base class now exposes mean(), variance(), skewness(), and kurtosis() methods. Each returns the theoretical value via a numerical fallback (tanh-sinh quadrature for continuous distributions, compensated summation for discrete) and can be overridden per-distribution with a closed-form formula. Cauchy overrides all four to return NaN (moments undefined) (#403).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the log-transformed distributions LogNormal, Gilbrat, LogGamma, and LogLaplace, overriding the numerical fallback from #403 and verified against mpmath (mp.dps = 50). LogNormal/Gilbrat use the standard log-normal formulas. LogGamma (Wolfram exp-gamma parameterization, X = e^Y + μ − 1 with Y ~ Gamma(α, β)) and LogLaplace (X = e^Y with Y ~ Laplace(μ, b)) derive their raw moments from the underlying gamma/Laplace MGF and return Infinity for the parameter regimes where a moment diverges (β ≤ k and kb ≥ 1 respectively for the k-th moment) (#572).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the exponential family: Exponential (1/λ, 1/λ², constants 2 and 6), Lindley (polynomial closed forms in θ), Hyperexponential (weighted raw moments E[X^r]=r!·Σwᵢ/λᵢ^r via Neumaier summation), QExponential (GP-parametrized threshold guards returning Infinity for ξ ≥ 1, ½, ⅓, ¼), and ExponentialLogarithmic (via polylogarithm Li_n(1−p)/β^(n−1)/ln(p), computed with Wynn-ε acceleration). Adds src/special/polylogarithm.js as a new special function (#574).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the Gamma/Chi family: Gamma (mean = α/β, var = α/β², skew = 2/√α, excess kurt = 6/α), InverseGamma (conditionally Infinity below shape thresholds 1/2/3/4), InverseChi2 (conditionally Infinity below ν thresholds 2/4/6/8), Chi (via Γ-function ratios, compact closed forms for skew/kurt), GeneralizedGamma and Nakagami (central moments from raw moments via logGamma), MaxwellBoltzmann (mean/var from encoded scale, skew/kurt are universal constants), and DoubleGamma (mean = 0, skewness = 0, variance = α(α+1)/β², excess kurtosis = (α+2)(α+3)/(α(α+1)) − 3). Erlang and Chi2 inherit from Gamma automatically. InverseGamma/InverseChi2 now correctly return Infinity instead of silently wrong finite values from the numerical fallback (#573).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the Beta distribution family: Beta (standard formulas), BetaPrime (with Infinity for regimes β≤1/2/3/4 respectively), BetaRectangular (mixture-of-components formula assembling central moments from Beta and Uniform components), Kumaraswamy (raw moments m_n = b·B(1+n/a, b) assembled into central moments), Arcsine (mean=(a+b)/2, var=(b−a)²/8, skewness=0, excess kurtosis=−3/2), BaldingNichols (mean=p, var=p(1−p)F, skewness/kurtosis inherited from Beta), PERT (mean=(a+4b+c)/6, var=(c−a)²αβ/252, skewness/kurtosis scale-invariant from Beta(α,β) where α+β=6 always), Wigner (mean=0, var=R²/4, skewness=0, excess kurtosis=−1), and R (mean=0, var=1/(c+1), skewness/kurtosis inherited from Beta(c/2,c/2)) (#575).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for bounded-support and uniform-shape distributions: Uniform (mean=(a+b)/2, var=(b−a)²/12, kurt=−6/5), Triangular (standard formulas, kurt=−3/5 always), Trapezoidal (raw moments from piecewise integration, assembled into central moments), UQuadratic (var=3(b−a)²/20, kurt=−38/21), RaisedCosine (var=s²(⅓−2/π²), kurt=6(90−π⁴)/(5(π²−6)²)), Anglit (var=β²(π²−8)/16, kurt=2(96−π⁴)/(π²−8)²), HyperbolicSecant (mean=0, var=1, kurt=2), Bradford (raw moments assembled via E[Xⁿ] closed form in ln(1+c)), and Reciprocal (raw moments E[Xⁿ]=(bⁿ−aⁿ)/(n·ln(b/a)), assembled into central moments) (#576).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the extreme-value and Weibull family: Weibull (λⁿ·Γ(1+n/k) raw moments assembled into central moments), DoubleWeibull (symmetric, mean=0, skewness=0, variance=λ²·Γ(1+2/k), kurtosis=Γ(1+4/k)/Γ(1+2/k)²−3), Rayleigh (special-case closed forms in σ and π, skewness and kurtosis are universal constants), ExponentiatedWeibull (generalized-binomial series E[Xⁿ]=λⁿ·α·Γ(1+n/k)·Σ(−1)ʲC(α−1,j)/(j+1)^(1+n/k), assembled into central moments), InvertedWeibull (raw moments Γ(1−n/c) with Infinity below thresholds c=1/2/3/4), Frechet (same Gamma-function thresholds on α=1/2/3/4, mean shifted by location m), Gumbel (mean=μ+βγ, variance=π²β²/6, skewness=12√6·ζ(3)/π³, excess kurtosis=12/5), and GeneralizedExtremeValue (GEV: g_r=Γ(1+r·c) with standard-GEV shape ξ=−c, Infinity below thresholds c=−1/−½/−⅓/−¼, sign(ξ) factor on skewness) (#578).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the Pareto/power-law family: Pareto and Lomax (standard tail-index threshold guards — mean for α>1, var for α>2, skew for α>3, kurt for α>4; Infinity below each threshold for var/mean, NaN for skew/kurt), GeneralizedPareto (same threshold structure against ξ: mean for ξ<1, var for ξ<½, skew for ξ<⅓, kurt for ξ<¼), BoundedPareto (all moments finite via raw-moment formula α·Lᵅ·(H^(r−α)−L^(r−α))/((r−α)·(Hᵅ−Lᵅ)) with a log limit when r≈α), PowerLaw (on [0,1], raw moments E[Xʳ]=a/(a+r) assembled into central moments), Benini (raw moments via σʳ·(r√(π/β)/2·exp(u²)·erfc(−u)+1) where u=(r−α)/(2√β)), and Champernowne (mean=x₀ by definition; skewness=0 from PDF symmetry about x₀; variance/kurtosis inherit the numerical fallback) (#577).
Analytical mean() and higher moments for six survival/reliability distributions, overriding the numerical fallback from #403: Gompertz (mean via E₁), GammaGompertz (mean via ₂F₁ inline series, falling back to quadrature for β≤0.5), Makeham (mean via Lentz CF for the upper incomplete gamma), Muth (mean=1 exact; variance via E₁), BenktanderII (mean=1+1/a; variance via regularized upper incomplete gamma), and BirnbaumSaunders (all four moments: polynomial closed forms in μ, β, γ). Adds src/special/e1.js (exponential integral E₁ via A&S 5.1.11 series for z≤1 and A&S 5.1.22 continued fraction for z>1) (#582).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the fundamental discrete distributions, overriding the numerical fallback from #403: Bernoulli (p, p(1−p), (1−2p)/√v, (1−6v)/v), Binomial (np, np(1−p), same Bernoulli skew/kurt in terms of v=np(1−p)), Rademacher (0, 1, 0, −2), Degenerate (x₀, 0, NaN, NaN — point mass has undefined higher moments), DiscreteUniform ((a+b)/2, (n²−1)/12, 0, −6(n²+1)/(5(n²−1)); NaN skewness/kurtosis for n=1), Categorical (direct O(n) sums over the normalized pdfTable), BetaBinomial (O(1) closed forms via factorial moments E[(X)_r]=(n)_r·(α)_r/(s)_r; overrides the Categorical sums), and HeadsMinusTails (E[X]=2nA, Var=2n−4n²A², μ₃/σ³ and μ₄/σ⁴ from Binomial(2n,½) cumulants where A=C(2n,n)/4^n) (#585).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the Poisson and compound-Poisson family, overriding the numerical fallback from #403: Poisson (λ, λ, 1/√λ, 1/λ), Skellam (μ₁−μ₂, μ₁+μ₂, (μ₁−μ₂)/(μ₁+μ₂)^{3/2}, 1/(μ₁+μ₂)), GeneralizedHermite (cumulant formula κᵣ=a1+mʳ·a2 from the compound-Poisson structure), NeymanA (compound-Poisson cumulants κᵣ=λ·Tᵣ(φ) via Touchard polynomials of Y~~Poisson(φ)), PolyaAeppli (compound-Poisson cumulants κᵣ=λ·Aᵣ(θ) via Eulerian polynomial moments of Y~~Geometric(1−θ) on {1,2,…}), Delaporte (κᵣ=λ+NegBin cumulants: κ₁=αβ, κ₂=αβ(1+β), κ₃=αβ(1+β)(1+2β), κ₄=αβ(1+β)(1+6β+6β²)), Borel (MGF-recurrence cumulants: var=μ/(1−μ)³, NaN skewness/kurtosis at μ=0 degenerate), BorelTanner (n-fold Borel convolution, n-scaled cumulants, NaN at μ=0), and ConwayMaxwellPoisson (raw-moment series via log-space recurrence, mode-guided stopping at λ^{1/ν}) (#586).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the geometric and discrete-series family, overriding the numerical fallback from #403: Geometric ((1−p)/p, (1−p)/p², (2−p)/√(1−p), 6+p²/(1−p); NaN at p=1), NegativeBinomial (rp/(1−p), rp/(1−p)², (1+p)/√(rp), 6/r+(1−p)²/(rp); NaN at p=0), LogSeries (closed forms assembled from raw moments E[Kⁿ]=−(1/ln(1−p))·Σk^n·p^k), FlorySchulz (sum-of-two-geometrics convolution: (2−a)/a, 2(1−a)/a², (2−a)/√(2(1−a)), 3+a²/(2(1−a))), YuleSimon (falling-factorial raw moments via Γ(ρ)/Γ(ρ−n), Stirling-2nd-kind conversion; Infinity below thresholds ρ≤1/2/3/4 for mean/var/skew/kurt), and Soliton (ideal soliton: harmonic-sum closed forms for E[Xⁿ] with n=1,2,3,4). DiscreteWeibull retains the numerical fallback — no elementary closed form exists for arbitrary β (#587).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the Logistic family: Logistic (mean=μ, var=π²s²/3, skewness=0, excess kurtosis=6/5), GeneralizedLogistic (Type I, CDF=(1/(1+e^{−z}))^c — cumulants from K(t)=log Γ(t+c)+log Γ(1−t)−log Γ(c) via digamma and Hurwitz-zeta for polygamma), HalfLogistic (raw moments E[X^k]=2·k!·(1−2^{1−k})·ζ(k) from series expansion), LogLogistic (raw moments α^n·nπ/β/sin(nπ/β) with Infinity/NaN existence thresholds at β≤1/2/3/4), ShiftedLogLogistic (B_k=Γ(1+kξ)Γ(1−kξ)=kπξ/sin(kπξ) with Infinity/NaN thresholds at |ξ|≥1/½/⅓/¼). LogisticExponential retains the numerical fallback (#579).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the F, Student-t and noncentral families, verified against mpmath (mp.dps = 50): F (existence thresholds d2 > 2/4/6/8, Infinity below — also fixes F wrongly inheriting Beta's moment overrides, e.g. mean d1/(d1+d2) instead of d2/(d2−2)), StudentT (0, ν/(ν−2), 0, 6/(ν−4) with t-style thresholds: Infinity for divergent even moments, NaN for sign-indeterminate odd ones), StudentZ (t formulas shifted to n via Z = T(n−1)/√(n−1)), NoncentralChi2 (k+λ, 2(k+2λ), √8(k+3λ)/(k+2λ)^{3/2}, 12(k+4λ)/(k+2λ)²), NoncentralChi (odd raw moments via the generalized-Laguerre/Kummer ₁F₁ form, evaluated through the Kummer transform for an all-positive series), NoncentralF (raw moments from noncentral-χ² cumulants with thresholds d2 > 2/4/6/8, Infinity below), NoncentralBeta (Poisson-weighted series of central Beta raw moments via recursiveSum), and NoncentralT (raw moments (ν/2)^{j/2}·Γ((ν−j)/2)/Γ(ν/2)·E[(Z+μ)^j] for ν > j; the ν = 3 skewness divergence carries the sign of μ, returning ±Infinity) (#581).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for InverseGaussian (mean=μ, var=μ³/λ, skew=3√(μ/λ), kurt=15μ/λ), ReciprocalInverseGaussian (via GIG(½,λ,λ/μ²) moments: mean=1/μ+1/λ, var=(1/λ)(1/μ+2/λ), closed-form skew/kurt in terms of a=1/μ, b=1/λ), Levy (all four return Infinity — every positive-order moment diverges), and Rice (mean via Laguerre L_{1/2}, variance from E[X²]=2σ²+ν²−mean², E[X³] via L_{3/2} recursion, E[X⁴] exact from noncentral-χ² connection) (#583).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for the power-law discrete family: Hypergeometric (falling-factorial moments E[(X)_m]=(n)_m·(K)_m/(N)_m assembled into central moments via Stirling-2nd-kind recurrence), NegativeHypergeometric (mixed factorial moments E[(X)_m]=r^(m)↑·K^(m)↓/(N−K+1)^(m)↑, rising r/base, falling K), Zeta (E[X^r]=ζ(s−r)/ζ(s); returns Infinity below thresholds s=2/3/4/5 for mean/var/skew/kurt), Zipf (E[X^r]=H(N,s−r)/H(N,s) via generalized harmonic numbers; h₀…h₄ precomputed in constructor), and ZipfMandelbrot (raw moments μ₁…μ₄ accumulated from the normalized pdfTable in constructor) (#588).
Analytical mean(), variance(), skewness(), and kurtosis() (excess) for Burr (raw moments k·B(k−n/c, 1+n/c) with ck>n existence thresholds), Dagum (raw moments b^r·p·B(p+r/a, 1−r/a) with a>r thresholds), Mielke (Dagum-reparametrized Beta moments (k/s)·B((k+r)/s, 1−r/s) with s>r thresholds), Davis (raw moments b^r·Γ(n−r)·ζ(n−r)/(Γ(n)·ζ(n)) with n>r+1 thresholds), FisherZ (exact cumulants from the log-MGF K(t)=½log(d2/d1)+logΓ(d1/2+t/2)+logΓ(d2/2−t/2)−… via digamma and Hurwitz-zeta polygamma), and Moyal (universal constants: mean=μ+σ(γ+ln 2), var=σ²π²/2, skewness=28√2·ζ(3)/π³, excess kurtosis=4), overriding the numerical fallback from #403 (#580).
BetaGeometric distribution: PMF f(k;α,β)=B(α+1,β+k−1)/B(α,β), support k∈{1,2,3,…}. Implements a closed-form O(1) CDF derived via a telescoping Beta-function identity (F(k)=1−B(α,β+k)/B(α,β)), replacing the previous PreComputed accumulation stub. Accessible as ran.dist.BetaGeometric(alpha, beta) (#703).
BetaNegativeBinomial distribution: PMF f(k;r,α,β)=Γ(r+k)/(Γ(k+1)Γ(r))·B(α+r,β+k)/B(α,β), support k∈{0,1,2,…}. Analytic CDF via forward recurrence, direct compound sampler (p~Beta(α,β), k|p~NegativeBinomial(r,p)). Accessible as ran.dist.BetaNegativeBinomial(r, alpha, beta) (#704).

Changed

Distribution.save() now includes type and k fields in the snapshot. Distribution.load() is now a static factory method — call ran.dist.Pareto.load(state) instead of new ran.dist.Pareto().load(state). The static form reconstructs an instance without a throw-away constructor call (#537).
fit() integer grid window for Chi2, Chi, InverseChi2, IrwinHall, UniformProduct, HeadsMinusTails, Soliton, Erlang, and F now adapts to the observed Fisher information at the seed via w = max(5, ⌈3/√I_obs⌉) instead of a fixed ±5 window. The window automatically widens for high-variance seeds (e.g., F(5, 300) from small samples) and narrows when the data strongly determines the integer parameter (#663).
Tightened distribution test tolerance to 1e-14 in test/dist.js and test/test-utils.js: refValTol now uses max(|expected|·1e-14, 1e-14) for normal-range values and a 1e-4 relative guard for sub-1e-14 reference values, the finite-difference pdf–cdf consistency check retains its own FD_FLOOR = 1e-9 floor, and reference values for Hoyt, Kolmogorov, NoncentralBeta, ConwayMaxwellPoisson, NegativeHypergeometric, and IrwinHall are updated from their pre-v1.27.0 external sources to the current double-precision computed values (#562).
Distribution._qEstimateRoot() now delegates initial-bracket selection to a new protected _qInitialGuess(p) method. For bounded-support distributions, the base implementation returns [lo, hi] (the full support range) — CDF(lo) = 0 and CDF(hi) = 1 guarantee a sign change for any 0 < p < 1, so no further expansion is needed. For unbounded/semi-bounded support a compact starting point near the support boundary is returned and _qEstimateRoot() expands it until a sign change is found (#540, #563).

Removed

src/algorithms/bracket.js (exponential bracket search). Its only caller was Distribution._qEstimateRoot(), which now handles bracket expansion inline. Removing it eliminates an API surface that was never intended for external use (#563).

Fixed

DoublyNoncentralT moment estimates (mean(), variance(), skewness(), kurtosis()) were non-deterministic across runs because the numerical integration bounds were seeded through _qInitialGuess, which consumed PRNG state. Adding _q(p) with a deterministic bracket (center at μ, spread scaled by √ν + √θ) bypasses the PRNG-seeded fallback entirely; moment test tolerances tightened from ±0.003/0.005/0.02/0.1 to ±0.001/0.002/0.01/0.05 (#800).
Distribution.load() now correctly restores PRNG state. Xoshiro128p.save() was returning a live reference to the internal _state array; any sampling after save() mutated that array and corrupted the snapshot before load() could read it. Both save() and load() now copy the state array so each saved snapshot and each loaded generator hold independent copies (#792).
BetaPrime.kurtosis() denominator factor corrected from (beta + 1) to (beta - 2), matching the Wikipedia formula. The typo caused silently wrong excess kurtosis for all valid inputs (β > 4); e.g. BetaPrime(2, 5).kurtosis() returned 66 instead of 54 (#752).
LogitNormal.mean(), .variance(), .skewness(), and .kurtosis() now override the Normal base-class fallback with explicit tanh-sinh quadrature over (0, 1) instead of inheriting Normal's hard-coded μ, σ², 0, 0 values. When μ = 0 the distribution is symmetric about x = 0.5, giving mean() = 0.5 and skewness() = 0 exactly (#756).
Rice.kurtosis() returned 0 in the Bessel overflow regime (ν/σ > ~53.3, z = ν²/(4σ²) > 709). The correct leading-order asymptotic value is −6σ²/ν² from the noncentral-χ² 4th-cumulant expansion around the Gaussian limit; e.g. ≈ −0.00167 at ν/σ = 60 (#766).
YuleSimon.skewness() now returns (ρ+1)²·√(ρ−2)/(ρ·(ρ−3)). The previous formula (ρ+1)/(ρ−3)·√((ρ−2)/ρ) was missing a factor of (ρ+1)/√ρ — e.g. at ρ=5 the old code returned 2.324 instead of 6.235 (#587).
YuleSimon.kurtosis() falling-factorial coefficients corrected: E[K^(n)] = n!·(n-1)!·ρ/∏(ρ−i), so f3 = 12ρ and f4 = 144ρ (previously 6ρ and 24ρ, missing the (n-1)! factor for n≥3). At ρ=5 the old code returned excess kurtosis ≈19 instead of 118.8 (#587).
FlorySchulz._fitInit seed corrected from 2/mean to 2/(mean+1), matching the distribution's mean formula E[X]=(2−a)/a (previously the comment and implementation used the wrong formula E[X]=2/a) (#587).
Zeta.kurtosis() now returns Infinity (not NaN) for 3 < s ≤ 4, where the variance is finite (requires only s > 3) but the 4th central moment diverges. The previous s ≤ 4 → NaN boundary conflated a divergent moment with an undefined one, contradicting the divergence convention and the sibling Zeta.skewness() logic.
Zeta.skewness() now returns NaN (instead of Infinity) when s ≤ 3 (variance is infinite, making the standardized third central moment undefined); Zeta.kurtosis() now returns NaN (instead of Infinity) when s ≤ 4 (variance is infinite or E[X³] diverges, making the fourth central moment an indeterminate ∞−∞ form). The Infinity returns for the divergent-but-determinate cases (3 < s ≤ 4 for skewness, 4 < s ≤ 5 for kurtosis) are unchanged (#769).
Alpha.mean(), .variance(), .skewness(), and .kurtosis() now correctly return Infinity/NaN (PDF f(x)~C/x² makes E[X] diverge; variance/skewness/kurtosis are undefined when the mean diverges). UniformRatio gets the same fix (PDF 1/(2x²) for x>1 gives a divergent mean). LogCauchy.mean() now returns Infinity — the integral ∫ x·f(x) dx substituting u=ln x converges only one-sidedly unlike the Cauchy case, so E[X]=+∞. DoublyNoncentralF.skewness() and .kurtosis() now guard on d2>6 / d2>8, returning Infinity when the moments do not exist — matching the analogous guards in F and NoncentralF (#570).
Slash.mean(), .variance(), .skewness(), .kurtosis() now return NaN (all moments of the Slash distribution are undefined); JohnsonSB overrides all four to use tanh-sinh quadrature over its bounded support rather than inheriting Normal's hard-coded μ=0, σ²=1, 0, 0 values. Also corrects JohnsonSB support from closed: true to closed: false (the distribution's true support is the open interval (ξ, ξ+λ)), which was the root cause of _numericalRawMoment returning NaN (tanhSinh evaluated the PDF at the exact boundary via floating-point saturation, producing 0/0) (#736).
params() now returns only natural (user-facing) parameters for Chi2, Erlang, MaxwellBoltzmann, Rayleigh, DoubleWeibull, HalfNormal, Slash, LogCauchy, and StudentZ. Previously these subclasses leaked their parent's reparametrized keys (e.g., Rayleigh(1.5).params() returned { lambda: 1, lambda2: 2.12, k: 2 } instead of { sigma: 1.5 }), violating the documented contract of params() from ADR-0014. save()/load() round-trips are preserved; all computed methods are functionally identical (#742).
logGamma now returns exact IEEE 754 results for positive integer arguments z ≤ 171 via a 171-entry LOG_FACTORIAL table (each entry independently rounded from mpmath at 50 decimal places, ≤ 0.5 ULP), eliminating the Lanczos drift that previously accumulated to 2–6 ULP when logBeta/logBinomial combined three calls. Combined with improved CDF summation strategies for BetaBinomial and NegativeHypergeometric — using the forward sum directly when CDF(x) < 0.25 (avoiding catastrophic cancellation in 1 − bwd) and Math.min(1, max(fwd, 1 − bwd)) near the midpoint — this lifts BetaBinomial and NegativeHypergeometric pmf/cdf precision from ~1e-12 to the arithmetic floor of ~2e-14 (#684).
BetaBinomial._cdf and NegativeHypergeometric._cdf midpoint path now clamps the return value to Math.min(1, Math.max(fwd, 1 − bwd)), restoring the safety clamp inadvertently removed in #684. Without this, IEEE 754 rounding in logBeta/logBinomial terms can push accumulated probability sums marginally above 1, causing survival() to return negative values (#701).
Distribution._qEstimateRoot() expansion loop no longer stalls on semi-bounded distributions (e.g. NoncentralChi2, NoncentralChi). The previous code compared |f(a)| ≤ |f(b)| to pick the expansion side, but when a was already clamped at the lower support boundary the expansion became a no-op and f(a) never changed, causing the loop to spin until exhaustion and return NaN. The fix computes newA and newB before deciding which side to expand, and only expands a side if it would actually change; the expansion logic for bracket.js is now inlined into _qEstimateRoot() and bracket.js is removed (#563).

[1.27.0] - 2026-06-05

Added

Precision regression gates for all distributions pinning pdf, cdf, and quantile accuracy against external high-precision reference values (mpmath 1.4.1, mp.dps = 50): test/precision-continuous.js covers 110 continuous distributions at 3 parameter sets × 5 interior x-values (F⁻¹(p) for p ∈ {0.1, 0.3, 0.53, 0.72, 0.9}) asserting relative error ≤ 1e-14 (1e-12 cap for series/cancellation-limited families), with quantile verified by round-trip q(cdf(x)) = x (#633); test/precision-discrete.js covers 28 discrete distributions at 3 parameter sets × up to 5 k-values with the same bounds (1e-12 cap for log-gamma/Bessel-based families) (#634). Reference values are generated by scripts/precision-refs-continuous.py and scripts/precision-refs-discrete.py.

Changed

Bernoulli, Binomial, and Rademacher now extend Distribution directly instead of Categorical, replacing alias-table construction with analytical _pdf, _cdf, and _generator implementations: Bernoulli also corrects this.k from 2 to 1, fixing aic()/bic() penalty counts (#669); Binomial computes _pmf analytically via log-space formula with explicit guards for p=0/p=1 and sample() sums n independent Bernoulli trials (#670); Rademacher eliminates unnecessary alias-table infrastructure for a two-point symmetric distribution (#668).
Eleven distributions now implement exact closed-form quantile functions (_q) instead of falling back to the non-deterministic _qEstimateRoot root-finder: Chi, DoubleGamma, GeneralizedGamma, LogGamma, MaxwellBoltzmann, GeneralizedNormal, and HalfGeneralizedNormal via gammaLowerIncompleteInv (with GeneralizedNormal/HalfGeneralizedNormal accounting for their respective CDF transforms) (#689); Levy and Moyal via erfinv, Lindley via the W₋₁ branch of Lambert W, and Reciprocal directly as a·(b/a)^p (#619). Quantile computation is now deterministic and O(1) for all eleven; Lindley, Moyal, and Reciprocal _generator() samplers also reuse _q.
Owen T iterative methods (T1–T5 in src/special/owen-t.js) now use convergence-checked loops (|term| < |sum| · ε) instead of fixed iteration counts. The ORDERS table remains as a safety cap, but series terminate early once machine precision is reached, improving accuracy near sector boundaries from ~10 to ~15 significant digits (#554).
generalizedHarmonic direct-sum path now uses Neumaier compensated summation instead of raw accumulation, and the zeta-function switch threshold is lowered from n ≥ 20 to n ≥ 10. For large m the terms drop off rapidly, causing raw accumulation to lose ~2 decimal digits; compensated summation eliminates this rounding error (#553).
RaisedCosine and Moyal _generator() now use exact inverse-CDF samplers instead of rejection sampling: RaisedCosine uses Chandrupatla's bracketed root-finder on the standardised support [−1, 1]; Moyal uses erfinv to invert the erfc-based CDF analytically. Eliminates the silent-failure risk of the 100-iteration rejection cap (#548).
lambertW1m Halley stopping criterion changed from pure-relative (|Δw/w| < ε) to hybrid absolute/relative (|Δw| < ε·max(|w|,1)) to avoid false convergence near w = 0; W₋₁ initial guess replaced with a region-adaptive approximation — Corless et al. (1996) branch-point series for z ∈ [−1/e, −0.1], logarithmic Laurent for z ∈ (−0.1, 0) — improving convergence from ~10 iterations to 2–3 (#550).
Replaced Cohen-Rodriguez-Zagier alternating-series acceleration (accelerated-sum) with Wynn's epsilon algorithm (wynn-epsilon), a general-purpose series accelerator that works on both alternating and monotone series with dynamic convergence detection and no fixed term limit. Callers in riemann-zeta.js and doubly-noncentral-t.js updated to pass signed terms (#545).
Root-finding in src/algorithms/ now uses Chandrupatla (1997) instead of the previous incomplete Brent implementation. The new algorithm correctly handles roots near zero with a mixed absolute/relative stopping criterion (|dx| < ε·max(|x|,1)), throws on invalid brackets (same-sign endpoints) instead of silently returning NaN, and returns the best approximation after iteration budget exhaustion (#541).
Replaced adaptive trapezoidal doubling (src/algorithms/trap.js) with double-exponential (tanh-sinh) quadrature (src/algorithms/tanh-sinh.js). The new algorithm achieves machine-epsilon precision in ~7 refinement levels (≤500 evaluations) versus the trapezoidal method's worst-case 2^24 evaluations, and naturally handles endpoint singularities (#542).
recursiveSum stopping criterion changed from pure relative (|Δ/sum| < ε) to hybrid absolute/relative (|Δ| < ε·max(|sum|, 1)), fixing stalled convergence when the series sum is near zero — as occurs in noncentral β, t, χ², Von Mises, and hypergeometric series at cancellation points (#543).
recursiveSum iteration cap raised from 100 (MAX_ITER) to 500 (MAX_SERIES_ITER), fixing silent truncation for series that need O(|z|) or O(λ) terms to converge — specifically the Kummer ₁F₁ Taylor series for |z| close to 50 and Poisson-weighted noncentral distribution CDFs at large noncentrality parameters (#564).
Newton's method stopping criterion changed from pure relative (|dx/x| < ε) to hybrid absolute/relative (|dx| < ε·max(|x|, 1)). The old criterion produced NaN when x = 0 (causing the loop to always run to MAX_ITER), and for |x| < 1 demanded convergence tighter than machine precision. Used by erfinv and the Marcum-Q truncation-number computation (#549).
Distribution.fit() now uses Powell's conjugate-direction optimizer instead of Nelder-Mead (src/algorithms/powell.js replaces nelder-mead.js), and returns the exact closed-form MLE directly — skipping the optimizer entirely — for 20 distributions whose estimator is closed-form: Exponential, Normal, Poisson, Bernoulli, DiscreteUniform, Pareto, LogNormal, Rayleigh, MaxwellBoltzmann, HalfNormal, Geometric, Laplace, Reciprocal, Lindley, Uniform, InverseGaussian, LogitNormal, PowerLaw, Borel, and BorelTanner. Uniform, InverseGaussian, and PowerLaw _fitInit were corrected to their exact MLEs ([min, max]; λ̂ = n/Σ(1/xᵢ − 1/x̄); â = −1/mean(log x)). See ADR-0016 (#546).
besselI(n, x) normalisation switched from the _I0(x) Taylor series (capped at 100 iterations, wrong for x ≳ 200) to an in-sweep accumulated sum S = f_0 + 2*(f_1+f_2+…) with log-exp normalisation I_n = f_n * exp(x − log(S)) (DLMF 10.35.3). Also fixes besselI(0, x) at large x, which previously returned values off by exponentially many orders of magnitude. Relative error now ≤ 1e-10 at x = 50, 100, 200 (#544).
besselInu(nu, x) now returns values accurate to near-machine precision (~1e-15 relative error) for large arguments (x = 50, 100, 200) at fractional orders ν = 0.5, 1.5, 2.3; explicit test assertions lock in that fix for the fractional-order path (#629).

Removed

Distribution.q(p) no longer returns undefined for p outside [0, 1]; it now throws Error('Invalid probability. p must be in [0, 1].'). The deprecation warning introduced in #592 (v1.26.0) is also removed (#594).
scripts/bench.js and the 11 jstat / @stdlib devDependencies that backed it. The one-time comparative benchmark (issue #114) has served its purpose; keeping the packages inflated npm install and triggered false-positive alerts on snyk scans of the repo. ADR-0011 documents the original decision and rationale.

Fixed

Normal._cdf now uses 0.5·erfc(−z/√2) instead of 0.5·(1+erf(z/√2)), eliminating catastrophic cancellation in the far left tail (≥12 digits lost at z=−7). Normal._q adds a third Newton step (was two), reducing round-trip error from ~1e-13 to machine precision at 7σ. LogNormal._q replaces erfinv(2p−1) with the same three-step Newton inversion: erfinv loses ~11 digits near p≈0, while Newton converges to machine precision even at 7σ. erfinv itself is also fixed: the Newton residual now uses a three-way split — erf(t)−x for |x|≤0.5 (no cancellation), (1−x)−erfc(t) for x>0.5, and erfc(−t)−(1+x) for x<−0.5 — keeping both operands small near the root across all |x|. Combined, precision-gate tolerances for Normal and LogNormal far-tail probes (p≈2.87e-7 and p≈1.28e-12) are tightened from 1e-3/1e-4 to 1e-14 (#808, #835).
InverseGaussian._cdf now uses erfc(−a) instead of 1 + erf(a) for the first CDF term, eliminating catastrophic cancellation in the lower tail; adds erfcx (scaled complementary error function) to ran.special to guard the second term exp(2λ/μ)·erfc(b) against overflow for large 2λ/μ; and fixes the Laplace continued-fraction iteration limit in _erfcCF and _erfcxCF from 100 to 250, ensuring convergence for all x ≥ 1. Combined, these lift InverseGaussian cdf/quantile precision from ~5e-12 to the 1e-14 gate for all parameter sets (#690).
betaIncomplete(a, b, x) backward branch (triggered when x >= (a+1)/(a+b+2) and b != 0) now returns B(a,b) − bt·CF/b instead of 1 − bt·CF/b, correctly applying the complement identity B(a,b,x) = B(a,b) − B(b,a,1−x) for the unnormalized function. Example: betaIncomplete(2, 3, 0.5) now returns 11/192 ≈ 0.0573 instead of ≈ 0.974 (#675).
Chi2, Chi, InverseChi2, IrwinHall, UniformProduct, HeadsMinusTails, Soliton, Erlang, F, and FisherZ fit() now uses a profile likelihood grid search over a ±5 integer neighbourhood of the moment seed instead of relying on Powell to cross integer step-function boundaries. For Erlang, Powell optimises the continuous rate parameter at each fixed integer shape; for F/FisherZ, all (d1, d2) pairs in an 11×11 grid are evaluated. FisherZ.fit inherits the F grid via const Cls = this subclass-safe dispatch (#624).
fit() on Beta-family distributions (Beta, BetaRectangular, BetaPrime, and other Beta subclasses) no longer converges to near-singular shape parameters via a Jeffreys-like log-barrier penalty −0.5·(log α + log β) added through a new static _fitPenalty(dist) hook on Distribution; the base-class default returns 0 (pure MLE). Five re-parametrizing Beta subclasses (F, R, PERT, BaldingNichols, FisherZ) override _fitPenalty to return 0, blocking the inherited log-barrier whose MAP bias was unintended in their native parameter spaces. See ADR-0017 (#625, #660).
BetaBinomial, Hypergeometric, and NegativeHypergeometric _cdf(x) now use bidirectional raw-PMF summation (Math.min(1, Math.max(fwd, 1 − bwd))) instead of the inherited Categorical prefix-sum table, fixing 1-ULP quantile overshoot at round probability boundaries caused by AliasTable normalisation bias (#658).
Binomial._cdf(x) now uses regularizedBetaIncomplete(n−x, x+1, 1−p) instead of the inherited Categorical prefix-sum table, fixing a 1-ULP rounding error that caused Binomial(25, 0.5).q(0.5) to return 13 instead of the correct 12 (#654).
beta(m, n) now returns exact IEEE 754 results for small positive integer arguments (min(m, n) ≤ 30) via a direct recurrence B(1,n)=1/n, B(m,n)=B(m−1,n)·(m−1)/(m+n−1), instead of routing through three logGamma Lanczos calls which accumulated a sub-ULP round-trip error. This fixes YuleSimon(3).q(0.75) returning 2 instead of the correct 1 (#653).
hurwitzZeta precision for s ∈ (1, 3] improved to within 1e-14 by two fixes: the Bernoulli correction series now uses the exact ratio recurrence T_k/T_{k-1} = (B_{2k}/B_{2k-2})·(s+2k-3)(s+2k-2)/((2k)(2k-1))/N² instead of a running product of log-Gamma factors that silently collapsed to only T₁, combined with raising n_min from 20 to 50 (#678); and the partial-sum length is now dynamic (n = max(20, min(100, ceil(1/(s−1)))) instead of fixed n = 20), eliminating 3–6 significant digit precision loss when s ∈ (1, 1.05), with an Infinity guard added for |s−1| < ε (#552).
gamma, logGamma, and digamma now return Infinity at their non-positive integer poles (previously a huge finite number, because floating-point sin(πz) / tan(πz) evaluated to ~1e-16 instead of exactly 0) and stay full-precision within 1e-6 of a pole by reducing the trigonometric argument modulo its period. logGamma(z) for z ≤ 0 is now defined as ln|Γ(z)| via the log-reflection formula instead of returning a meaningless value (#555).
riemannZeta near s=1 now uses a Stieltjes-corrected Laurent expansion to eliminate catastrophic cancellation in 1−2^(1−s): a two-term approximation (1/(s−1) + γ − γ₁·(s−1)) for |s−1| < 0.01 (#551), extended to a three-term expansion (1/(s−1) + γ₀ − γ₁·(s−1) + γ₂·(s−1)²/2, DLMF 25.2.8) for s ∈ (0.99, 1.1001) (#642). Accuracy improves from ~3e-8 relative (Wynn-epsilon) to <1e-14 for s ∈ (1.01, 1.1].

[1.26.0] - 2026-05-31

Added

Distribution.params() public method returning the natural parameters of a distribution (#516). All nine Categorical subclasses (Bernoulli, Binomial, Hypergeometric, Soliton, Zipf, ZipfMandelbrot, BetaBinomial, NegativeHypergeometric, Rademacher) now expose their own named parameters (e.g. new Bernoulli(0.7).params() → { p: 0.7 }) instead of inheriting the internal lookup state from Categorical.
Interactive distribution demo page (docs/demo.html): select from 15 curated distributions, adjust parameters, visualise histogram+PDF and empirical+theoretical CDF side-by-side (rendered with dalian), and run MLE fitting via fit() to see how closely the fitted parameters match the planted values (#504).
Distribution.fit(data) static method for maximum-likelihood parameter estimation via Nelder-Mead simplex optimizer (#404). Covers all 140 exported distributions. The static _fitInit(data) hook lets each distribution seed the optimizer from a data-aware method-of-moments estimate; the base-class fallback draws random positive values until the parameter constraints pass. Together with sample() and test(), fit() closes the full statistical cycle: define → sample → fit → test.
static _fitInit(data) data-aware seeds added to 36 continuous distributions, replacing the base-class random-retry fallback: InverseGaussian, ReciprocalInverseGaussian, Nakagami, Hoyt, Lindley, Alpha, QExponential via method-of-moments (#486); Gompertz, Makeham, Muth, BenktanderII, BirnbaumSaunders, Davis, GeneralizedExponential, Rice via best-effort data-aware seeds (#488); JohnsonSU/JohnsonSB via the Slifker-Shapiro (1980) quantile method (#440); 9 noncentral distributions seeded from central-case moment equations (#439); 10 bounded-support distributions (Bates, Triangular, Trapezoidal, PERT, BetaRectangular, UQuadratic, Uniform, UniformProduct, Anglit, RaisedCosine) seeded from sample extremes (#433).
_fitInit data-aware seeds added to 34 discrete and extreme-value distributions: 22 discrete distributions (#438), 7 Weibull/extreme-value family distributions (#434), and StudentT, StudentZ, Degenerate, Soliton, IrwinHall from single-moment inversions (#441).
gaussLegendre(f, a, b, n) algorithm: fixed-order Gauss-Legendre quadrature with precomputed nodes and weights for n=5, 10, and 20; exact for polynomials of degree ≤ 2n−1 (#402).
fit() now works on zero-parameter distributions (Gilbrat, HalfLogistic, HyperbolicSecant, Kolmogorov, Rademacher, Slash, UniformRatio): Cls.fit(data) returns a fresh instance without optimization (#427).
Categorical.fit(data) and Hyperexponential.fit(data) now return fitted instances instead of throwing; Categorical uses closed-form empirical frequencies, Hyperexponential defaults to a two-component mixture initialised by a median split (#428).
ConwayMaxwellPoisson distribution: two-parameter count distribution generalizing Poisson (ν=1), supporting both overdispersion (ν<1) and underdispersion (ν>1). Normalizing constant Z uses a log-space recurrence with running log-sum-exp to prevent silent overflow for λ ≥ ~710 (#420).
ZipfMandelbrot distribution: three-parameter finite discrete distribution with PMF (k+q)^{-s} / H_{N,s,q}, generalizing Zipf by the shift parameter q ≥ 0 (#398).

Changed

npm test now runs under nyc with coverage thresholds enforced (branches ≥ 92%, lines ≥ 98%, functions = 100%, statements ≥ 98%); the suite exits non-zero if coverage drops below these baselines (#400).

Deprecated

Distribution.q(p) called with p outside [0, 1] now emits a one-time console.warn; the current behavior (returning undefined) is unchanged this release. The method will throw in v1.27.0 (see #594).

Fixed

bracket() now returns the caller-supplied initial [a0, b0] when no root is found and either boundary is 0 (#604).
bracket, brent, and newton in src/algorithms/ now return NaN (instead of undefined) on failure paths; quickselect now throws for an out-of-range index. Distribution._qEstimateRoot propagates NaN on bracket failure (#589).
All statistics modules (location, dispersion, shape, dependence, ts) now comply with ADR-0015: mismatched-length array arguments throw Error (caller error), indeterminate results (empty sample, zero-variance) return NaN, divergent results (KL divergence with Q=0,P>0, zero-denominator odds ratio) return Infinity. undefined is no longer returned as a failure sentinel from any public function (#593).
Davis._cdf(x) now uses a dual Bose-Einstein series (upper incomplete gamma series for x near μ, Bernoulli/Laurent series for large x) instead of Romberg integration; per-call cost drops from ~100ms to ~10µs, enabling full goodness-of-fit coverage with the standard sample size (#451).
Bates.fit() now uses a profile likelihood grid search over integer n instead of the inherited 3-parameter Nelder-Mead, which stalled on the staircase likelihood surface caused by integer rounding; n is now reliably recovered (#481).
Bradford._fitInit: small-c mean approximation coefficient corrected from 3·(1−2·mean) to 6·(1−2·mean), matching the correct first-order expansion E[X] ≈ ½ − c/12; the previous coefficient underestimated the starting value by a factor of 2 (#498).
NoncentralChi: .p.lambda now stores λ instead of λ²; both construction and .fit(data) now correctly return the user-facing noncentrality parameter (#491).
Mielke: aic()/bic() now use the correct parameter count of 2 (not 3 inherited from Dagum), and .p now correctly exposes { k, s } instead of Dagum's { p, a, b } bag (#480, #505).
Multi-level Distribution subclasses now report the correct free-parameter count k, fixing aic()/bic() for Weibull (2), ExponentiatedWeibull (3), Chi2 (1), MaxwellBoltzmann (1), GeneralizedGamma (3), LogGamma (3), Rayleigh (1), and HalfGeneralizedNormal (2) (#510).
R, F, FisherZ, and BaldingNichols .fit() no longer inherit Beta's wrong-arity initializer; each now seeds Nelder-Mead from a distribution-correct method-of-moments estimate (#441).
besselISpherical(n, x) now uses a Taylor series for |x| < 1, eliminating catastrophic cancellation in the closed-form expressions for n ≥ 1; relative error bounded by 2ε (#425).
Davis distribution: sample() now produces genuine random variates via an exact Zeta-Gamma mixture sampler instead of always returning 1; parameter constraint tightened to n > 1; _pdf NaN guard added for the lower support boundary (#447, #448).
neumaier() now returns ±Infinity (instead of NaN) when the input array contains ±Infinity; fixes lnL(), aic(), and bic() silently returning NaN when any observation falls outside a distribution's support (#442).

[1.25.0] - 2026-05-25

Added

bench/ directory with bench/index.js: a performance comparison script benchmarking ranjs against jStat and @stdlib/stats/base/dists across Normal, Gamma, Beta, Poisson, and Exponential distributions for sample, pdf, cdf, and quantile operations. Run with npm run bench. Closes #114.
TypeScript declarations are now generated from JSDoc via tsc --allowJs --declaration --emitDeclarationOnly as part of npm run build. The generated dist/index.d.ts replaces the hand-written dist/ranjs.d.ts, making type drift structurally impossible. Includes @overload annotations for sample(), float(), int(), choice(), shuffle(), and coin(). Closes #170.
DoublyNoncentralChi2 distribution: the law of X = U + V with U ~ ncχ²(k1, λ1) and V ~ ncχ²(k2, λ2) independent. Because the non-central chi-square is closed under addition, DoublyNoncentralChi2(k1, k2, λ1, λ2) is exactly ncχ²(k1 + k2, λ1 + λ2); it is implemented in that collapsed closed form rather than via a double Poisson series. Closes #228.
Distribution._qEstimateWalk(p, start) protected helper: deterministic linear walk from a caller-supplied integer start toward the infimum discrete quantile. Exits when cdf(k) >= p and cdf(k-1) < p. Provides a non-random alternative to _qEstimateRoot for infinite-support discrete distributions with analytically-known parameters. Closes #284.
Property tests for all distributions: cdfMonotonicity now asserts cdf(x₂) >= cdf(x₁) across a deterministic grid (it was previously a no-op that only asserted scalar arithmetic ordering). A new Tests.quantileRoundtrip helper asserts |cdf(q(p)) − p| < 1e-6 for continuous distributions and the two-sided infimum definition for discrete distributions across the fixed probability grid {0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95}. Closes #212.
R distribution generator, PDF, and CDF now produce a symmetric distribution matching the documented formula f(x; c) = (1 - x²)^(c/2-1) / B(1/2, c/2) on [-1, 1]. The previous implementation configured Beta(0.5, c/2) (the squared-variable parent) but then applied the affine substitution y = (x+1)/2 and squared it, breaking x → −x symmetry. For c=4, pdf(-0.95) returned ~0.7495 instead of the correct ~0.0731. Reduced to the affine U = (X+1)/2 ~ Beta(c/2, c/2), which is one-to-one and avoids the 0·∞ corner at x=0 for c<2. refVals for R(4) (previously deferred) added to the test suite. Closes #261.
Soliton distribution support truncation fixed: the weight array was built with length: N-2, silently omitting k=N and causing the Categorical base class to renormalize the remaining weights upward. Changed to length: N-1 so pmf(1) returns the correct 1/N and pmf(N) returns 1/(N(N-1)). Closes #263.
Catastrophic cancellation in _cdf near the lower support boundary fixed for 20 distributions: FlorySchulz (naive 1 − (1−a)^k·(1+ka) rewritten with expm1/log1p, #248); Moyal (Q(½,z) now routed through gammaUpperIncomplete directly, #247); Rice and NoncentralChi2 (complementary Marcum Q computed via new marcumP export instead of 1 - marcumQ, #246, #245); InverseGamma and InverseChi2 (upper tail via gammaUpperIncomplete directly, #244, #243); and 13 distributions using Math.expm1/Math.log1p/Math.tanh builtins: Exponential, Benini, Gompertz, Hyperexponential, GeneralizedPareto, Lomax, Burr, Pareto, GammaGompertz, Makeham, GeneralizedExponential, HalfLogistic, LogisticExponential, Muth (#214). The checkRefVals test helper now falls back to relative tolerance for sub-precision expected values. Boundary-region reference values added for all fixed distributions.
scipy/numpy/R reference values (refVals) added to 55 distributions across all families, computed from scipy 1.17.1 / R and never transcribed from the ranjs source: 10 chi/t/F/gamma-family (#126), 9 normal and logistic variants (#127), 9 extreme-value and heavy-tail (#128), 12 miscellaneous continuous including InverseGaussian, LogGamma, Skellam, FisherZ (#130), 7 noncentral (#133), and 8 core discrete (#134). Three non-trivial scipy parameterization differences (Geometric indexing, NegativeBinomial success/failure swap, Hypergeometric name collision) are documented in solutions/testing/2026-05-18-1443-discrete-refvals-scipy-parameterization-traps.md.
eslint-plugin-jsdoc added as a devDependency with a jsdoclint npm script and a new parallel CI job. Enforces JSDoc presence on public Distribution methods and exported namespace functions; catches stale @param/@returns after signature changes. Fixes @return → @returns in 12 files and adds JSDoc to two previously undocumented internal helpers. Closes #202.
Docs site now includes a scipy.stats → ranjs porting guide (porting-scipy.html) with side-by-side Python/JavaScript examples for the 20 most-used scipy distributions, a method-mapping table, and callouts for non-trivial parameter differences (LogNormal, Weibull, Triangular, Geometric, Hypergeometric).
Per-distribution subpath exports: import Normal from 'ranjs/dist/normal' now resolves correctly in Node.js ESM, browsers, and bundlers (Vite, Webpack, esbuild). Each of the 134 exported distributions has a corresponding self-contained ESM bundle at dist/<name>.esm.js (≈8–15× smaller than importing the full library). See ADR-0005.
Automated npm publish workflow (.github/workflows/release.yml): pushing a v* tag now runs lint, typecheck, and tests before publishing to npm with provenance attestation. Requires an NPM_TOKEN secret in repository settings.
TypeScript type declarations added (dist/ranjs.d.ts). All 135 distribution classes, the Distribution base class (17 public methods), and the core, location, dispersion, shape, dependence, and test namespaces are now fully typed. "types": "./dist/ranjs.d.ts" added to package.json at the top level and inside "exports" for full compatibility with all TypeScript moduleResolution modes. See ADR-0003.
Build now produces three artifacts: dist/ranjs.esm.js (ES module), dist/ranjs.cjs.js (CommonJS), and dist/ranjs.min.js (UMD, minified, CDN). "exports" field added to package.json routing import to ESM and require to CJS. "sideEffects": false added to enable tree-shaking across the 130+ distribution classes.

Changed

Distribution JSDoc references de-duplicated: removed 135 @see <url> tags whose URL was identical to the inline {@link <url>} already present in the distribution's description. The rendered "References" section now appears only when it carries non-trivial information (papers, books, algorithm citations) rather than repeating the distribution-name link. Four // Source: code comments (Alpha, Anglit, Arcsine, BaldingNichols) were promoted to class-level @see entries so the cited works show up in the rendered docs.
Gamma.q(p) (and Chi2.q(p), InverseGamma.q(p)) Wilson-Hilferty seed now uses an A&S §26.2.17 rational approximation instead of erfinv, eliminating nested Newton iteration inside the outer Halley loop and accelerating quantile throughput. Closes #384.
Normal.q(p) now uses an A&S §26.2.17 rational approximation seed with two fixed Newton refinement steps instead of erfinv, removing the convergent Newton loop and reducing from 3–5 erf evaluations to exactly 2. Closes #385.
StudentT.q(p) now uses a 4-term Cornish-Fisher seed (A&S §26.7.8) + Halley refinement instead of the previous 2-term seed + Newton, reducing CDF evaluations from ~5 to ~2 and yielding ~2.4× additional throughput improvement. Closes #381.
StudentT.q(p) now uses a Cornish-Fisher seed + Newton refinement instead of Brent root-finding, reducing CDF evaluations from ~20 to ~3 and yielding ~7× throughput improvement. Closes #368.
Normal variate generator now uses the Improved Ziggurat algorithm (Marsaglia-Tsang 2000 / Doornik 2005) instead of Box-Muller, achieving 2–4× throughput improvement on Normal.sample() and all distributions that derive normal variates (Levy, NoncentralT, DoublyNoncentralT, InverseGaussian, Gamma-family). Public API unchanged. Also eliminates a latent NaN when the PRNG returned exactly 0 (log(0) in Box-Muller). Closes #369.
Distribution internal field this.t renamed to this._type for readability. No behavioural change. Closes #205 (PR 1/6).
marcumQ and marcumP now evaluate the transition band y ≈ x + μ with μ ≥ 135 via the large-μ uniform asymptotic expansion (Section 4.2 of Gil, Segura & Temme, arXiv:1311.0681) instead of the O(μ) three-term recurrence. This is the fifth and final computation branch of the source algorithm and removes the recurrence's accumulated rounding in that regime. Closes #315.
DoublyNoncentralChi2 now extends NoncentralChi2 instead of reimplementing its PDF and CDF. The _pdf, _cdf, and _generator are fully inherited. The model complexity used by aic() and bic() changes from 4 to 2, reflecting that the distribution has 2 identifiable parameters in its collapsed form NoncentralChi2(k1+k2, λ1+λ2). NoncentralChi2 now also accepts lambda = 0 (was lambda > 0), degenerating correctly to a central chi-squared. Closes #316.
BREAKING: All 131 distribution constructors now require their parameters — default values have been removed. Constructing a distribution with no arguments (e.g. new Normal(), new Exponential()) now throws Error('Invalid parameters. Required parameters missing or not a number: ...') instead of silently using arbitrary defaults. The 7 distributions that genuinely have no parameters (Gilbrat, HalfLogistic, HyperbolicSecant, Kolmogorov, Rademacher, Slash, UniformRatio) still construct with no arguments. Distribution.validate() now rejects undefined and NaN parameter values as a centralized fail-fast guard. See ADR-0004 and #50.
docs/index.html is no longer tracked in the repository; it is now built and deployed to GitHub Pages automatically via actions/deploy-pages@v4 on every push to main. The docs-build CI job now uploads via actions/upload-pages-artifact@v3 instead of actions/upload-artifact@v4.
CI now runs npm run build on every push to main and every pull request; a build badge scoped to the build job was added to README.md.
Replaced hand-rolled SVG pixel math in .github/scripts/gen-badge.js with badge-maker; removed legacy .circleci/config.yml.
Upgraded rollup from ^2.64.0 to ^4.x. Replaced unmaintained rollup-plugin-terser with @rollup/plugin-terser. Upgraded @rollup/plugin-node-resolve from ^13.x to ^16.x.
Docs build (npm run docs) is now driven by a pages array in docs/index.js; adding a page is one array entry plus one Pug template that extends the new shared layout docs/templates/_layout.pug. The compiled SCSS is written once to docs/styles/style.css and linked externally from every page (previously inlined into each rendered HTML). See ADR-0002.
Removed dead coveralls devDependency and its coveralls npm script (was never wired into CI).
Removed npm from devDependencies (unconventional; runner's npm is used directly).
Upgraded nodemon from ^2.0.15 to ^3.0.0 to fix a semver ReDoS vulnerability in simple-update-notifier.
Fixed 24 of 41 npm audit vulnerabilities via npm audit fix.

Deprecated

ran.dist.Hoyt is deprecated. It was implementing the Nakagami-m distribution under the wrong name; ran.dist.Nakagami is the canonical, correctly-named class. new Hoyt(q, omega) now emits a console.warn and delegates entirely to Nakagami(q, omega). The parameter constraint changes from 0 < q ≤ 1 to q ≥ 0.5 (the Nakagami-m domain); computed values are identical for all previously valid q ∈ [0.5, 1]. Hoyt will be removed in a future major release. Closes #226.

Removed

Hand-written dist/ranjs.d.ts removed from version control (now a build artifact).
scripts/check-declarations.js deleted (structural completeness now guaranteed by tsc).
Distribution base class now exposes bounded(), returning 'bounded', 'lower', 'upper', or 'unbounded' based on whether the support endpoints are finite. type() and support() are documented as stable public API. TypeScript declarations updated accordingly. Closes #119.
GeneralizedPareto, ShiftedLogLogistic, and TukeyLambda GoF sampling tests now cover the boundary branches (xi=0 / lambda=0) in _q, exercising the −log(1−p), logistic, and log(p/(1−p)) code paths respectively. Closes #270.

Fixed

Docs build now locates the ran module entry in documentation's output by kind/name instead of root[0], restoring the API documentation section, sidebar menu, and search list on docs/index.html.
Gamma.q(p) now uses a dedicated Wilson-Hilferty initial estimate + Halley refinement algorithm instead of the generic Brent root-finder, eliminating the 2.5–3× quantile overhead. Chi2.q(p) and InverseGamma.q(p) benefit automatically. gammaLowerIncompleteInv is now exported from ran.special. Closes #367.
Quantile throughput restored for 10 derived distributions (BirnbaumSaunders, DoubleWeibull, ExponentiatedWeibull, JohnsonSB, JohnsonSU, LogCauchy, LogLaplace, LogNormal, LogitNormal, TruncatedNormal): super._q calls replaced with inlined closed-form formulas, eliminating the V8 megamorphic deoptimization that caused up to 56× slowdown. Closes #366.
HeadsMinusTails now rejects n = 0: constraint tightened from n >= 0 to n > 0, matching the documented domain $n \in \mathbb{N}^+$. Closes #363.
InverseGamma: removed unused this.c.betaAlpha pre-computation (Math.pow(beta, alpha)) that was computed on every construction but never read. Closes #373.
docs/porting-scipy.html styling now matches the API page: h2 section headings, standalone table layout, side-by-side .code-pair code blocks (two columns ≥ 800 px, stacked on mobile), .callout warning blocks, and sidebar parity for the static (non-checkbox) jump menu. Closes #209.
NegativeBinomial _pdf(0) returned NaN at p=0 (0 * -Infinity), and _generator() returned undefined at p=1 (Poisson(Infinity)). Added degenerate-case guards in _pdf, _cdf, and _generator for p=0 (all mass at k=0), and tightened the parameter constraint from p ≤ 1 to p < 1 (p=1 yields an all-zero PMF and no valid distribution). Closes #145.
Champernowne distribution was a non-functional stub: _generator() returned undefined, _cdf(x) always returned 1, and _pdf(x) lacked its normalization constant. Fixed all three: normalization constant is now alpha * sqrt(1 - lambda²) / (2 * arccos(lambda)), CDF uses the closed-form arctan(k * tanh(...)) formula, and _generator() uses inverse-transform sampling via a new closed-form _q(p). The class is now exported from src/dist/index.js and declared in dist/ranjs.d.ts. Closes #337.
BenktanderII near-boundary refVals at x = 1+1e-6 and x = 1+1e-4 (params [2, 0.9995]) were replaced with values derived independently via Python Decimal at 60 decimal places using the direct mathematical formula, not the expm1-based implementation formula. Closes #295.
romberg returned the silent sentinel 0 when the 20-step budget was exhausted without convergence — indistinguishable from a genuine zero integral. It now returns the best Richardson extrapolate accumulated so far, consistent with how trap returns its last estimate on timeout. The stray console.log in Davis._cdf (which exposed this bug during development) has also been removed. Closes #312.
marcumQ and marcumP were accurate only for x < 30 — the asymptotic, recurrence and quadrature computation branches existed only as commented-out scaffolding. Activated all four methods of Gil, Segura & Temme (arXiv:1311.0681) behind a regime dispatcher: series expansion (§3), large-ξ asymptotic expansion (§4.1), three-term recurrence relation (Eq. 14) and trapezoidal quadrature (§5). The Marcum functions are now accurate across the full μ ≥ 1, x > 0, y > 0 domain, restoring full-range CDF precision for the Rice, NoncentralChi2, DoublyNoncentralChi2 and Skellam distributions. Closes #253.
neumaier sorted its input array in place, silently reordering the caller's array as a side effect. It now sorts a shallow copy, leaving the original array untouched. Closes #313.
DoublyNoncentralBeta._pdf and ._cdf returned NaN when lambda1 or lambda2 was 0, because the outward-summation Poisson-weight initialisation evaluated 0 * Math.log(0) = NaN (IEEE 754). Added early-return guards: when lambda1 = 0 the double sum collapses to NoncentralBeta(beta, alpha, lambda2) at (1-x); when lambda2 = 0 it collapses to NoncentralBeta(alpha, beta, lambda1) at x. DoublyNoncentralF (which inherits both methods) is fixed implicitly. Closes #304.
NoncentralBeta._pdf and ._cdf returned NaN for lambda = 0 because the Poisson weight computation evaluated 0 * Math.log(0) = NaN (IEEE 754). Added a guard: when lambda / 2 === 0, the weight for the sole k=0 term is 1. Closes #267.
besselI(n=1, x) had only ~8 significant digits of accuracy due to a polynomial approximation (_I1, from Numerical Recipes) with limited-precision coefficients. Replaced with Miller's backward recurrence — the same algorithm used for n≥2 — and extended the loop upper-bound by ⌈2|x|⌉ to ensure the recurrence contracts the K_n component before reaching n=1 (required when |x| > n). Added odd-function sign correction for negative arguments. Fixes ~3.7e-10 CDF error in VonMises at intermediate x; expands VonMises refVals from 3 to 11 reference points. Closes #255.
BenktanderII._cdf lost precision near the lower support boundary (x ≈ 1) due to catastrophic cancellation in 1 − exp(arg) and 1 − xᵇ⁻¹·exp(u) when their arguments approach zero. Rewrote using Math.expm1 for the b=1 branch and a split (1−xᵇ⁻¹) − xᵇ⁻¹·expm1(u) decomposition for the general branch, eliminating the cancellation. Closes #242.
Bernoulli._q returned 0 for all p > 0.5 because this.p.p is undefined after the Categorical parent constructor overwrites this.p with { n, weights, min }. Fixed by using this.p.weights[0] (the CDF at k=0) as the threshold. Closes #212.
DiscreteUniform._q, Geometric._q, and DiscreteWeibull._q returned a quantile one too large when p landed exactly on a CDF step (e.g., Geometric(0.5).q(0.5) returned 1 instead of 0). Each used Math.floor on the algebraic inverse k+1; changed to Math.ceil(…) - 1 which is identical for non-integer arguments but correct at exact integers. Closes #212.
Skellam._q applied Math.floor to the result of _qEstimateRoot, which finds a continuous root of CDF(x) − p. For a step function the root lands just below the integer boundary, causing Math.floor to undershoot by 1. Added a one-step correction: if (this.cdf(k) < p) k++. Closes #212.
Skellam._q could silently return NaN in the extreme tails: _qEstimateRoot uses a random bracket initialisation and returns undefined when its 100-iteration cap is exhausted, and Math.floor(undefined) is NaN. Replaced with _qEstimateWalk(p, Math.floor(μ₁ − μ₂)), anchored at the distribution mean; the walk is fully deterministic and always returns a valid integer. Closes #283.
NoncentralBeta._pdf(0) and ._cdf(0) now return 0 instead of NaN at the closed-support lower boundary. The series in recursiveSum hits a 0/0 indeterminate form at exact x = 0, but the mathematical limit is 0 for alpha > 1. Added boundary guard and { x: 0, pdf: 0, cdf: 0 } to NoncentralBeta refVals. Closes #230.
DoublyNoncentralT._pdf(0) returned NaN when mu !== 0. Added an x === 0 guard analogous to the one already present in NoncentralT._pdf; the j=0-only closed form exp(c[0]) · Γ((ν+1)/2) · ₁F₁((ν+1)/2, ν/2; θ/2) is now returned directly. Closes #229.
FisherZ constructor now passes (d1, d2) to the F base class instead of (d1/2, d2/2). Previously, new FisherZ(d1, d2) silently produced Fisher's z with degrees of freedom (round(d1/2), round(d2/2)) (e.g. FisherZ(5, 5) was internally F(3, 3)). Internal-consistency tests passed because all distributional self-checks used the same wrong d.o.f. The bug was surfaced by adding scipy-independent reference values per #130. Closes #130.
MaxwellBoltzmann constructor now passes rate = 1/(2a²) to the Gamma base class (was incorrectly 2a²). The distribution was previously self-consistent but produced values from the wrong density; the correct PDF is f(x;a) = sqrt(2/π) x² exp(−x²/2a²) / a³. Closes #126.
erf and erfc in src/special/error.js now use a hybrid Taylor series (|x| ≤ 2) / Laplace continued fraction (|x| > 2) instead of delegating to gammaLowerIncomplete/gammaUpperIncomplete. This fixes relative precision loss in the tails (5σ+) and resolves the // TODO Replace with continued fraction comments. Adds far-tail Normal(0, 2) reference values at x = ±10 and ±14 to the test suite. Closes #211.
npm test now works on Node 20+ by replacing the unmaintained esm loader with @babel/register, aligning the test and coverage execution paths.
Kolmogorov.pdf(0) and Kolmogorov.cdf(0) now correctly return 0. Previously the lower support bound was declared closed: true (contradicting the documented support x > 0), causing cdf(0) to evaluate a non-convergent Grandi's series and return −1.
FisherZ.pdf(x) no longer returns Infinity for right-tail values at default parameters (d1=1, d2=1). Replaced delegating computation through F→Beta with a direct log-space formula that avoids float64 precision loss when the Beta argument rounds to 1.0.
NegativeBinomial constructor now correctly rejects out-of-range parameters: r ≤ 0, p < 0, and p > 1. Previously some values slipped through validation.
Gamma sampler now runs Marsaglia-Tsang directly at shape α = 1 instead of routing through the Gamma(α+1) · U^(1/α) boost branch. The boost is mathematically exact but consumes an extra PRNG draw per sample, which pushed the seed-42 KS statistic just over the p=0.01 critical value at N=10000. Fix transitively repairs sampling-test failures for Gamma, Chi, Chi2, Erlang, InverseGamma, LogGamma, Nakagami, and GeneralizedGamma at their default parameters (#193).
SkewNormal sampler now draws both Box-Muller outputs from a single uniform pair instead of calling _normal twice (which discarded one branch per call). Halves PRNG consumption per sample and resolves the seed-12345 KS failure for the positive-shape-parameter case (#195).
NoncentralF.pdf(0) and NoncentralF.cdf(0) now correctly return 0. Previously the delegating computation through NoncentralBeta produced NaN at the closed-support boundary x = 0 (NoncentralBeta.{pdf,cdf}(0) returns NaN; tracked separately as #230). Added closed-form guard in _pdf/_cdf. This also eliminates a flaky quantile-test failure where _qEstimateRoot's bracket-search probed cdf(0) during expansion and propagated the NaN through Brent's method (#233).

Accepted risks (pending follow-up issues)

mathjax-node-page ≥1.4.1 and its transitive chain (form-data, mathjax, qs, tough-cookie, yargs-parser) retain known vulnerabilities. The fix requires downgrading to mathjax-node-page@2.0.0 (breaking change). These are docs-only tools with no impact on library users; issue #116 will replace this toolchain.
serialize-javascript ≤7.0.4 remains in mocha and rollup-plugin-terser. The mocha fix requires a major upgrade (tracked in #99); the rollup-plugin-terser fix is a downgrade and tracked under #107.
vue-template-compiler ≥2.0.0 remains in documentation. The fix would downgrade to documentation@6.2.0; issue #116 will replace this tool.

FilesExpand file tree

CHANGELOG.md

Latest commit

History

CHANGELOG.md

File metadata and controls

Changelog

[Unreleased]

Added

Changed

Fixed

[1.32.0] - 2026-08-09

Added

Changed

Deprecated

Fixed

[1.31.0] - 2026-07-20

Added

Changed

Deprecated

Removed

Fixed

Security

[1.30.0] - 2026-07-09

Added

Changed

Fixed

[1.29.0] - 2026-07-04

Added

Changed

Removed

Fixed

[1.27.0] - 2026-06-05

Added

Changed

Removed

Fixed

[1.26.0] - 2026-05-31

Added

Changed

Deprecated

Fixed

[1.25.0] - 2026-05-25

Added

Changed

Deprecated

Removed

Fixed

Accepted risks (pending follow-up issues)