rán · 然 — "so; correct"
A comprehensive JavaScript library for probability distributions, random variate generation, and statistical analysis.
- Features
- Installation
- Usage
- API Overview
- Distribution API
- Process API
- MC API
- Return values and errors
- Numerical precision
- Bundle size budget
- Documentation
- Support policy
- License
- 146 probability distributions — continuous and discrete, each with PDF/PMF, CDF, quantile (
q), hazard, survival, log-likelihood (lnL), AIC/BIC, goodness-of-fit testing, and MLE fitting (fit) - Statistical measures — location (mean, median, mode, …), dispersion (variance, IQR, Gini, …), shape (skewness, kurtosis, …), and dependence (Pearson, Spearman, Kendall, …)
- Hypothesis tests — Bartlett, Levene, Brown–Forsythe, Cramér-von Mises, Mann–Whitney U, HSIC
- Reproducible sampling — every distribution accepts an optional seed for deterministic output
- TypeScript support — declaration files generated from JSDoc, covering all public APIs
- Tree-shakeable — import individual distributions without pulling in the full bundle
npm install ranjs<script src="https://cdn.jsdelivr.net/npm/ranjs@1.31.0/dist/ranjs.min.js"></script>Pin an exact version — ranjs follows a numpy/scipy-style versioning policy where breaking changes can ship in
minor releases, so an unpinned or range-based CDN URL can change behavior under you. Bump 1.31.0 to whichever
version you're targeting; the full API docs let you browse any past release.
The library is exported globally as ranjs.
Import a single distribution for minimal bundle size:
import Normal from 'ranjs/dist/normal'
const n = new Normal(0, 1)
n.pdf(0) // => 0.3989422804014327
n.cdf(1.96) // => 0.9750021048517796
n.sample(5) // => [0.42, -1.03, 0.17, 1.81, -0.55]const ran = require('ranjs')
const skellam = new ran.dist.Skellam(1, 3)
const values = skellam.sample(1e4)
skellam.test(values)
// => { statistics: 14.025360669436635, passed: true }
for (let k = -4; k <= 4; k++) {
console.log(k, skellam.pdf(k), skellam.cdf(k))
}
// => -4 0.10963424740027695 0.21542206959904264
// -3 0.16622843570192460 0.38165050508716936
// -2 0.20277318483535026 0.58442368966117290
// ...Every distribution can be individually seeded:
import Gamma from 'ranjs/dist/gamma'
const g = new Gamma(2, 1)
g.seed(42)
g.sample(3) // always produces the same sequenceconst ran = require('ranjs')
const data = new ran.dist.Skellam(1, 3).sample(1e4)
const fitted = new ran.dist.Skellam(1, 3)
const misfit = new ran.dist.Skellam(1.2, 7.5)
console.log(fitted.aic(data)) // => 41937.67252974663
console.log(misfit.aic(data)) // => 66508.74299363888ranjs closes the full statistical cycle — define a model, generate data, fit parameters from data via MLE, then verify the fit:
import { dist } from 'ranjs'
// 1. Define and sample
const model = new dist.Normal(3, 1).seed(42)
const data = model.sample(500)
// 2. Fit parameters from data via MLE
const fitted = dist.Normal.fit(data)
console.log(fitted.p) // => { mu: 3.000, sigma: 1.000 }
// 3. Test goodness of fit
console.log(fitted.test(data)) // => { statistics: 0.42, passed: true }fit() is a static method called on the class, not on an instance: dist.Normal.fit(data), not model.fit(data). All 146 exported distributions support fit(). Most have a data-aware initial guess for reliable MLE convergence; zero-parameter distributions skip optimization and return a fresh instance.
Full walkthrough (chaining aic/bic/sample off a fitted instance, a discrete-distribution example): Parameter estimation guide.
| Namespace | Contents |
|---|---|
ran.dist |
146 probability distributions |
ran.process |
Stochastic processes: AR(1), Brownian motion, Brownian bridge, Cox–Ingersoll–Ross, geometric Brownian motion, Ornstein–Uhlenbeck, Poisson process, random walk |
ran.mc |
MCMC samplers (random-walk & adaptive Metropolis, slice, HMC, NUTS, MALA, Gibbs, adaptive rejection sampling), parallel tempering, multi-chain runner, Gelman–Rubin convergence diagnostic |
ran.location |
Mean, median, mode, geometric mean, harmonic mean, trimean, midrange |
ran.dispersion |
Variance, standard deviation, IQR, Gini coefficient, entropy, CV, … |
ran.shape |
Skewness, kurtosis, quantiles, moments, min, max, rank |
ran.dependence |
Pearson, Spearman, Kendall, distance correlation, Kullback–Leibler, … |
ran.test |
Bartlett, Levene, Brown–Forsythe, Cramér-von Mises, Mann–Whitney U, HSIC |
ran.core |
Seeded PRNG (xoshiro128+), uniform float/int/bool generators |
Every distribution exposes a consistent interface:
const d = new ran.dist.Gamma(2, 1)
d.type() // 'continuous' or 'discrete'
d.params() // current parameter object, e.g. { alpha: 2, beta: 1 }
d.support() // [{ value, closed }, { value, closed }] — lower/upper bounds
d.sample(n) // generate n random variates
d.pdf(x) // probability density / mass function
d.cdf(x) // cumulative distribution function
d.q(p) // inverse CDF (quantile function)
d.survival(x) // complementary CDF (1 − CDF)
d.hazard(x) // hazard rate (pdf / survival)
d.cHazard(x) // cumulative hazard (−log survival)
d.lnPdf(x) // log probability density / mass
d.lnL(data) // log-likelihood over an array of observations
d.aic(data) // Akaike information criterion
d.bic(data) // Bayesian information criterion
d.test(data) // Anderson-Darling test (continuous) or chi-squared test (discrete)
d.seed(value) // set PRNG seed; returns the instance
ran.dist.Gamma.fit(data) // static — MLE fit; returns a new instancesave() and load(state) let you snapshot and restore the exact PRNG state and parameters of a distribution instance, so a sequence of samples can be reproduced exactly across sessions or process restarts.
const d = new ran.dist.Gamma(2, 1)
d.seed(42)
d.sample(10) // advance the internal PRNG
const state = d.save() // plain object: { type, params, prngState, ... }
const d2 = ran.dist.Gamma.load(state) // new instance with identical state
d.sample(5) // some sequence of variates
d2.sample(5) // identical sequence — same PRNG position, same parametersEvery process in ran.process extends a common Process base class and exposes the same interface:
const bm = new ran.process.BrownianMotion(0, 1, 0.1) // mu=0, sigma=1, dt=0.1
bm.next() // advance one step; returns the new state
bm.reset() // reset to initial state
bm.state() // current state value
bm.path(100) // generate a path of 100 steps; returns array of 101 states
bm.ensemble(5, 100) // generate 5 independent paths of 100 steps each
bm.pdf(x, t) // marginal density at state x and time t
bm.mean(t) // theoretical mean at time t
bm.variance(t) // theoretical variance at time t
bm.covariogram(s, t) // theoretical covariance Cov(X(s), X(t))
bm.marginal(2) // Distribution instance representing the marginal at time t
bm.lnL(path) // transition log-likelihood of an observed path under this process
bm.seed(42) // seed the PRNG for reproducible paths; returns the instanceAvailable processes:
| Class | Description |
|---|---|
ran.process.AR1(phi, sigma) |
First-order autoregressive process; stationary for |φ| < 1 |
ran.process.BrownianMotion(mu, sigma, dt) |
Brownian motion with drift; exact discrete-time sampler |
ran.process.OrnsteinUhlenbeck(theta, mu, sigma, dt) |
Mean-reverting process; exact discrete-time sampler |
ran.process.GeometricBrownianMotion(mu, sigma, dt) |
Multiplicative Brownian motion; log-normal increments |
ran.process.BrownianBridge(sigma, T, dt) |
Brownian bridge pinned to 0 at time T |
ran.process.Poisson(lambda, dt) |
Counting process with Poisson(λ·dt) increments per step |
ran.mc provides Markov chain Monte Carlo sampling and convergence diagnostics for targets whose density is known only up to a normalizing constant.
const rwm = new ran.mc.RWM({ logDensity, config: { dim: 1 } }) // logDensity: unnormalized log target density
rwm.warmUp() // tune proposal step size and the thinning interval
rwm.sample(null, 1000) // draw 1000 (thinned) samples
rwm.statistics() // per-dimension { mean, std, cv } since the last reset
rwm.ar() // acceptance rate over the most recent config.arWindow iterations (sliding window)
rwm.ac() // autocorrelation vs. lag, per dimension
rwm.ess() // effective sample size, per dimension (Geyer's initial positive monotone sequence estimator, IPSM)
rwm.state() // snapshot: { x, samplingRate, internal } for resuming a chain
rwm.seed(42) // seed the PRNG for reproducible chains; returns the instanceAvailable samplers:
| Class | Description |
|---|---|
ran.mc.RWM({ logDensity, config, initialState }) |
Random-walk Metropolis-Hastings sampler with Robbins-Monro step-size adaptation during warm-up |
ran.mc.AdaptiveMetropolis({ logDensity, config, initialState }) |
Full-covariance adaptive Metropolis (Haario-Saksman-Tamminen 2001); adapts the joint proposal covariance from the chain's own history during warm-up, then freezes it for sampling |
ran.mc.Gibbs({ conditionals, config, initialState }) |
Component-wise (systematic-scan) Gibbs sampler; draws each dimension directly from its full conditional, so every iteration is accepted (ar() is always 1.0). Each conditional is called as conditionals[d](x, rng) — seed() only reproduces a conditional's draws if it consumes rng.next() for its own randomness instead of an independently-seeded generator |
ran.mc.HMC({ logDensity, gradLogDensity, config, initialState }) |
Hamiltonian Monte Carlo sampler: proposes distant moves via a leapfrog integrator over config.pathLength steps of size config.stepSize, with Metropolis accept/reject on the augmented (position, momentum) system; step size is adapted during warm-up via Robbins-Monro dual averaging and jittered per iteration to avoid periodicity artifacts. config.metric adapts a Euclidean mass matrix during warm-up: 'diag' (default) estimates a per-dimension variance; 'dense' estimates the full covariance matrix via Matrix.ldl() |
ran.mc.MALA({ logDensity, gradLogDensity, config, initialState }) |
Metropolis-Adjusted Langevin Algorithm: proposes a single gradient-informed Langevin step per iteration (x' = x + (stepSize² / 2) · ∇log p(x) + stepSize · z), with a Metropolis-Hastings correction for the proposal's asymmetry; step size is adapted during warm-up via batch Robbins-Monro toward the MALA-optimal 0.574 acceptance rate |
ran.mc.NUTS({ logDensity, gradLogDensity, config, initialState }) |
No-U-Turn Sampler (Hoffman & Gelman 2014): extends HMC with a doubling-tree trajectory that automatically stops at a U-turn, eliminating the need to hand-tune pathLength; the transition is selected via slice sampling over the tree, and step size is adapted during warm-up via the same Robbins-Monro dual averaging as HMC. config.metric adapts a Euclidean mass matrix during warm-up: 'diag' (default) estimates a per-dimension variance; 'dense' estimates the full covariance matrix via Matrix.ldl() |
ran.mc.Slice({ logDensity, config, initialState }) |
Coordinate-wise slice sampler (Neal 2003) using stepping-out and shrinkage; no proposal tuning or gradient required, interval width w is adapted per dimension during warm-up, and ar() is always 1.0 |
ran.mc.ParallelTempering({ logDensity, ...options }) |
Parallel Tempering / Replica Exchange MCMC (Geyer 1991) for multimodal targets; not a subclass of MCMC, it coordinates an array of independent replica samplers (default RWM) at descending inverse temperatures, periodically swapping adjacent replicas' positions so the cold (β = 1) replica inherits the hot replicas' mode-crossing moves |
Work through these questions in order — each one narrows the field:
- Do you know the full conditional distribution of every parameter, given the rest? →
Gibbs. Every draw comes directly from an exact conditional, so there's no accept/reject step andar()is always 1. - Is the target univariate and log-concave on a known finite bracket? →
ARS. Every draw is an exact, independent sample from the target — no warm-up, no burn-in, and no chain to check for convergence. - Do you have the gradient of the log-density? →
NUTSis the default (self-tuning trajectory length via the doubling-tree/U-turn criterion, plus the sameconfig.metricEuclidean mass-matrix adaptation asHMC);HMCif you want the trajectory length under direct control;MALAif a single Langevin step per iteration is enough. BothNUTSandHMCacceptconfig.metric('diag'by default,'dense'for correlated parameters), adapted during warm-up, so a poorly-scaled or correlated target mixes well under either. - No gradient — low or high dimensional? →
Slicefor low-d (no proposal to tune, interval width self-adapts);RWMorAdaptiveMetropolisfor higher-d (AdaptiveMetropoliswhen parameters are correlated, since it learns the full proposal covariance instead of a diagonal one).
Gibbs, NUTS, HMC, MALA, Slice, RWM, and AdaptiveMetropolis all explore the target locally around their current state, so a target with well-separated modes can trap any of them in a single mode. ParallelTempering is the library's remedy: wrap your chosen sampler in it whenever you suspect (or know) the target is multimodal. (ARS sidesteps the question entirely — it isn't a Markov chain, and only accepts log-concave, hence unimodal, targets in the first place.)
Whichever sampler you pick, use runChains with gelmanRubin to check convergence — no signal computable from a single chain can distinguish "converged" from "stuck".
ran.mc.gelmanRubin(samples, maxLength) computes the R-hat convergence diagnostic across two or more independent chains (each an array of states returned by sample()):
const chain1 = new ran.mc.RWM({ logDensity, config: { dim: 1 } }).seed(1)
chain1.warmUp()
const chain2 = new ran.mc.RWM({ logDensity, config: { dim: 1 } }).seed(2)
chain2.warmUp()
const rHat = ran.mc.gelmanRubin([chain1.sample(null, 500), chain2.sample(null, 500)])
// rHat[0] → R-hat vs. iteration count for dimension 0; values near 1 indicate convergenceran.mc.runChains(Sampler, samplerOptions, runOptions) mechanizes the pattern above: it constructs multiple independently-seeded chains of the given MCMC subclass, warms up and samples each, and returns their samples plus the gelmanRubin() diagnostic in one call — the recommended workflow for gating MCMC convergence (single-chain diagnostics cannot distinguish "converged" from "stuck"):
const { samples, rhat } = ran.mc.runChains(ran.mc.RWM, { logDensity, config: { dim: 1 } })
// samples → one sample array per chain (2 chains, seeded [1, 2] by default)
// rhat → the gelmanRubin() diagnostic across those chainsFor multimodal targets that a single chain cannot mix across, ran.mc.ParallelTempering runs several replicas at a geometric ladder of inverse temperatures and swaps their positions so the cold replica benefits from the hot replicas' freer exploration:
const pt = new ran.mc.ParallelTempering({ logDensity, nReplicas: 4, tempMax: 100 })
// or: { logDensity, temperatures: [1, 0.5, 0.25, 0.1] } for an explicit ladder
pt.warmUp(null, 20) // warms up every replica independently, no swaps
const samples = pt.sample(null, 2000) // lockstep sampling with swap proposals; returns the cold replica's draws
pt.swapRate() // accepted/attempted swap fraction per adjacent replica pairranjs signals an unusual result through one of four channels, chosen by the kind of situation:
| Situation | What you get |
|---|---|
| Invalid input — missing/NaN parameters, broken constraints, wrong arity, mismatched dimensions | a thrown Error |
| A valid query whose answer is mathematically undefined (e.g. the mean of a Cauchy distribution) | NaN |
| A valid query whose answer diverges (e.g. the variance of a Pareto with shape ≤ 2, any moment of a Lévy) | Infinity (or -Infinity) |
| A correct value that simply equals zero (e.g. a density evaluated outside the support) | 0 |
Functions never return undefined to mean "failed" or "does not exist" — numeric results stay numbers (NaN/Infinity), and genuine misuse throws. NaN and Infinity are kept distinct on purpose: NaN means no value exists, Infinity means the value grows without bound. This mirrors the conventions of SciPy and R.
const ran = require('ranjs')
new ran.dist.Cauchy(0, 1).mean() // => NaN (undefined moment)
new ran.dist.Pareto(1, 2).variance() // => Infinity (divergent moment)
new ran.dist.Normal(0, 1).pdf(-Infinity) // => 0 (outside support)ranjs targets ≤ 1e-14 relative error for all public outputs in non-degenerate parameter regions. Outputs involving deeply composed operations (quantile inversion, extreme parameter regimes) have a documented floor of ~1e-12, looser still for a handful of quantiles computed by numerical root-finding or near-boundary asymptotics (see below).
All reference values in test/dist-cases-continuous.js and test/dist-cases-discrete.js are sourced from external tools — mpmath at mp.dps = 50, scipy.stats, or Wolfram Alpha — never computed from ranjs itself. Use scripts/gen-dist-refs.py to generate reference values when adding a new distribution, and verify at least one value per distribution against an independent source. pdf, cdf, and pmf reference-value assertions enforce 1e-14 relative tolerance by default; distributions that cannot reach 1e-14 in specific regimes use 1e-12 with an explanatory comment.
All 31 discrete distributions are verified against mpmath references at 50 decimal places. BetaBinomial and NegativeHypergeometric sit at the ~2e-14 float64 arithmetic floor. The following distributions cap at 1e-12 at certain parameter settings: Binomial, Hypergeometric, NegativeBinomial, Poisson. Skellam caps at 1e-12 for most parameter settings, loosening to 1e-11 at highly asymmetric, large means (e.g. mu1=5000, mu2=1) where log-domain cancellation grows with mu1.
All 115 continuous distributions are likewise verified against mpmath references at 50 decimal places (three parameter sets each). pdf/cdf cap at 1e-12–1e-13 at certain parameter settings for: Bates, IrwinHall, Levy, NoncentralBeta, NoncentralChi, NoncentralT, DoublyNoncentralT, SkewNormal, Rice, Tweedie, and R. Quantiles with a closed-form or Halley-refined inverse round-trip to 1e-14; those computed by numerical root-finding (BaldingNichols, Bates, BetaPrime, Davis, FisherZ, Muth, NoncentralChi2, NoncentralF, DoublyNoncentralChi2, DoublyNoncentralT, SkewNormal, Student's t/z, UniformProduct, R) round-trip to ~1e-13–1e-10, and BenktanderII's near-boundary asymptotic branch (b → 1) to ~1e-9.
The summary-statistics functions in src/location/, src/dispersion/, src/shape/, and src/dependence/ have an analogous precision gate in test/precision-summary-stats.js, generated by scripts/precision-refs-summary-stats.py from mpmath references at mp.dps = 50 and checked via scripts/eval-summary-stats.js (see that script's module docstring for the full methodology).
The stochastic processes in src/process/ have the same gate in test/precision-process.js, generated by scripts/precision-refs-process.py. All 9 processes with a closed-form time-t marginal — AR1, BrownianBridge, BrownianMotion, CompoundPoisson, CoxIngersollRoss, GeometricBrownianMotion, OrnsteinUhlenbeck, Poisson, and RandomWalk — are verified across three parameter sets × three times × five interior points each. Every reference gates three independent code paths: pdf(x, t), marginal(t).pdf(x), and marginal(t).cdf(x). Seven of the nine hold at 1e-14 with no exception; RandomWalk at p = 0.3 caps at 3e-14 (log-gamma ULP amplification at t = 30) and CompoundPoisson's pdf at 6e-14 (its marginal is a Tweedie, whose density is an infinite series — its cdf still holds at 1e-14). Each marginal law is re-derived from the process's SDE or update rule inside the generator, which self-checks those derivations against externally sourced values before emitting a single literal.
The tests above pin exact values at a curated set of points; docs/accuracy.md instead reports measured error bounds — in ULP against mpmath at mp.dps = 50 — from a differential-testing harness that sweeps randomized inputs across a documented parameter and argument domain for each special function and distribution pdf/pmf/cdf.
A bound in that table holds only over the domain its row states as swept — it is not a claim about behavior outside that domain. A function or distribution the harness has not yet swept appears as not yet measured rather than being left out, and a known accuracy defect is reported with its actual bad value and a link to the tracking issue rather than rounded away. Treat an unmeasured entry as an open question, not a passing one.
CI fails the build job if dist/ranjs.min.js exceeds 350 KiB raw (358400 bytes) or 96 KiB gzipped (98304 bytes), checked with plain wc -c (raw) and gzip -9 -c | wc -c (gzipped) against the RAW_BUDGET/GZIP_BUDGET values in .github/workflows/ci.yml — no size-analysis dependency needed. Both budgets carry roughly 25-30% headroom over the current size, enough to catch an accidental tree-shaking break (e.g. a helper imported unconditionally by every distribution, defeating the per-distribution subpath exports) without needing a bump on every ordinary PR. When growth is intentional — a new distribution, process, or MCMC sampler — raise the relevant budget in the same PR and note the new size here.
The gzipped figure is the one that reflects real download cost, since npm/CDN delivery is compressed. The full bundle's ~74 KiB gzipped is not what most consumers actually ship, though: this project's own guidance (see ESM — single distribution import) is to import individual distributions via their subpath exports. Those subpath builds are minified with terser, and gzip to roughly 10 KiB each (e.g. dist/beta.esm.js and dist/normal.esm.js are both ~10 KiB gzipped) — well inside the range the community treats as unremarkable for a library that does real work (for reference, lodash's full, non-tree-shaken build is ~24 KiB gzipped; moment.js's oft-cited ~67 KiB gzipped is the anchor most commonly pointed to as "too big" and a reason to look elsewhere).
Full API reference and distribution catalogue: https://synesenom.github.io/ran/
The public API is considered stable: a breaking change first ships a deprecation warning and stays in place for at least one released minor version before removal, documented in CHANGELOG.md as it happens. Bug reports and security issues are triaged through GitHub issues. Releases are batched, not continuous — changes land on main under [Unreleased] and ship together in a dedicated release PR once enough have accumulated, not on every merge. There is no guaranteed response-time SLA.