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
- 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.
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.
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/