Release v1.32.0 · synesenom/ran

This release adds three new distributions (Tweedie, ExponentiallyModifiedGaussian, WrappedCauchy), three new goodness-of-fit tests (cramerVonMises, kolmogorovSmirnov, andersonDarling), a major expansion of the ran.process.Process API (fit(), lnL(), marginal(), params(), plus a new stochastic-process precision gate), and ran.dist.guess() for automatic candidate-distribution selection by BIC weight. ran.dist.VonMises gains a location parameter (mu) — the old single-argument constructor is deprecated and will be removed in v1.33.0. The remainder is a large batch of precision, overflow/underflow, and convergence fixes across Bessel/Marcum-Q/incomplete-gamma special functions and several noncentral distributions (NoncentralT, DoublyNoncentralT, DoublyNoncentralBeta, Skellam, and others), plus hot-path performance caching across ~35 distributions.

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.

Navigation Menu

v1.32.0

Choose a tag to compare

Sorry, something went wrong.