math/stats
Package math/stats.
Statistics over lists of any numeric type, and streaming accumulators for numbers that never all exist at once. Written in slang over math.
import "math/stats";
let xs = [12.5, 9.0, 14.25, 11.0];
stats.mean(xs); // 11.6875
stats.median(xs); // 11.75
stats.stddev(xs); // sample (n - 1); pstddev is population
stats.quantile(xs, 0.9); // R type 7 / NumPy's default
stats.correlation(xs, ys); // Pearson's r; spearman for ranks
let fit = stats.linear_regression(xs, ys); // fit.slope, fit.intercept, fit.r2, fit.at(x)
let h = stats.histogram(xs, 10); // h.edges, h.counts
// Streaming, O(1) per value, mergeable across tasks:
let lat = stats.running(); // count, mean, variance, min, max
lat.add(elapsed_ms);
let p = stats.digest(); // quantiles of a stream (t-digest)
p.add(elapsed_ms);
p.quantile(0.99);
stats.normal(100.0, 15.0).cdf(130.0); // 0.9772498680518208
stats.normal(0.0, 1.0).quantile(0.975); // 1.9599639845400538
stats.poisson(3.0).pmf(5);- Center and spread:
sum(compensated),product,mean,weighted_mean,geometric_mean,harmonic_mean,median,mode,min,max,spread,variance/pvariance,stddev/pstddev,mad,zscores,skewness,kurtosis(excess). - Order:
quantile,percentile,quartiles,iqr,ranks;sorted_floatsandquantile_sortedto sort once and ask many times. - Two lists:
covariance,correlation,spearman,linear_regression. - Streaming:
running()(Welford),ewma(alpha)andewma_half_life(n),digest()(a merging t-digest: p99 and p999 stay accurate in about 100 centroids). Each hasadd, andrunninganddigestmerge. - Distributions:
normal,log_normal,uniform,exponential(pdf,cdf,quantile),poissonandbinomial(pmf,cdf).
A result with no value for its input -- the mean of no numbers, the variance of one, a correlation with a constant list -- is NaN, which propagates and tests true with math.is_nan; functions that return one of the elements (min, max, mode) return an opt of the element type. Sums are compensated (Neumaier) and variances are two-pass around the mean, so neither drifts on long or offset data: stats.sum([0.1, 0.2, 0.3]) is 0.6.
API#
The normal (Gaussian) distribution with mean mu and stddev sigma.
Through erfc, which keeps precision far into the lower tail where 1 + erf would round to 0.
The x with cdf(x) = p; -inf at 0, inf at 1, NaN outside [0, 1].
The distribution of exp(X) for X normal(mu, sigma): response times and sizes, which are positive and skewed right.
The continuous uniform distribution on [lo, hi].
The waiting time between events arriving at rate per unit time.
1 - exp(-rate x), through expm1 so a tiny x keeps its precision.
The number of events in an interval where lambda are expected.
P(X = k): directly while k! is exact (k <= 20) and lambda^k cannot overflow, otherwise in log space, where lgamma costs a few ulps.
P(X <= k).
The number of successes in n independent trials, each succeeding with probability p.
P(X = k): directly while the binomial coefficient is exact (n up to 1000, as math.binomial computes it), otherwise in log space.
P(X <= k).
The sum of xs, compensated: the exact rounding error of each addition (Knuth's TwoSum) is carried and added back at the end. 0 for an empty list; an infinity or NaN in xs gives the plain sum's inf or NaN. Two independent running sums take alternate elements, so consecutive additions do not wait on each other: one chain held a sum to the latency of a float add per element, about half the speed of a plain loop. Every list function here sums the same way, in place -- they used to collect the terms in a list first and sum that, which for variance on a million values made the pushes 90% of the time. A for x in xs loop, not a while over indices: the for loop takes its safepoint once around the loop and reads elements unchecked, where a while loop brackets every iteration and bounds-checks every index -- twice the time per element.
The product of xs; 1 for an empty list.
The arithmetic mean; NaN for an empty list.
The mean with weights ws (one per value); NaN when the lengths differ, the list is empty, or the weights sum to 0.
The geometric mean, exp(mean(log x)): for growth rates and ratios. 0 if any value is 0; NaN for an empty list or a negative value.
The harmonic mean, n / sum(1/x): for averaging rates. 0 if any value is 0; NaN for an empty list or a negative value.
The smallest element; none for an empty list. NaN never compares smaller, so a float list's NaNs are passed over.
The largest element; none for an empty list. NaNs are passed over.
max - min, as a float; NaN for an empty list.
xs as floats, sorted ascending (NaN first, as sort.floats orders it): what median and quantile work on. Compute it once and use the _sorted variants when asking for several quantiles of one list.
The middle value, or the mean of the two middle values; NaN for an empty list. Copies and sorts: O(n log n).
The q-quantile (q in [0, 1]) by linear interpolation between the two nearest ranks -- type 7 in Hyndman and Fan, the default of R, NumPy and Excel's PERCENTILE.INC. NaN for an empty list or q outside [0, 1].
quantile on a list already sorted ascending (sorted_floats), without copying it.
quantile with p in percent: percentile(xs, 95) is quantile(xs, 0.95).
The three quartiles, by quantile's method.
The interquartile range, q3 - q1.
The most frequent value, the smallest of them on a tie; none for an empty list. Floats compare exactly, so 0.1 + 0.2 and 0.3 are different values.
The sample variance (divided by n - 1): the variance of the population the values were drawn from, estimated without bias. NaN for fewer than two values. pvariance is the variance of exactly these values.
The population variance (divided by n); NaN for an empty list.
The sample standard deviation, sqrt(variance).
The population standard deviation, sqrt(pvariance).
The median absolute deviation, median(|x - median(xs)|): a spread that one wild value cannot move, unlike stddev. Unscaled; multiply by 1.4826 to estimate a normal distribution's stddev.
Each value's distance from the mean in population standard deviations (SciPy's zscore). All NaN when every value is the same.
The skewness, the population third standardized moment (Fisher-Pearson g1): 0 for a symmetric distribution, positive with a long right tail. NaN for an empty list or when every value is the same.
The excess kurtosis, m4 / m2^2 - 3: 0 for a normal distribution, positive for heavier tails. NaN for an empty list or constant values.
The sample covariance of xs and ys (divided by n - 1); NaN when the lengths differ or there are fewer than two pairs.
Pearson's correlation coefficient r, in [-1, 1]; NaN when the lengths differ, there are fewer than two pairs, or either list is constant.
Spearman's rank correlation: Pearson's r of the values' ranks, so any monotonic relationship scores 1, not just a straight line. Tied values share the mean of their ranks.
Each value's rank, 1 for the smallest; tied values share the mean of the ranks they span ([10, 20, 20, 30] ranks as [1, 2.5, 2.5, 4]).
The fitted line at x.
The least-squares line y = intercept + slope * x through the points. Every field is NaN when the lengths differ, there are fewer than two points, or every x is the same.
Counts xs into bins equal-width bins over [min, max]; the last bin also holds max itself. NaNs are not counted. bins below 1 is taken as 1, and a list of one repeated value gets one bin of width 1 around it.
Count, mean, variance, min and max of a stream (Welford's algorithm): O(1) per value, numerically stable, no samples kept. let r = stats.running(); r.add(12.5).add(9.0); r.mean(); r.stddev(); r.max();
Adds one value; returns the accumulator so adds chain.
Folds another accumulator's values into this one (Chan et al.), as if every value added to other had been added here.
NaN before the first value.
mean * count.
The sample variance; NaN before the second value.
The population variance; NaN before the first value.
NaN before the first value.
NaN before the first value.
An exponentially weighted moving average: each value moves the average a fraction alpha of the way toward it, so recent values count most and old ones fade. The first value is taken as is.
alpha in (0, 1]: 0.1 forgets slowly, 0.5 quickly. Clamped into range.
The average whose weight halves every n values: an alpha of 1 - 0.5^(1/n), for a decay stated in samples rather than as a fraction.
NaN before the first value.
Streaming quantiles: a merging t-digest (Dunning and Ertl, 2019). The stream is summarized as at most about compression weighted centroids, small near the extremes and large in the middle, so p99 and p999 stay accurate (to well under 1% of rank at the default 100) where a fixed-width histogram would need to know the range up front. Digests merge, so each worker can keep its own. let d = stats.digest(); d.add(latency_ms); d.quantile(0.5); d.quantile(0.99);
A digest at the default compression, 100.
A digest keeping about compression centroids: more is more accurate and larger (the memory is about 16 bytes a centroid). At least 20.
Adds one value. NaN is ignored: it has no rank.
Adds a value standing for w observations (w > 0).
Folds another digest into this one, as if its values had been added here (other is not changed).
The number of observations, as the sum of their weights.
NaN when empty.
NaN when empty.
The estimated q-quantile, q in [0, 1]: min at 0, max at 1, linear between centroid centers. NaN when empty or q is outside [0, 1].
The estimated fraction of observations at or below x, in [0, 1]; NaN when empty.