Benchmarks, dashboards and experiments keep asking the same questions of a list of numbers: its mean and spread, its median and 99th percentile, the line through a set of points and how well it fits. When the data fits in memory, an exact answer is simpler and faster than an estimate. This library gives those answers, with the quantile and regression conventions of NumPy and SciPy.

ckms estimates quantiles over an unbounded stream in bounded memory, and hdr buckets values at a fixed relative precision; use those when the data does not fit in memory or does not end. Stats.Running covers the case between: a mean and a variance over a stream, in constant memory and one pass. It uses Welford's update (1962): on each new x it keeps the count n, the mean m and M2, the sum of squared deviations from the mean, with m' = m + (x − m)/n and M2' = M2 + (x − m)(x − m'), and the variance is M2/(n − 1). No large sums are ever subtracted, so it avoids the cancellation described under Conventions.

Absent data and invalid data #

The library reports the two differently. A question about an empty sample, such as the mean of no numbers, has no answer and returns None. An input that cannot produce a statistic, such as a NaN among the observations or a quantile outside [0, 1], raises Invalid_argument. Returning a number in either case would let a wrong answer travel downstream unnoticed.

Sample.t is abstract because it carries an invariant: it is sorted, and every observation in it is finite. A quantile read from an unsorted array looks like a quantile but is not one.

Checked against NumPy and SciPy #

The interface promises Hyndman and Fan type 7 quantiles, NumPy's inverted_cdf for weighted quantiles, and the four numbers scipy.stats.linregress reports. test/interop/numpy and test/interop/scipy replay committed traces from NumPy 2.0.2 and SciPy 1.13.1 and compare. REGEN=1 dune build @traces regenerates them; the tests themselves need no Python.

The traces also record two deliberate differences. For a two-point fit SciPy reports a standard error of 0, and this library reports infinity. The slope's standard error divides the residual sum of squares by n − 2 degrees of freedom; with two points the line passes through both, so the estimate is 0/0 and says nothing about the spread, which infinity reports and 0 hides. For the variance of one observation NumPy returns nan, and this library returns None.

Install #

$ opam install stats

If opam cannot find the package, add the overlay repository first:

$ opam repo add samoht https://tangled.org/gazagnaire.org/opam-overlay.git
$ opam update
$ opam install stats

Conventions #

Every regression is one weighted least-squares fit underneath. Libraries disagree on two choices, so both are stated here.

The fit is computed from centred sums. The raw-sum form computes the slope as (n Σxy − Σx Σy) / (n Σx² − (Σx)²), which subtracts two large numbers to get a small one; with x near 2000 (a year), x² is near 4,000,000 and the difference cancels most of the 53-bit mantissa. The centred form first subtracts the means and computes Σ(x − x̄)(y − ȳ) / Σ(x − x̄)², whose terms stay small.

Quantiles are linearly interpolated between the two neighbouring order statistics, Hyndman and Fan type 7, the default of NumPy, R and pandas. For a sorted sample x[0] ≤ ... ≤ x[n − 1] and a probability p, let h = (n − 1)p and k = ⌊h⌋; the quantile is x[k] + (h − k)(x[k + 1] − x[k]). In the example below the sorted sample is 1 1 2 3 4 5 6 9, so p = 0.99 gives h = 6.93 and 6 + 0.93 × 3 = 8.79. Taking the nearer of the two instead discards the fraction and moves the answer by a whole rank, which on a short series separates a median from its neighbour.

The slope's standard error is the textbook ordinary least-squares one, which assumes independent observations. Spatial and time-series data usually are not independent, so on that data it understates the true spread; the interface documentation says so.

Usage #

# let s = Stats.Sample.of_list [ 3.; 1.; 4.; 1.; 5.; 9.; 2.; 6. ]
val s : Stats.Sample.t = <abstr>
# Stats.Sample.median s
- : float option = Some 3.5
# Stats.Sample.quantiles s [ 0.25; 0.75; 0.99 ]
- : float list option = Some [1.75; 5.25; 8.79]
# Stats.Sample.mean s, Stats.Sample.stddev s
- : float option * float option = (Some 3.875, Some 2.74837614393871288)

A line, and what to do with it:

# let f = Option.get (Stats.fit [ (1., 2.); (2., 4.); (3., 5.); (4., 9.) ])
val f : Stats.fit =
  {Stats.n = 4; slope = 2.2; intercept = -0.5; r = 0.964763821237732189;
   slope_stderr = 0.424264068711928177}
# f.Stats.slope, f.Stats.intercept
- : float * float = (2.2, -0.5)
# Stats.r_squared f
- : float = 0.930769230769230904
# Stats.predict f 5.
- : float = 10.5

Pearson correlation measures how close the points are to a straight line, and Spearman correlation only whether y rises with x:

# let curve = [ (1., 1.); (2., 4.); (3., 9.); (4., 16.); (5., 25.) ]
val curve : (float * float) list =
  [(1., 1.); (2., 4.); (3., 9.); (4., 16.); (5., 25.)]
# Stats.pearson curve
- : float = 0.981104910251592854
# Stats.spearman curve
- : float = 1.

And over a stream, in constant memory:

# let r = Stats.Running.v ()
val r : Stats.Running.t = <abstr>
# List.iter (Stats.Running.add r) [ 1.; 2.; 3.; 4. ]
- : unit = ()
# Stats.Running.mean r, Stats.Running.stddev r
- : float option * float option = (Some 2.5, Some 1.2909944487358056)

A test against a null distribution drawn by resampling, with the false discovery rate held over many such tests:

# let null = [| 1.; 2.; 5.; 3.; 0.5; 6.; 1.; 2.; 3. |]
val null : float array = [|1.; 2.; 5.; 3.; 0.5; 6.; 1.; 2.; 3.|]
# Stats.empirical_p Upper 5. null
- : float = 0.3
# Stats.q_values [ 0.01; 0.04; 0.03; 0.2 ]
- : float list = [0.04; 0.0533333333333333368; 0.0533333333333333368; 0.2]
# Stats.wilson 5 10
- : float * float = (0.236589593615487315, 0.76341040638451263)

API #

Stats (lib/stats.mli) is the whole interface: Sample for order statistics and descriptives, fit/weighted/predict/r_squared for least squares, pearson/spearman for correlation, Running for one pass over a stream, and empirical_p/q_values/permutation_null/bootstrap/ percentile_interval/wilson for tests and intervals by resampling.

License #

ISC. See LICENSE.md.

Clone this repository

stats #