Science.js

Science.js is a JavaScript library for scientific and statistical computing: summary statistics, distributions, kernel density estimation, clustering, smoothing and linear algebra, with methods similar to those provided by R.

Installation

npm install science

Science.js is an ES module, and works in Node.js 18 or later and in modern browsers. Import the functions you need:

import { mean, loess } from "science";

CommonJS code can load it on Node.js 20.19+ or 22.12+, where require supports ES modules:

const { mean, loess } = require("science");

In a browser, you can instead load dist/science.min.js with a <script> tag, which defines a global science with the same functions, such as science.mean.

Type definitions for TypeScript are included.

API Reference

• Summary statistics
• Distributions
• Kernel density estimation
• Smoothing
• Clustering
• Distances
• Linear algebra
• Other

Functions take arrays (or typed arrays) of numbers, and never modify their inputs. Optional settings are passed as an options object, and an option that the function does not have throws, so that typos are caught. Invalid input throws a TypeError (for a value of the wrong type, or an unknown option) or a RangeError (for a value out of range, such as arrays of different lengths).

Summary statistics

These functions take an array of numbers. NaN and infinities are allowed, and propagate as in ordinary arithmetic.

# mean(values)

Returns the arithmetic mean of values, or NaN if it is empty. It is computed incrementally, so that large values do not overflow.

mean([1, 2, 3, 4]); // 2.5

# median(values)

Returns the median of values, or NaN if it is empty.

median([3, 1, 4, 2]); // 2.5

# quantile(values, p)

Returns the p-quantile of values, where p is between 0 and 1, using linear interpolation between the sorted values (type 7 in R's quantile, its default). Returns NaN if values is empty.

quantile([1, 2, 3, 4], 0.25); // 1.75

# quantiles(values, probabilities)

Returns the quantile of values for each of the probabilities, sorting values only once.

quantiles([1, 2, 3, 4], [0, 0.5, 1]); // [1, 2.5, 4]

# iqr(values)

Returns the interquartile range of values: the difference between its 0.75 and 0.25 quantiles.

iqr([1, 2, 3, 4]); // 1.5

# variance(values)

Returns the unbiased sample variance of values (with denominator n − 1), or NaN if there are fewer than two values. It uses Welford's algorithm, which is accurate even when the values are large compared with their spread.

variance([1, 2, 3, 4]); // 1.6666666666666667

# standardDeviation(values)

Returns the sample standard deviation of values, the square root of its variance.

standardDeviation([2, 4, 4, 4, 5, 5, 7, 9]); // ≈ 2.1380899

# modes(values)

Returns all of the most frequent values in values, in the order they first appear, or an empty array if values is empty. Values are compared as in a Map, so 1 and "1" are different values.

modes([1, 2, 2, 3]); // [2]
modes([3, 3, 1, 1, 2]); // [3, 1]

Distributions

# normal([options])

Returns a normal distribution with the given options:

• mean: the mean, which defaults to 0.
• sd: the standard deviation, which must be positive and defaults to 1.

The distribution is an object with these properties and methods:

• mean and sd: the parameters.
• pdf(x): the probability density at x.
• cdf(x): the cumulative probability at x, with a relative error of about 1e-14 within 8 standard deviations of the mean, and about 1e-13 even far into the tails.
• quantile(p): the inverse of cdf, for a probability p between 0 and 1.
• sample([random]): a random number drawn from the distribution, using the Marsaglia polar method. random is a function returning uniform random numbers in [0, 1), which defaults to Math.random; pass a seeded generator for reproducible samples.

const heights = normal({mean: 170, sd: 10});
heights.cdf(180); // ≈ 0.8413447
heights.quantile(0.975); // ≈ 189.5996398
normal().pdf(0); // ≈ 0.3989423

# erf(x)

Returns the error function at x, accurate to about 1e-15 relative error.

erf(1); // ≈ 0.8427008

# erfc(x)

Returns the complementary error function, 1 − erf(x), computed directly so that it stays accurate for large x, where 1 − erf(x) would round to 0.

erfc(10); // ≈ 2.0884876e-45

Kernel density estimation

# kde(sample[, options])

Returns the kernel density estimate for the array sample, as a function that takes a value x and returns the estimated density there. Its bandwidth property is the bandwidth used. The options are:

• kernel: the kernel function, which defaults to kernels.gaussian.
• bandwidth: a positive number, or a function that is passed the sample and returns one, such as one of the bandwidth selectors. Defaults to bandwidths.nrd0, as in R's density.

const density = kde([1, 2, 3, 4, 5], {bandwidth: 1});
density(3); // ≈ 0.1981731
density.bandwidth; // 1

# kernels

Kernel functions, each of which takes a value u and returns its weight. All of them integrate to 1, and all but the Gaussian are zero outside [−1, 1].

• kernels.uniform: ½.
• kernels.triangular: 1 − |u|.
• kernels.epanechnikov: ¾(1 − u²).
• kernels.quartic: (15/16)(1 − u²)².
• kernels.triweight: (35/32)(1 − u²)³.
• kernels.gaussian: e−u²/2 / √(2π).
• kernels.cosine: (π/4) cos(πu/2).

kernels.epanechnikov(0); // 0.75

# bandwidths

Bandwidth selectors for Gaussian kernels, matching R's bw.nrd0 and bw.nrd. Each takes an array of at least two values and returns a bandwidth.

• bandwidths.nrd0: Silverman's rule of thumb, 0.9 × min(s, IQR / 1.34) × n−1/5, where s is the standard deviation and IQR the interquartile range. If that minimum is 0, it falls back to s, then to the absolute value of the first value, then to 1.
• bandwidths.nrd: Scott's variation, 1.06 × min(s, IQR / 1.34) × n−1/5. Unlike nrd0, it is 0 when the interquartile range is 0, such as when most values are repeated.

bandwidths.nrd0([1, 4, 3, 2]); // ≈ 0.7635139
bandwidths.nrd([1, 4, 3, 2]); // ≈ 0.8992498

Smoothing

# loess(x, y[, options])

Smooths the points (x[i], y[i]) with LOESS, returning the smoothed y value for each point. x and y must be arrays of finite numbers of the same length, and x must be sorted in ascending order (tied values are allowed). It follows R's lowess, and is accurate even for large x values such as timestamps. The options are:

• bandwidth: the fraction of the points used for each local regression, which must be greater than 0 and at most 1, and defaults to 0.3. Lower values follow the data more closely; higher values give a smoother curve. With three or more points, floor(bandwidth × n) must be at least 2: for example, the default bandwidth needs at least 7 points.
• robustnessIterations: the number of iterations after the initial fit that reduce the weight of points with large residuals, making the curve less sensitive to outliers. Must be a non-negative integer, and defaults to 2; use 0 for a plain locally weighted regression.
• accuracy: a relative tolerance, which defaults to 1e-12. As in R's lowess, a local regression whose x values spread less than accuracy × the range of x fits a flat line, and the robustness iterations stop once the median residual falls below accuracy × the mean absolute y value, so the results do not depend on the units of x and y.
• weights: an array with a non-negative weight for each point, scaling its influence. Points with a weight of 0 have no influence, but are still given a smoothed value from their neighbours. At least one weight must be positive. Defaults to all 1.

loess([1, 2, 3, 4, 5], [2, 4, 6, 8, 10], {bandwidth: 0.6}); // [2, 4, 6, 8, 10]

Clustering

# kmeans(vectors, options)

Clusters the array of vectors (arrays of numbers, all the same length) with k-means, choosing the initial centroids with k-means++ and iterating until the assignments stop changing. The options are:

• k: the number of clusters, which is required. There must be at least k distinct vectors.
• distance: the distance function used to assign vectors to centroids, called as distance(vector, centroid), which defaults to distances.euclidean. Every distance must be finite.
• maxIterations: the maximum number of iterations, which defaults to 1,000.
• random: a function returning uniform random numbers in [0, 1) for the initial centroids, which defaults to Math.random; pass a seeded generator for reproducible results.

Returns an object with these properties:

• assignments: the index of the cluster of each vector.
• centroids: the mean of the vectors in each cluster.
• iterations: the number of iterations performed.

const {assignments} = kmeans([[0], [1], [10], [11]], {k: 2});
assignments[0] === assignments[1]; // true
assignments[0] === assignments[2]; // false

# hcluster(vectors[, options])

Clusters the array of vectors (arrays of the same length) with agglomerative hierarchical clustering, repeatedly merging the two closest clusters. The vectors can contain any values that the distance function accepts, such as booleans with distances.jaccard. The options are:

• linkage: how the distance between two clusters is computed from the distances between their vectors: "single" (the smallest, the default), "complete" (the largest) or "average" (the mean).
• distance: the distance function between vectors, which defaults to distances.euclidean. It must be symmetric, and every distance must be finite.

Returns the root of the binary tree of clusters (a single leaf if there is only one vector). Each node has these properties:

• left and right: the two merged clusters, or null for a leaf.
• height: the linkage distance between left and right, or 0 for a leaf.
• size: the number of vectors in the cluster.
• depth: the height of the node in merges: 0 for a leaf, otherwise 1 more than the larger depth of its children.
• index: for a leaf, the index of its vector in vectors. Merged clusters are numbered from n (the number of vectors) in the order they are merged, as in SciPy.
• centroid: the mean of the vectors in the cluster, or null unless every value in vectors is a finite number.

const root = hcluster([[1], [2], [5]], {linkage: "average"});
root.height; // 3.5
root.left.height; // 1
root.size; // 3

Distances

# distances

Distance functions, each of which takes two arrays of the same length, for use with kmeans and hcluster.

• distances.euclidean: the Euclidean distance.
• distances.manhattan: the Manhattan distance, the sum of the absolute differences.
• distances.chebyshev: the Chebyshev distance, the largest absolute difference.
• distances.minkowski(p): returns the Minkowski distance function of order p. Order 1 is the Manhattan distance, 2 the Euclidean and Infinity the Chebyshev.
• distances.hamming: the Hamming distance, the number of positions at which the arrays differ.
• distances.jaccard: the Jaccard distance between the arrays treated as sets, where a truthy value at position i means that i is a member: 1 − |a ∩ b| / |a ∪ b|, or 0 if both are empty.
• distances.braycurtis: the Bray–Curtis dissimilarity, Σ|ai − bi| / Σ|ai + bi|, for non-negative values.

distances.euclidean([0, 0], [3, 4]); // 5
distances.jaccard([1, 0, 1], [1, 1, 0]); // ≈ 0.6666667

Linear algebra

Vectors are arrays of numbers, and matrices are arrays of rows.

# dot(a, b)

Returns the dot product of the vectors a and b, which must have the same length.

dot([1, 2, 3], [4, 5, 6]); // 32

# cross(a, b)

Returns the cross product of the three-dimensional vectors a and b.

cross([1, 0, 0], [0, 1, 0]); // [0, 0, 1]

# norm(a)

Returns the Euclidean length of the vector a, without overflow or underflow for very large or small values.

norm([3, 4]); // 5

# normalize(a)

Returns the vector in the direction of a with length 1. Throws for a zero vector.

normalize([3, 4]); // [0.6, 0.8]

# transpose(a)

Returns the transpose of the matrix a.

transpose([[1, 2, 3], [4, 5, 6]]); // [[1, 4], [2, 5], [3, 6]]

# multiply(a, b)

Returns the product of the matrix a with b, which is either a matrix (giving a matrix) or a vector (giving a vector). The number of columns of a must equal the number of rows of b.

multiply([[1, 2], [3, 4]], [[5, 6], [7, 8]]); // [[19, 22], [43, 50]]
multiply([[1, 2], [3, 4]], [5, 6]); // [17, 39]

# determinant(a)

Returns the determinant of the square matrix a, computed by LU decomposition with partial pivoting.

determinant([[1, 2], [3, 4]]); // -2

# solve(a, b)

Solves the linear system a x = b for x, where a is a square matrix and b is a vector or a matrix of right-hand sides, by LU decomposition with partial pivoting. Throws a RangeError if a is singular, treating pivots no larger than n × Number.EPSILON × the largest magnitude in a as zero.

// 2x + y = 5 and x + 3y = 10.
solve([[2, 1], [1, 3]], [5, 10]); // [1, 3]

# inverse(a)

Returns the inverse of the square matrix a. Throws a RangeError if a is singular, as for solve. To solve a linear system, use solve instead, which is faster and more accurate.

inverse([[4, 7], [2, 6]]); // [[0.6, -0.7], [-0.2, 0.4]]

# solveTridiagonal(lower, diagonal, upper, rhs)

Solves a tridiagonal system of n linear equations with the Thomas algorithm, where equation i is lower[i − 1] x[i − 1] + diagonal[i] x[i] + upper[i] x[i + 1] = rhs[i]. diagonal and rhs have n elements, and lower and upper have n − 1. No pivoting is done, so the system should be, for example, diagonally dominant.

solveTridiagonal([1, 1], [2, 2, 2], [1, 1], [3, 4, 3]); // [1, 1, 1]

# eigen(a)

Returns the eigendecomposition of the square matrix a, based on JAMA, as an object with these properties:

• values: the real parts of the eigenvalues.
• imaginary: the imaginary parts of the eigenvalues.
• vectors: a matrix whose columns are the eigenvectors, each of unit length.

If a is symmetric, the eigenvalues are real and in ascending order, and the eigenvectors are orthonormal. Otherwise, for a complex conjugate pair of eigenvalues at positions i and i + 1, columns i and i + 1 of vectors hold the real and imaginary parts of the eigenvector for the eigenvalue at position i, scaled together to unit length. In the rare case that the iteration does not converge, an Error is thrown.

const {values, vectors} = eigen([[2, 0], [0, 1]]);
values; // [1, 2]
vectors; // [[0, 1], [1, 0]]

Other

# quadratic(a, b, c[, options])

Returns the real roots of a**x² + b**x + c = 0 in ascending order: two roots, one (a repeated root), or none. If a is 0, the equation is solved as a linear equation. The roots are computed so as to avoid loss of precision. With the option complex set to true, the roots are returned as objects {re, im} and complex roots are included, in ascending order of their imaginary parts.

quadratic(1, 1, -2); // [-2, 1]
quadratic(1, 0, 1); // []
quadratic(1, 0, 1, {complex: true}); // [{re: 0, im: -1}, {re: 0, im: 1}]

# EULER_GAMMA

The Euler–Mascheroni constant, γ ≈ 0.5772156649015329.

# version

The version of Science.js, such as "2.0.0".

Development

To help develop Science.js, you need to have Node.js installed. Then run the following from the root directory of this repository:

npm install
npm test

npm test builds the browser bundles in dist/ and runs the tests. The type definitions in index.d.ts are maintained by hand; the tests check that they declare every export.

Each example in examples/ is a single HTML page that loads D3 from a CDN and Science.js from src/ as ES modules, so they must be served over HTTP. For example, run python3 -m http.server in the root directory and open http://localhost:8000/examples/loess/.