Statistics
Descriptive statistics of a column or a data matrix: moments, covariance and correlation, ranks, order statistics, quantiles, sorting and weighted forms.
This module describes a column of numbers, or a matrix of columns: its location, spread and shape; the
covariance and correlation between columns; ranks, order statistics, the median, quantiles and sorting;
rank correlations; and weighted versions of the moments. Every definition is written in Hilbert over
einsum, a few elementwise functions and rows/columns, so each one compiles to a single tensor
program.
Conventions. A column is a vector n. A data matrix is n×p: one row per observation, one column per
variable. The singular names work on one column (mean(x), variance(x)) and the plural names on every
column of a matrix at once (means(x), variances(x)). variance, std, covariance and the matrix
forms are the sample statistics, dividing by \(n-1\); the population forms have names of their own
(populationVariance, populationStd) rather than a flag. The higher moments, skewness and kurtosis use
the population moments (dividing by \(n\)). Ranks and order-statistic positions count from zero.
Ties and cost. Ranks and order statistics are built from an n×n comparison matrix, so they cost \(O(n^2)\) in time and memory — comfortable at a few thousand rows, not at a million. Ties are handled by averaging: tied values share a midrank, and an order statistic, median or mode that falls on several equal values is their average.
open Prelude.Statistics
Location
total
total(x) · x column n → scalar
Let \(x = (x_1, \dots, x_n) \in \mathbb{R}^n\) be a column of \(n\) values. Then
\[\operatorname{total}(x) = \sum_{i=1}^{n} x_i .\]The sum of the values.
Intuition The grand sum: one contraction collapses the column’s single index into a scalar. Most location and spread statistics are a total of something divided by a count, so this is the first brick of the module.
Use it when you need an amount rather than a typical value: total revenue, total events in a period, or the denominator of a share (x / total(x) turns counts into proportions). It expects one column. For per-column averages of a matrix use means, and for the sum of a weight vector use weight.
data-summaryaggregation
def total(x) ≔ einsum("i->", x)
count
count(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values (or \(X \in \mathbb{R}^{n \times p}\) a matrix). Then
\[\operatorname{count}(x) = n .\]The number of values \(n\) (for a matrix, the number of rows), read at run time.
Intuition \(n\) is read from the tensor’s shape at run time, not counted from its contents: every row counts, whatever it holds, zeros and NaNs included. For a matrix it is the number of rows, that is, the number of observations.
Use it when a formula needs the sample size: a divisor, the \(n-1\) degrees of freedom, or the \(\sqrt{n}\) in a standard error. It does not skip missing or zero values. To count events, sum a 0/1 flag vector with total instead.
data-summarysample-size-planning
def count(x) ≔ rows(x)
mean
mean(x) · x column n → scalar
Let \(x = (x_1, \dots, x_n) \in \mathbb{R}^n\) be a column of \(n\) values. Then
\[\operatorname{mean}(x) = \bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i .\]The arithmetic mean.
Intuition The balance point of the data: the value around which the deviations sum to zero, and the constant that minimises the sum of squared deviations. Every value pulls on it in proportion to its size, so a single extreme value can move it a long way.
Use it when you summarise a roughly symmetric column, estimate an expected value, or need a centre to subtract. For skewed or outlier-prone data (incomes, latencies) report median as well. Use means for every column of a matrix and weightedMean when rows count unequally.
central-tendencydata-summary
def mean(x) ≔ total(x) / count(x)
centre
centre(x) · x column n → column n
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values with mean \(\bar{x}\). Then
\[\operatorname{centre}(x)_i = x_i - \bar{x} , \qquad i = 1, \dots, n .\]The centred column. Every moment below is a contraction of this with itself.
Intuition Moves the origin to the mean, so each entry says how far above or below typical that value is, and the column now sums to zero. Variances, covariances and the higher moments are all products of centred values, which is why they are built on this.
Use it when you need deviations rather than levels: before a dot product that should measure co-movement rather than size, or as the first step of standardising (zscore divides it by std). For a whole data matrix use centred.
feature-scalingdata-preprocessing
def centre(x) ≔ x − mean(x)
Spread
The sample variance divides by \(n - 1\) and the population variance by \(n\). The two differ by 41% at three values and by a fraction of a percent at a few thousand, so the choice is made by name, visibly at the call site.
\[s^2 = \frac{1}{n-1}\sum_i (x_i - \bar{x})^2\]variance
variance(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n \ge 2\) values with mean \(\bar{x}\). Then
\[\operatorname{variance}(x) = s^2 = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})^2 .\]The sample variance.
Intuition The average squared distance from the mean, in squared units. The mean was fitted to these same values, which pulls the deviations in. Dividing by \(n-1\) instead of \(n\) corrects for that, so the result is an unbiased estimate of the variance of the population the sample came from.
Use it when the column is a sample and you reason beyond it: confidence intervals, t-tests, or comparing variability between groups. It needs \(n \ge 2\), since one value gives zero over zero. For spread in the data’s own units use std. When the column is the entire population, use populationVariance.
spreaduncertainty-quantificationdata-summary
def variance(x) ≔ squares(x) / (count(x) − 1)
std
std(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n \ge 2\) values with mean \(\bar{x}\). Then
\[\operatorname{std}(x) = s = \sqrt{\frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})^2} .\]The sample standard deviation.
Intuition The square root undoes the squaring, so spread is back in the units of the data: a typical distance from the mean. For roughly normal data about 68% of values lie within one standard deviation of the mean and about 95% within two.
Use it when you report variability next to a mean, set thresholds such as “three standard deviations”, or scale features (zscore divides by it). It is the sample form (\(n-1\)), suited to inference. For a purely descriptive spread of a complete population use populationStd. Outliers inflate it, and iqr is the robust alternative.
spreaddata-summaryoutlier-detection
def std(x) ≔ sqrt(variance(x))
populationVariance
populationVariance(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values with mean \(\bar{x}\). Then
\[\operatorname{populationVariance}(x) = \sigma^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^2 .\]The population variance.
Intuition The same average squared distance, but over exactly \(n\) values. It is the right number when the column is the whole population rather than a sample, and it is the maximum-likelihood variance of a normal fit. It is always \((n-1)/n\) times variance.
Use it when the data are complete (every account, every day of the period), when you must match numpy’s default var, or when a formula calls for the maximum-likelihood estimate. It is defined at \(n = 1\), where it gives 0. To infer the variance of a larger population from a sample, use variance.
spreaddata-summary
def populationVariance(x) ≔ squares(x) / count(x)
populationStd
populationStd(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values with mean \(\bar{x}\). Then
\[\operatorname{populationStd}(x) = \sigma = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^2} .\]The population standard deviation.
Intuition The root of populationVariance: the typical distance from the mean over a complete population, in the data’s units. It is smaller than std by the factor \(\sqrt{(n-1)/n}\), about 11% at five values and negligible at a few hundred.
Use it when you describe a fixed, complete group (all grades in one class, all servers in a fleet), when you must match numpy’s default std, or when you standardise within that group only. When the values are a sample and the spread is used for inference, use std.
spreaddata-summary
def populationStd(x) ≔ sqrt(populationVariance(x))
Building blocks
squares(x)— Let \(x \in \mathbb{R}^n\) with mean \(\bar{x}\). Then \(\operatorname{squares}(x) = \sum_{i=1}^{n} (x_i - \bar{x})^2\). The sum of squared deviations from the mean.
Shape
The k-th central moment is \(m_k = \frac{1}{n}\sum_i (x_i - \bar{x})^k\). The third measures asymmetry and the fourth the weight of the tails.
moment2
moment2(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values with mean \(\bar{x}\). Then
\[\operatorname{moment2}(x) = m_2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^2 .\]The second central moment (equal to the population variance).
Intuition The second central moment, equal in value to populationVariance. It exists under this name so that the moment family \(m_2, m_3, m_4\) is normalised the same way (all divided by \(n\)), which is what the shape ratios built from them require.
Use it when you compose your own standardised moments or a method-of-moments fit and want the normalisation to match moment3 and moment4. To report spread, use std or variance instead.
distribution-shapeparameter-estimation
def moment2(x) ≔ squares(x) / count(x)
moment3
moment3(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values with mean \(\bar{x}\). Then
\[\operatorname{moment3}(x) = m_3 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^3 .\]The third central moment.
Intuition Cubing keeps the sign, so large deviations above the mean add and large deviations below subtract. A long right tail makes it positive, and a symmetric sample gives about zero. It is in cubed units, so its size depends on the scale of the data.
Use it when you need raw asymmetry for a method-of-moments fit or a statistic of your own. To compare asymmetry between data sets in different units, use the scale-free skewness.
distribution-shapeparameter-estimation
def moment3(x) ≔ einsum("i,i,i->", centre(x), centre(x), centre(x)) / count(x)
moment4
moment4(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values with mean \(\bar{x}\). Then
\[\operatorname{moment4}(x) = m_4 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^4 .\]The fourth central moment.
Intuition The fourth power lets the few largest deviations dominate, so \(m_4\) measures how much of the spread comes from rare extreme values rather than many moderate ones. Its units are the data’s units to the fourth power.
Use it when you build a tail diagnostic or a method-of-moments fit yourself. For a scale-free number use kurtosis or excess. A single outlier can dominate it, so it needs many values to be stable.
distribution-shapeparameter-estimation
def moment4(x) ≔ einsum("i,i,i,i->", centre(x), centre(x), centre(x), centre(x)) / count(x)
skewness
skewness(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values, not all equal, with central moments \(m_2\) and \(m_3\). Then
\[\operatorname{skewness}(x) = g_1 = \frac{m_3}{m_2^{1.5}} .\]The moment coefficient of skewness: zero for a symmetric sample, positive when the right tail is longer. No small-sample correction is applied.
Intuition \(m_3\) divided by the cube of the population standard deviation, so it has no units: how lopsided is the distribution? Zero means symmetric. Positive means a few large values stretch the right tail (incomes, waiting times), and negative means a long left tail.
Use it when you check whether a mean or a normal assumption is reasonable, choose a transform (strong positive skew often calls for a log), or monitor a feature’s distribution over time. It is the uncorrected \(g_1\), as in scipy’s skew by default; bias-corrected versions (Excel’s SKEW, pandas’ skew) differ noticeably at small \(n\). A constant column divides by zero.
distribution-shapedata-summary
def skewness(x) ≔ moment3(x) / moment2(x)^1.5
kurtosis
kurtosis(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values, not all equal, with central moments \(m_2\) and \(m_4\). Then
\[\operatorname{kurtosis}(x) = \frac{m_4}{m_2^2} .\]The kurtosis, which is 3 for a normal distribution.
Intuition \(m_4 / m_2^2\) measures tail weight: how much of the variance comes from rare, large deviations. A normal distribution scores 3. Heavier tails, with more outliers than a normal model predicts, score above 3, and light-tailed shapes score below (a uniform distribution gives 1.8).
Use it when you judge whether extreme values are more frequent than a normal model allows, for instance before trusting a normal-based risk figure or a standard-deviation threshold. Most software reports excess kurtosis instead, so use excess when comparing numbers. No small-sample correction is applied, and it needs a large sample to be stable.
distribution-shaperisk-measurementoutlier-detection
def kurtosis(x) ≔ moment4(x) / moment2(x)^2
excess
excess(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values, not all equal, with central moments \(m_2\) and \(m_4\). Then
\[\operatorname{excess}(x) = \frac{m_4}{m_2^2} - 3 .\]The excess kurtosis, which is 0 for a normal distribution.
Intuition Kurtosis with its zero moved to the normal distribution: 0 means normal-like tails, positive means heavier tails (Student’s t, daily asset returns), negative means lighter ones. It carries the same information as kurtosis, measured from a more readable origin.
Use it when you want “how far from normal are the tails” at a glance, or you compare against libraries that report excess kurtosis. It matches scipy’s kurtosis with its defaults. pandas’ kurt applies a small-sample correction, so the two agree only for large \(n\).
distribution-shaperisk-measurement
def excess(x) ≔ kurtosis(x) − 3
error
error(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n \ge 2\) values with sample standard deviation \(s = \operatorname{std}(x)\). Then
\[\operatorname{error}(x) = \frac{s}{\sqrt{n}} .\]The standard error of the mean.
Intuition std describes how spread out the values are; the standard error describes how uncertain their mean is. Averaging independent values cancels noise, so the uncertainty shrinks with \(\sqrt{n}\): four times the data halves it.
Use it when you put error bars or a rough 95% interval (mean ± 2 standard errors, for moderate \(n\)) on a mean, compare two group means, or plan how many samples you need. It assumes independent observations; for an autocorrelated series it understates the uncertainty (check autocorrelation in Prelude.Series). To describe the spread of the values themselves, use std.
uncertainty-quantificationhypothesis-testingsample-size-planning
def error(x) ≔ std(x) / sqrt(count(x))
variation
variation(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n \ge 2\) values with mean \(\bar{x} \neq 0\) and sample standard deviation \(s\). Then
\[\operatorname{variation}(x) = \frac{s}{\bar{x}} .\]The coefficient of variation. It is meaningful for positive data only.
Intuition Spread relative to level. A standard deviation of 10 is large for a mean of 20 and tiny for a mean of 10 000. Dividing by the mean removes the units, so the variability of quantities on different scales can be compared.
Use it when you compare the relative variability of positive quantities: assay repeatability, demand volatility across products of different volume, or timing jitter. With a mean near zero it blows up, and with data that can be negative it can change sign, so do not use it there. The divisor is not guarded.
spreadquality-control
def variation(x) ≔ std(x) / mean(x)
Standardising
zscore
zscore(x) · x column n → column n
Let \(x \in \mathbb{R}^n\) be a column of \(n \ge 2\) values, not all equal, with mean \(\bar{x}\) and sample standard deviation \(s\). Then
\[\operatorname{zscore}(x)_i = \frac{x_i - \bar{x}}{s} , \qquad i = 1, \dots, n .\]The standard scores, using the sample standard deviation.
Intuition Each value re-expressed as “how many standard deviations from the mean”. The result has mean 0 and sample standard deviation 1, so columns measured in different units become directly comparable.
Use it when you scale features for distance-based or gradient-trained models (k-means, SVMs, regularised regression), flag outliers with a rule such as \(\lvert z \rvert > 3\), or add scores from different tests. The mean and std are themselves pulled by outliers, which can mask them; on contaminated data use robust. For a fixed 0-to-1 range use minmax. A constant column divides by zero.
feature-scalingoutlier-detection
def zscore(x) ≔ centre(x) / std(x)
Two columns
covariance
covariance(x, y) · x column n · y column n → scalar
Let \(x, y \in \mathbb{R}^n\) be two columns of \(n \ge 2\) paired values with means \(\bar{x}\) and \(\bar{y}\). Then
\[\operatorname{covariance}(x, y) = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y}) .\]The sample covariance.
Intuition The average product of paired deviations: positive when x and y tend to be above their means together, negative when one is above while the other is below. Its size mixes the strength of the relationship with the units of both columns, so it is hard to read on its own.
Use it when you need raw co-movement for further algebra: a regression slope (covariance(x, y) / variance(x)), the variance of a sum or a portfolio, or one entry of a covariance matrix. To judge strength use correlation; for every pair of columns at once use covariances. It is the sample form (\(n-1\)).
correlation-analysisregressionrisk-measurement
def covariance(x, y) ≔ einsum("i,i->", centre(x), centre(y)) / (count(x) − 1)
correlation
correlation(x, y) · x column n · y column n → scalar
Let \(x, y \in \mathbb{R}^n\) be two columns of \(n \ge 2\) paired values, each not constant, with sample standard deviations \(s_x\) and \(s_y\). Then
\[\operatorname{correlation}(x, y) = r = \frac{\operatorname{covariance}(x, y)}{s_x s_y} = \frac{\sum_{i} (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_i (x_i - \bar{x})^2} \sqrt{\sum_i (y_i - \bar{y})^2}} .\]Pearson’s correlation coefficient, between −1 and 1.
Intuition Covariance with the units divided out: the cosine of the angle between the two centred columns. +1 means the points lie on a rising straight line, −1 on a falling one, and 0 means no linear relationship, though a curved one may still exist.
Use it when you measure how strongly two variables move together linearly, screen features against a target, or check two predictors for collinearity. One outlier can swing it, and it misses curved relationships; for monotone or outlier-prone data use spearman or kendall. For every pair of columns use correlations. Nothing is guarded, so a constant column gives NaN.
correlation-analysisfeature-selection
def correlation(x, y) ≔ covariance(x, y) / (std(x) · std(y))
norm
norm(x) · x vector n → scalar
Let \(x \in \mathbb{R}^n\). Then
\[\operatorname{norm}(x) = \lVert x \rVert = \sqrt{\sum_{i=1}^{n} x_i^2} .\]The Euclidean norm (not centred).
Intuition The length of the vector, Pythagoras in \(n\) dimensions. It is not centred, so it measures size from zero, not spread about the mean; for a centred column it equals \(\sqrt{n-1}\) times std.
Use it when you normalise a vector to unit length (x / norm(x)), measure a distance as norm(x − y), track the size of a residual or a gradient, or need the denominator of cosine. For spread about the mean use std.
similarity-searchfeature-scaling
def norm(x) ≔ sqrt(einsum("i,i->", x, x))
cosine
cosine(x, y) · x vector n · y vector n → scalar
Let \(x, y \in \mathbb{R}^n\) be non-zero vectors. Then
\[\operatorname{cosine}(x, y) = \frac{x^{\top} y}{\lVert x \rVert \, \lVert y \rVert} = \frac{\sum_{i} x_i y_i}{\sqrt{\sum_i x_i^2} \sqrt{\sum_i y_i^2}} .\]The cosine similarity: the cosine of the angle between the two vectors.
Intuition Only direction matters: scaling either vector leaves the cosine unchanged. It is 1 for vectors pointing the same way, 0 for orthogonal ones and −1 for opposite ones. It is correlation without the centring, so an offset shared by both vectors counts as agreement.
Use it when you compare embeddings, TF-IDF document vectors or rating profiles, where a vector’s length reflects size (document length, user activity) rather than meaning. Nearest neighbours are the highest cosines. When a shared baseline should not count as similarity, use correlation. A zero vector divides by zero.
similarity-searchrecommendation
def cosine(x, y) ≔ einsum("i,i->", x, y) / (norm(x) · norm(y))
A matrix of columns
The same statistics over every column of an n×p data matrix at once. Centring subtracts the vector of column means from every row by broadcasting.
means
means(x) · x data matrix n×p → vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be a data matrix whose rows are observations. Then
The mean of each column.
Intuition mean for every column at once: one contraction collapses the row index and leaves a vector of \(p\) column centres, the “average observation” of the table.
Use it when you summarise a feature table, centre a data matrix (as centred does), or store training-set means so the same shift can be applied to new data. For one column use mean; when rows carry weights, such as cluster responsibilities, use weightedMeans.
central-tendencydata-summaryfeature-scaling
def means(x) ≔ einsum("ij->j", x) / rows(x)
centred
centred(x) · x data matrix n×p → matrix n×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be a data matrix with column means \(\bar{x}_j\). Then
The data matrix with each column’s mean subtracted.
Intuition Subtracts each column’s mean from every row by broadcasting, so every feature has mean zero and the origin sits at the centre of the data cloud. Inner products between columns then measure how they vary together instead of how large they are.
Use it when you prepare data for PCA or an SVD, form a covariance matrix by hand, or remove per-feature offsets. When new data arrive later, subtract the training set’s means, not the new data’s own. For one column use centre; to centre with row weights use weightedCentred.
feature-scalingdata-preprocessingdimensionality-reduction
def centred(x) ≔ x − means(x)
covariances
covariances(x) · x data matrix n×p → matrix p×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be a data matrix with \(n \ge 2\) rows, and \(X_c = \operatorname{centred}(X)\). Then
The sample covariance matrix.
Intuition Every column’s covariance with every other in one symmetric p×p matrix: the diagonal holds the variances and the off-diagonal entries say how pairs move together. Together they describe the shape and orientation of the data cloud as an ellipsoid.
Use it when you run PCA (its eigenvectors are the principal axes), compute Mahalanobis distances, fit a Gaussian model, or measure portfolio risk as \(w^\top C w\). It is the sample form (\(n-1\)), and it can only be invertible when there are more rows than columns. For unit-free strengths use correlations; with row weights use weightedCovariances.
correlation-analysisdimensionality-reductionrisk-measurement
def covariances(x) ≔ einsum("ij,ik->jk", centred(x), centred(x)) / (rows(x) − 1)
gram
gram(x) · x matrix n×p → matrix p×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x). Then
The Gram matrix of the columns, without centring.
Intuition All inner products between columns, with no centring and no division. It is the raw ingredient of least squares: the normal equations read \(X^\top X \beta = X^\top y\).
Use it when you solve ordinary least squares or ridge regression, especially with an intercept column of ones that must not be centred away. On already-centred columns, gram(x) divided by \(n-1\) equals covariances. It is p×p, so it stays small when there are many rows and few columns. For inner products between observations instead, use rowGram.
regressioncorrelation-analysis
def gram(x) ≔ einsum("ij,ik->jk", x, x)
variances
variances(x) · x data matrix n×p → vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be a data matrix with \(n \ge 2\) rows and column means \(\bar{x}_j\). Then
The sample variance of each column (the diagonal of covariances(x)).
Intuition The diagonal of covariances, computed directly without building the p×p matrix: each column’s sample variance, side by side.
Use it when you only need per-feature spread, for example to drop near-constant features, rank features by variability, or scale columns. The cost grows with \(np\) rather than \(np^2\), which matters when there are many columns. For spread in the data’s units use deviations; for a single column use variance. It is the sample form (\(n-1\)).
spreadfeature-selectiondata-summary
def variances(x) ≔ einsum("ij,ij->j", centred(x), centred(x)) / (rows(x) − 1)
deviations
deviations(x) · x data matrix n×p → vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be a data matrix with \(n \ge 2\) rows. Then
The sample standard deviation of each column.
Intuition std for every column: the typical distance from each column’s mean, in that column’s own units.
Use it when you standardise a matrix (centred(x) / deviations(x) is a column-wise zscore), compare how variable the features are, or store training-set scales to apply to new data. It is the sample form. A constant column gives 0, which a later division turns into NaN or infinity, so drop such columns first.
spreadfeature-scaling
def deviations(x) ≔ sqrt(variances(x))
correlations
correlations(x) · x data matrix n×p → matrix p×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be a data matrix with \(n \ge 2\) rows and no constant
column, \(C = \operatorname{covariances}(X)\) and \(s = \operatorname{deviations}(X)\). Then
The Pearson correlation matrix: each covariance divided by the product of the two columns’ standard deviations. The diagonal is 1.
Intuition The covariance matrix with every entry divided by the two columns’ standard deviations, so all features stand on the same footing: 1 on the diagonal and Pearson’s \(r\) between −1 and 1 elsewhere. It is the covariance matrix of the z-scored data.
Use it when you scan a feature table for redundant, collinear pairs, draw a correlation heatmap, or run PCA on features measured in different units. It captures only linear association; for monotone relationships use spearman pair by pair. A constant column fills its row and column with NaN.
correlation-analysisfeature-selectiondimensionality-reduction
def correlations(x) ≔ covariances(x) / spreads(x)
rowGram
rowGram(x) · x matrix n×p → matrix n×n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) have rows \(x_i^{\top}\). Then
The Gram matrix of the rows: every observation’s inner product with every other, which is what distances and cosines between observations are computed from.
Intuition The n×n matrix of inner products between observations. All pairwise geometry follows from it: squared Euclidean distances are \(G_{ii} + G_{jj} - 2G_{ij}\) and cosines are \(G_{ij}/\sqrt{G_{ii} G_{jj}}\).
Use it when you need all pairwise distances or similarities between rows, for k-nearest neighbours, clustering, or the linear kernel of an SVM or kernel ridge regression. Memory grows as \(n^2\): fine for thousands of rows, not for millions. For relations between columns use gram.
similarity-searchclustering
def rowGram(x) ≔ einsum("ik,jk->ij", x, x)
Building blocks
spreads(x)— Let \(X \in \mathbb{R}^{n \times p}\) and \(s = \operatorname{deviations}(X)\). Then \(\operatorname{spreads}(X)_{jk} = s_j s_k\), that is \(s s^{\top}\). The matrix of products of column standard deviations, used to turn covariances into correlations.
Order
Sorting is expressed without a sort. Comparing every value against every other gives an n×n matrix of ones and zeros, and the number of values below \(x_i\) is a row sum: that is the rank of \(x_i\), and every order statistic follows from it. Where several values are tied, a definition selects all of them and divides by how many there are.
rank
rank(x) · x column n → column n
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values. Then
\[\operatorname{rank}(x)_i = \frac{1}{2} \Big( \sum_{j=1}^{n} [\,x_j < x_i\,] + \sum_{j=1}^{n} [\,x_j \le x_i\,] - 1 \Big) , \qquad i = 1, \dots, n .\]The zero-based midrank of each value: the number of values below it, with tied values sharing the average of the positions they occupy. This is the ranking Spearman’s correlation is defined on.
Intuition Each value is replaced by its position in sorted order, counted from zero: distances are thrown away and only order is kept. Tied values share the average of the positions they occupy, so the ranks always add up to the same total, \(n(n-1)/2\), with or without ties.
Use it when you want a scale-free, outlier-proof version of a column: for rank correlations (spearman), non-parametric tests such as Mann–Whitney, percentile features (rank(x) / (count(x) − 1)), or leaderboard positions. It builds an n×n comparison matrix, so it costs \(O(n^2)\). For the values themselves in sorted order use sorted.
rankingrobust-statisticshypothesis-testing
def rank(x) ≔ (below(x) + atMost(x) − 1) / 2
order
order(x, k) · x column n · k zero-based position, scalar → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values and \(k \in \{0, 1, \dots, n-1\}\) a position counted from zero. Let
\[b_i = \sum_{j=1}^{n} [\,x_j < x_i\,], \qquad a_i = \sum_{j=1}^{n} [\,x_j \le x_i\,], \qquad \pi_i = [\,b_i \le k\,] \, [\,a_i > k\,] .\]Then
\[\operatorname{order}(x, k) = x_{(k)} = \frac{\sum_{i=1}^{n} x_i \, \pi_i}{\sum_{i=1}^{n} \pi_i} .\]The k-th smallest value (k = 0 is the minimum). If the position falls on several tied values, they are
averaged, which returns that shared value.
Intuition The value that sits at position k once the column is sorted, found without sorting: it is the value with at most k values strictly below it and more than k at or below it. When several tied values share that position they are averaged, which just returns their common value.
Use it when you need one specific position: the second-smallest price, the largest value but one (order(x, count(x) − 2)), or a trimmed bound. k counts from 0. Each call costs \(O(n^2)\), so for many positions call sorted once instead. For fractional positions with interpolation use quantile.
rankingdata-summary
def order(x, k) ≔ einsum("i,i->", x, picked(x, k)) / einsum("i->", picked(x, k))
least
least(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values. Then
\[\operatorname{least}(x) = x_{(0)} = \min_i x_i .\]The smallest value, as the order statistic at position 0.
Intuition The minimum, read as the order statistic at position 0: the value that nothing lies strictly below.
Use it when you need the smallest value, for example the lower edge of a range check or the best time in a set of runs. It goes through order, so it builds an n×n comparison and costs \(O(n^2)\); range and minmax reach the minimum through the built-in elementwise minimum instead. A single bad record sets it, so for a minimum that ignores the odd glitch use a low quantile such as 0.01.
data-summarydata-validation
def least(x) ≔ order(x, 0)
greatest
greatest(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values. Then
\[\operatorname{greatest}(x) = x_{(n-1)} = \max_i x_i .\]The largest value, as the order statistic at position \(n - 1\).
Intuition The maximum, read as the order statistic at the last position, \(n-1\): the value with every other value at or below it.
Use it when you need the largest value, such as a peak load, a record high, or the upper edge of a range check. Like least it costs \(O(n^2)\), where range uses the built-in elementwise maximum. The maximum grows with sample size and is set by a single record, so for reporting worst cases prefer a high quantile such as 0.99, which is far less at the mercy of one sample.
data-summarydata-validation
def greatest(x) ≔ order(x, count(x) − 1)
median
median(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values, \(h = (n-1)/2\), and \(x_{(k)} = \operatorname{order}(x, k)\). Then
\[\operatorname{median}(x) = \frac{x_{(\lfloor h \rfloor)} + x_{(\lceil h \rceil)}}{2} .\]The middle order statistic, or the mean of the two middle ones when \(n\) is even.
Intuition The middle value: half the data lie below it and half above. It depends only on order, so the largest value can grow without limit and the median does not move. It is also the constant that minimises the sum of absolute deviations.
Use it when you summarise skewed or outlier-prone data (house prices, incomes, response times), as the robust counterpart of mean, or as the centre for robust scaling. With an even count it averages the two middle values. It calls order twice, each \(O(n^2)\).
central-tendencyrobust-statisticsdata-summary
def median(x) ≔ (order(x, floor(middle(x))) + order(x, ceil(middle(x)))) / 2
quantile
quantile(x, q) · x column n · q probability in [0, 1], scalar → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values, \(0 \le q \le 1\), \(h = q(n-1)\), \(f = h - \lfloor h \rfloor\), and \(x_{(k)} = \operatorname{order}(x, k)\). Then
\[\operatorname{quantile}(x, q) = Q(q) = x_{(\lfloor h \rfloor)} (1 - f) + x_{(\lceil h \rceil)} \, f .\]The q-quantile by linear interpolation between the two order statistics around position \(q(n-1)\) (the same rule as numpy’s default).
Intuition The value below which a fraction q of the data lie. Sorted positions run from 0 to \(n-1\), so q maps to position \(q(n-1)\); when that falls between two values the answer interpolates linearly, which makes the result change smoothly with q.
Use it when you report percentiles (p95 latency, a 5% value-at-risk), set data-driven thresholds or clipping bounds, or cut a feature into equal-count bins. The rule matches numpy’s default and R’s type 7; other conventions give slightly different answers on small samples. Each call runs order twice, \(O(n^2)\) each.
data-summaryrisk-measurementrobust-statistics
def quantile(x, q) ≔ order(x, beneath(x, q)) · (1 − fraction(x, q)) + order(x, beyond(x, q)) · fraction(x, q)
q1
q1(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values. Then
\[\operatorname{q1}(x) = Q_1 = \operatorname{quantile}(x, 0.25) .\]The first quartile, quantile(x, 0.25).
Intuition The 25th percentile: a quarter of the values lie below it. With median and q3 it splits the data into four equal-count parts, and it is the bottom edge of the box in a box plot.
Use it when you build a five-number summary, describe the lower part of a skewed quantity, or compute Tukey’s lower outlier fence \(Q_1 - 1.5\,\mathrm{IQR}\). Quartile conventions differ between tools and textbooks; this one interpolates at position \(0.25(n-1)\), so hand-computed quartiles may differ on small samples.
data-summaryoutlier-detectionrobust-statistics
def q1(x) ≔ quantile(x, 0.25)
q3
q3(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values. Then
\[\operatorname{q3}(x) = Q_3 = \operatorname{quantile}(x, 0.75) .\]The third quartile, quantile(x, 0.75).
Intuition The 75th percentile: three quarters of the values lie below it, the top edge of the box in a box plot.
Use it when you describe the upper spread of a skewed quantity or compute Tukey’s upper outlier fence \(Q_3 + 1.5\,\mathrm{IQR}\). Together with q1 it defines iqr. For other tail levels, such as the 90th or 99th percentile, call quantile directly.
data-summaryoutlier-detection
def q3(x) ≔ quantile(x, 0.75)
iqr
iqr(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values. Then
\[\operatorname{iqr}(x) = Q_3 - Q_1 = \operatorname{quantile}(x, 0.75) - \operatorname{quantile}(x, 0.25) .\]The interquartile range.
Intuition The width of the middle half of the data. The top and bottom quarters do not enter, so up to a quarter of the values at either end can be arbitrarily wild without changing it. For normal data it is about 1.35 standard deviations.
Use it when you measure spread on skewed or contaminated data, flag outliers with Tukey’s fences, choose histogram bin widths (the Freedman–Diaconis rule), or scale with robust. Dividing it by 1.349 gives a robust estimate of a normal \(\sigma\). It needs four order calls, \(O(n^2)\) each.
spreadrobust-statisticsoutlier-detection
def iqr(x) ≔ q3(x) − q1(x)
range
range(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values. Then
\[\operatorname{range}(x) = \max_i x_i - \min_i x_i .\]The difference between the largest and smallest values.
Intuition The full width of the data, from the smallest value to the largest. It depends only on the two most extreme values, so one outlier sets it entirely.
Use it when you check data against physical or business limits, validate a sensor, or need the denominator of minmax. It uses the built-in elementwise maximum and minimum, so it costs \(O(n)\), unlike greatest minus least. As a measure of spread it tends to grow with sample size, so compare groups with iqr or std instead.
spreaddata-validation
def range(x) ≔ maximum(x) − minimum(x)
mode
mode(x) · x column n → scalar
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values, and let
\[t_i = \sum_{j=1}^{n} [\,x_j = x_i\,], \qquad c_i = \big[\, t_i = \max_{j} t_j \,\big] .\]Then
\[\operatorname{mode}(x) = \frac{\sum_{i=1}^{n} x_i \, c_i}{\sum_{i=1}^{n} c_i} .\]The most frequent value. When several values are equally frequent, their average is returned, since picking one of them would be arbitrary. The comparison is exact equality, so the mode is meaningful for discrete or rounded data.
Intuition The value that occurs most often, found by counting, for each value, how many entries equal it exactly. When several values tie for most frequent, their average is returned, which need not be one of the data values.
Use it when you summarise categorical codes, integer counts or rounded measurements: the most common basket size, rating or shoe size. On continuous data every value is unique, every count is 1, and the result is just the mean, so round or bin first. A sample with two equally common peaks returns the average of the two. It costs \(O(n^2)\).
central-tendencycategorical-data
def mode(x) ≔ einsum("i,i->", x, commonest(x)) / einsum("i->", commonest(x))
Building blocks
ones(x)— Let \(x\) be a tensor of any shape. Then \(\operatorname{ones}(x)\) has the shape of \(x\) and every entry \(x_i^0 = 1\). A tensor of ones shaped likex(writtenx^0, which is 1 even at 0 and at NaN).held(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{held}(x)_{ij} = x_i\) for \(i, j = 1, \dots, n\), that is \(x \, \mathbf{1}_n^{\top}\). The column spread into an n×n matrix with the value of row i everywhere in row i.laid(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{laid}(x)_{ij} = x_j\) for \(i, j = 1, \dots, n\), that is \(\mathbf{1}_n x^{\top}\). The column spread into an n×n matrix with the value of column j everywhere in column j.below(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{below}(x)_i = \sum_{j=1}^{n} [\,x_j < x_i\,]\). For each value, how many values are strictly smaller.atMost(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{atMost}(x)_i = \sum_{j=1}^{n} [\,x_j \le x_i\,]\). For each value, how many values are smaller or equal (itself included).picked(x, k)— Let \(x \in \mathbb{R}^n\) and \(k\) a position counted from zero. Then \(\operatorname{picked}(x, k)_i = [\,\operatorname{below}(x)_i \le k\,] \cdot [\,\operatorname{atMost}(x)_i > k\,]\). A mask of the values that occupy zero-based positionk.middle(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{middle}(x) = (n - 1)/2\). The halfway position.position(x, q)— Let \(x \in \mathbb{R}^n\) and \(q \in [0, 1]\). Then \(\operatorname{position}(x, q) = h = q(n - 1)\). The fractional position of a quantile.beneath(x, q)— Let \(x \in \mathbb{R}^n\) and \(q \in [0, 1]\). Then \(\operatorname{beneath}(x, q) = \lfloor q(n - 1) \rfloor\). The position rounded down.beyond(x, q)— Let \(x \in \mathbb{R}^n\) and \(q \in [0, 1]\). Then \(\operatorname{beyond}(x, q) = \lceil q(n - 1) \rceil\). The position rounded up.fraction(x, q)— Let \(x \in \mathbb{R}^n\), \(q \in [0, 1]\) and \(h = q(n - 1)\). Then \(\operatorname{fraction}(x, q) = h - \lfloor h \rfloor\). The fractional part of the position, the interpolation weight.repeats(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{repeats}(x)_i = \sum_{j=1}^{n} [\,x_j = x_i\,]\). For each value, how many values equal it (itself included).commonest(x)— Let \(x \in \mathbb{R}^n\) and \(t = \operatorname{repeats}(x)\). Then \(\operatorname{commonest}(x)_i = [\,t_i = \max_j t_j\,]\). A mask of the values that occur most often.
Robust scaling and rank correlation
robust
robust(x) · x column n → column n
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values with \(\operatorname{iqr}(x) \neq 0\). Then
\[\operatorname{robust}(x)_i = \frac{x_i - \operatorname{median}(x)}{\operatorname{iqr}(x)} , \qquad i = 1, \dots, n .\]Robust standardisation, which outliers affect far
less than zscore.
Intuition The robust twin of zscore: centre on the median and divide by the IQR instead of the mean and standard deviation. Outliers can barely move either reference, so ordinary points keep sensible scores and outliers stand out clearly instead of inflating the scale that measures them.
Use it when you scale features with heavy tails or bad records (the same as scikit-learn’s RobustScaler with its defaults), or score anomalies. On normal data its scores are about 1.35 times smaller than z-scores. When the middle half of the values is constant the IQR is 0 and the division fails. It costs \(O(n^2)\).
feature-scalingrobust-statisticsoutlier-detection
def robust(x) ≔ (x − median(x)) / iqr(x)
minmax
minmax(x) · x column n → column n
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values, not all equal. Then
\[\operatorname{minmax}(x)_i = \frac{x_i - \min_j x_j}{\max_j x_j - \min_j x_j} , \qquad i = 1, \dots, n .\]Min–max scaling onto [0, 1].
Intuition A straight-line map that sends the smallest value to 0 and the largest to 1. Relative distances are kept, so the shape of the distribution is unchanged; only the ruler changes.
Use it when a model or a display needs bounded inputs: neural-network inputs, image intensities, a colour scale, or indicators combined into an index. One outlier squeezes every other value into a narrow band, so on dirty data apply winsorised first or use robust. New data outside the original minimum and maximum land outside 0 to 1. It costs \(O(n)\), and a constant column divides by zero.
feature-scalingdata-preprocessing
def minmax(x) ≔ (x − minimum(x)) / range(x)
winsorised
winsorised(x, p) · x column n · p tail fraction in [0, 0.5], scalar → column n
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values and let \(0 \le p \le 0.5\). Then
\[\operatorname{winsorised}(x, p)_i = \min\Big( \max\big( x_i, \operatorname{quantile}(x, p) \big), \; \operatorname{quantile}(x, 1 - p) \Big) , \qquad i = 1, \dots, n .\]The column with every value below the p-quantile raised to it and every value above the (1 − p)-quantile lowered to it.
Intuition Pull the extremes in instead of deleting them: every value beyond the p and 1 − p quantiles is set to that quantile. The number of observations stays the same and no two values swap order, but no single value can dominate a mean or a variance any more.
Use it when you compute means, regressions or correlations on data with a few wild values you do not want to drop, such as financial returns, survey incomes or sensor spikes. p of 0.01 or 0.05 is typical. It changes the data, so report the p you used. It makes two quantile calls, \(O(n^2)\).
robust-statisticsdata-cleaning
def winsorised(x, p) ≔ min(max(x, quantile(x, p)), quantile(x, 1 − p))
spearman
spearman(x, y) · x column n · y column n → scalar
Let \(x, y \in \mathbb{R}^n\) be two columns of \(n \ge 2\) paired values, each not constant. Then
\[\operatorname{spearman}(x, y) = \rho = \operatorname{correlation}\big( \operatorname{rank}(x), \operatorname{rank}(y) \big) .\]Spearman’s rank correlation: Pearson’s correlation of the midranks. It measures monotone rather than linear association and handles ties.
Intuition Pearson’s correlation computed on ranks, so it asks “when x goes up, does y tend to go up?” whatever the shape of the curve. Any strictly increasing relationship, linear or exponential or logarithmic, scores 1, and an outlier can do no more than take the top or bottom rank.
Use it when a relationship is monotone but not straight, the data are ordinal (ratings, rankings), or outliers would distort correlation. Ties receive midranks. It costs \(O(n^2)\) through rank. kendall measures the same kind of association on a different scale, usually smaller in magnitude.
correlation-analysisrobust-statisticsranking
def spearman(x, y) ≔ correlation(rank(x), rank(y))
Weighted statistics
Each weighted form is the unweighted contraction with the weights added as one more operand. The weights
w are non-negative and need not sum to one; they are normalised by their total. The weighted variance
and covariances divide by the total weight, with no small-sample correction.
weight
weight(w) · w weights n → scalar
Let \(w \in \mathbb{R}^n\) be a vector of non-negative weights. Then
\[\operatorname{weight}(w) = \sum_{i=1}^{n} w_i .\]The total weight.
Intuition The total weight, which plays the role of \(n\) in every weighted formula below: with all weights equal to 1 it is exactly the count.
Use it when you normalise weights into proportions (w / weight(w)), write a weighted statistic of your own, or check that the weights are not all zero before calling weightedMean and its relatives, whose divisions are not guarded.
weighted-datadata-summary
def weight(w) ≔ einsum("i->", w)
weightedMean
weightedMean(w, x) · w weights n · x column n → scalar
Let \(w \in \mathbb{R}^n\) be non-negative weights with \(\sum_i w_i > 0\), and \(x \in \mathbb{R}^n\) a column of \(n\) values. Then
\[\operatorname{weightedMean}(w, x) = \bar{x}_w = \frac{\sum_{i=1}^{n} w_i x_i}{\sum_{i=1}^{n} w_i} .\]The weighted mean.
Intuition Each value counts in proportion to its weight, so a weight of 2 counts a value twice. Only relative weights matter: multiplying every weight by the same constant changes nothing.
Use it when rows stand for different amounts (survey weights, a frequency table whose w holds counts, position sizes in a portfolio, group means combined by group size), or when more reliable observations should count more, as with inverse-variance weights. With equal weights it is mean. Weights must be non-negative with a positive total.
central-tendencyweighted-data
def weightedMean(w, x) ≔ einsum("i,i->", w, x) / weight(w)
weightedVariance
weightedVariance(w, x) · w weights n · x column n → scalar
Let \(w \in \mathbb{R}^n\) be non-negative weights with \(\sum_i w_i > 0\), and \(x \in \mathbb{R}^n\) a column with weighted mean \(\bar{x}_w\). Then
\[\operatorname{weightedVariance}(w, x) = \frac{\sum_{i=1}^{n} w_i (x_i - \bar{x}_w)^2}{\sum_{i=1}^{n} w_i} .\]The weighted variance.
Intuition The weighted average squared distance from the weighted mean. Integer weights behave like repeating each value that many times, and the result then equals the population variance of the expanded data, because there is no \(n-1\) correction.
Use it when you measure spread in a frequency table, an importance-weighted sample or one component of a mixture. It is the population (maximum-likelihood) form: with frequency weights and a small total \(W\), multiply by \(W/(W-1)\) to get the sample variance of the expanded data. Scaling all weights leaves it unchanged. With equal weights it equals populationVariance, not variance.
spreadweighted-data
def weightedVariance(w, x) ≔ einsum("i,i,i->", w, x − weightedMean(w, x), x − weightedMean(w, x)) / weight(w)
weightedMeans
weightedMeans(w, x) · w weights n · x data matrix n×p → vector p
Let \(w \in \mathbb{R}^n\) be non-negative weights with \(\sum_i w_i > 0\), and \(X \in \mathbb{R}^{n \times p}\) (the argument x) a data matrix. Then
The weighted mean of each column, one weight per row.
Intuition weightedMean for every column at once, with one weight per row: the centre of the data cloud when some observations count more than others.
Use it when you compute centroids with soft memberships, as in the M-step of a Gaussian mixture or fuzzy k-means (where w is one cluster’s column of responsibilities), survey-weighted feature means, or the centre to subtract before weighted PCA. With equal weights it is means.
central-tendencyweighted-dataclustering
def weightedMeans(w, x) ≔ einsum("i,ij->j", w, x) / weight(w)
weightedCentred
weightedCentred(w, x) · w weights n · x data matrix n×p → matrix n×p
Let \(w \in \mathbb{R}^n\) be non-negative weights with \(\sum_i w_i > 0\), and \(X \in \mathbb{R}^{n \times p}\) (the argument x) a data matrix with weighted column means \(\bar{x}_{w,j}\). Then
The data matrix with each column’s weighted mean subtracted.
Intuition Subtracts the weighted column means, so that the weighted average of every column becomes zero. The plain, unweighted average of the result is generally not zero.
Use it when you prepare a data matrix for a weighted covariance, weighted least squares or weighted PCA, where the centring must use the same weights as the step that follows. With equal weights it is centred.
data-preprocessingweighted-data
def weightedCentred(w, x) ≔ x − weightedMeans(w, x)
weightedCovariances
weightedCovariances(w, x) · w weights n · x data matrix n×p → matrix p×p
Let \(w \in \mathbb{R}^n\) be non-negative weights with \(\sum_i w_i > 0\), and \(X \in \mathbb{R}^{n \times p}\) (the argument x) a data matrix with rows \(x_i^{\top}\) and weighted column means \(\bar{x}_w \in \mathbb{R}^p\). Then
The weighted covariance matrix.
Intuition How every pair of columns moves together, with each row’s contribution scaled by its weight. It is the covariance of a distribution that puts probability \(w_i / \sum_j w_j\) on row \(i\), which is why it divides by the total weight with no \(n-1\) correction.
Use it when you estimate a component covariance in a Gaussian mixture’s M-step, build a risk model that weights recent returns more (exponential decay weights), or run survey-weighted PCA. With equal weights it equals \((n-1)/n\) times covariances.
correlation-analysisweighted-datarisk-measurement
def weightedCovariances(w, x) ≔ einsum("i,ij,ik->jk", w, weightedCentred(w, x), weightedCentred(w, x)) / weight(w)
Sorting
order(x, k) answers one position; sorted answers every position at once, by the same selection with
the position as one more index.
sorted
sorted(x) · x column n → column n
Let \(x \in \mathbb{R}^n\) be a column of \(n\) values, with \(b_i = \sum_j [\,x_j < x_i\,]\) and \(a_i = \sum_j [\,x_j \le x_i\,]\), and let
\[S_{ki} = [\,b_i \le k - 1\,] \, [\,a_i > k - 1\,] , \qquad k, i = 1, \dots, n .\]Then
\[\operatorname{sorted}(x)_k = \frac{\sum_{i=1}^{n} S_{ki} \, x_i}{\sum_{i=1}^{n} S_{ki}} = x_{(k-1)} , \qquad k = 1, \dots, n ,\]where \(x_{(k-1)} = \operatorname{order}(x, k - 1)\) is the order statistic at position \(k - 1\) counted from zero.
The values in ascending order. Tied values each appear as many times as they occur.
Intuition Every order statistic at once. Row \(k\) of an n×n selection mask marks the value that belongs at sorted position \(k\), using the same below and at-most counts as rank, so no comparison sort ever runs. Tied values appear as many times as they occur.
Use it when you need many positions together, such as an empirical CDF, a Q–Q plot, a trimmed mean or a top-k readout, rather than calling order repeatedly. Time and memory are \(O(n^2)\), so keep \(n\) in the thousands. For each original entry’s position in the order, use rank instead.
rankingdata-summary
def sorted(x) ≔ einsum("ki,i->k", selects(x), x) / einsum("ki->k", selects(x))
Building blocks
positions(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{positions}(x)_j = \sum_{i=1}^{n} [\,j \ge i\,] - 1 = j - 1\) for \(j = 1, \dots, n\). The positions \(0, 1, \dots, n-1\) as a vector, built from a triangle of ones.sought(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{sought}(x)_{ki} = k - 1\) for \(k, i = 1, \dots, n\). The positions spread into an n×n matrix, one row per position.fewer(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{fewer}(x)_{ki} = \operatorname{below}(x)_i = \sum_j [\,x_j < x_i\,]\). The strictly-below counts spread into an n×n matrix, one column per value.upto(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{upto}(x)_{ki} = \operatorname{atMost}(x)_i = \sum_j [\,x_j \le x_i\,]\). The at-most counts spread into an n×n matrix, one column per value.selects(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{selects}(x)_{ki} = [\,\operatorname{below}(x)_i \le k - 1\,] \cdot [\,\operatorname{atMost}(x)_i > k - 1\,]\), so row \(k\) equals \(\operatorname{picked}(x, k - 1)\). An n×n mask: row k marks the values that occupy position k − 1 counted from zero.
Kendall’s tau
For every ordered pair of observations, the product of the signs of the differences is +1 on a concordant pair, −1 on a discordant pair and 0 where either variable is tied, so one sum counts concordant minus discordant pairs.
kendall
kendall(x, y) · x column n · y column n → scalar
Let \(x, y \in \mathbb{R}^n\) be two columns of \(n \ge 2\) paired values, each not constant. Let
\[T_x = \sum_{i=1}^{n} \sum_{j=1}^{n} [\,x_i = x_j\,] - n , \qquad T_y = \sum_{i=1}^{n} \sum_{j=1}^{n} [\,y_i = y_j\,] - n .\]Then
\[\operatorname{kendall}(x, y) = \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} \operatorname{sign}(x_j - x_i) \operatorname{sign}(y_j - y_i)}{\sqrt{\big(n(n-1) - T_x\big)\big(n(n-1) - T_y\big)}} = \frac{n_c - n_d}{\sqrt{(n_0 - t_x)(n_0 - t_y)}} = \tau_b ,\]where \(n_c\) and \(n_d\) count concordant and discordant pairs, \(n_0 = n(n-1)/2\), and \(t_x = T_x / 2\), \(t_y = T_y / 2\) count the unordered pairs tied in each variable. The sums run over ordered pairs, so numerator and denominator both carry a factor of two, which cancels.
Kendall’s rank correlation \(\tau_b\), which corrects the denominator for pairs tied in each variable.
Intuition Look at every pair of observations: if x and y put the pair in the same order it is concordant, if in opposite orders discordant. \(\tau\) is concordant minus discordant over the comparable pairs, so without ties \(\tau = 0.5\) means a random pair is three times as likely to agree as to disagree. The b-variant removes tied pairs from the denominator.
Use it when you compare two rankings (search results, judges’ scores), have ordinal data with many ties, or a small sample, where its sampling behaviour is better than spearman’s. It is usually smaller in magnitude than Spearman’s on the same data. It compares all pairs, \(O(n^2)\), and a constant column divides by zero.
correlation-analysisrankingrobust-statistics
def kendall(x, y) ≔ concord(x, y) / sqrt((pairs(x) − bound(x)) · (pairs(y) − bound(y)))
Building blocks
concord(x, y)— Let \(x, y \in \mathbb{R}^n\). Then \(\operatorname{concord}(x, y) = \sum_{i=1}^{n} \sum_{j=1}^{n} \operatorname{sign}(x_j - x_i) \operatorname{sign}(y_j - y_i) = 2(n_c - n_d)\). Twice the number of concordant pairs minus discordant pairs, summed over ordered pairs.pairs(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{pairs}(x) = n(n - 1)\). The number of ordered pairs.bound(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{bound}(x) = \sum_{i=1}^{n} \sum_{j=1}^{n} [\,x_i = x_j\,] - n\). The number of ordered pairs \((i, j)\) with \(i \neq j\) tied inx.
Examples
Summarise a column
open Prelude.Statistics
fn averageOf(x) ↦ mean(x)
fn sampleVariance(x) ↦ variance(x)
fn sampleStd(x) ↦ std(x)
fn popVariance(x) ↦ populationVariance(x)
fn skew(x) ↦ skewness(x)
using Examples;
using MLambda.Hilbert.Runtime;
var x = new[] { 2d, 4d, 4d, 4d, 5d, 5d, 7d, 9d };
var mean = Stats.AverageOf((Tensor)x);
var variance = Stats.SampleVariance((Tensor)x);
var std = Stats.SampleStd((Tensor)x);
var population = Stats.PopVariance((Tensor)x);
var skew = Stats.Skew((Tensor)x);
Output (verified, rounded to 6 significant digits):
mean = 5
variance = 4.57143
std = 2.13809
population = 4
skew = 0.65625
The squared deviations sum to 32, so the sample variance is 32/7 and the population variance 32/8.
Median, quartiles and mode
open Prelude.Statistics
fn middleOf(x) ↦ median(x)
fn quantileOf(x, q) ↦ quantile(x, q)
fn spreadOf(x) ↦ iqr(x)
fn commonestValue(x) ↦ mode(x)
using Examples;
using MLambda.Hilbert.Runtime;
var x = new[] { 2d, 4d, 4d, 4d, 5d, 5d, 7d, 9d };
var median = Stats.MiddleOf((Tensor)x);
var lower = Stats.QuantileOf((Tensor)x, (Tensor)0.25);
var iqr = Stats.SpreadOf((Tensor)x);
var mode = Stats.CommonestValue((Tensor)x);
Output (verified):
median = 4.5
lower = 4
iqr = 1.5
mode = 4
With eight values the median averages the 4th and 5th (4 and 5). The third quartile sits at position \(0.75 \cdot 7 = 5.25\), a quarter of the way from 5 to 7, which is 5.5.
Ranks and sorting, with ties
open Prelude.Statistics
fn ranks(x) ↦ rank(x)
fn ascending(x) ↦ sorted(x)
fn kth(x, k) ↦ order(x, k)
using Examples;
using MLambda.Hilbert.Runtime;
var r = Stats.Ranks((Tensor)new[] { 2d, 4d, 4d, 4d, 5d, 5d, 7d, 9d });
var s = Stats.Ascending((Tensor)new[] { 3d, 1d, 2d, 1d });
var third = Stats.Kth((Tensor)new[] { 5d, 1d, 4d, 2d }, (Tensor)2.0);
Output (verified):
r = [0, 2, 2, 2, 4.5, 4.5, 6, 7]
s = [1, 1, 2, 3]
third = 4
The three 4s occupy positions 1, 2 and 3 and share the midrank 2. Position 2 (zero-based) of
{ 5, 1, 4, 2 } is the third smallest value, 4.
Covariance and correlation matrices
open Prelude.Statistics
fn columnMeans(x) ↦ means(x)
fn covarianceMatrix(x) ↦ covariances(x)
fn correlationMatrix(x) ↦ correlations(x)
using Examples;
using MLambda.Hilbert.Runtime;
Tensor data = new double[,] { { 1, 2 }, { 2, 4 }, { 3, 5 }, { 4, 9 } };
var means = Stats.ColumnMeans(data);
var cov = Stats.CovarianceMatrix(data);
var cor = Stats.CorrelationMatrix(data);
Output (verified, rounded to 6 significant digits):
means = [2.5, 5]
cov = [[1.66667, 3.66667], [3.66667, 8.66667]]
cor = [[1, 0.964764], [0.964764, 1]]
Three correlation coefficients
open Prelude.Statistics
fn pearson(x, y) ↦ correlation(x, y)
fn rankCorrelation(x, y) ↦ spearman(x, y)
fn tau(x, y) ↦ kendall(x, y)
using Examples;
using MLambda.Hilbert.Runtime;
var x = new[] { 2d, 4d, 4d, 4d, 5d, 5d, 7d, 9d };
var y = new[] { 1d, 3d, 2d, 5d, 4d, 6d, 8d, 7d };
var r = Stats.Pearson((Tensor)x, (Tensor)y);
var rho = Stats.RankCorrelation((Tensor)x, (Tensor)y);
var t = Stats.Tau((Tensor)x, (Tensor)y);
var simple = Stats.Tau((Tensor)new[] { 1d, 2d, 3d, 4d }, (Tensor)new[] { 1d, 3d, 2d, 4d });
Output (verified, rounded to 6 significant digits):
r = 0.845594
rho = 0.883852
t = 0.771517
simple = 0.666667
In the last case one of the six pairs (2 and 3) is discordant, so \(\tau = (5 - 1)/6\).
Weighted mean and variance
open Prelude.Statistics
fn weightedAverage(w, x) ↦ weightedMean(w, x)
fn weightedSpread(w, x) ↦ weightedVariance(w, x)
using Examples;
using MLambda.Hilbert.Runtime;
var w = new[] { 1d, 1d, 2d };
var x = new[] { 1d, 2d, 4d };
var mean = Stats.WeightedAverage((Tensor)w, (Tensor)x);
var variance = Stats.WeightedSpread((Tensor)w, (Tensor)x);
Output (verified):
mean = 2.75
variance = 1.6875
The weight of 2 on the value 4 counts it twice: the result equals the population mean and variance of
{ 1, 2, 4, 4 }.
Tame outliers
open Prelude.Statistics
fn clipped(x, p) ↦ winsorised(x, p)
fn scaledRange(x) ↦ minmax(x)
fn standardised(x) ↦ zscore(x)
using Examples;
using MLambda.Hilbert.Runtime;
var w = Stats.Clipped((Tensor)new[] { 1d, 2d, 3d, 4d, 100d }, (Tensor)0.25);
var m = Stats.ScaledRange((Tensor)new[] { 2d, 4d, 6d });
var z = Stats.Standardised((Tensor)new[] { 1d, 2d, 3d });
Output (verified):
w = [2, 2, 3, 4, 4]
m = [0, 0.5, 1]
z = [-1, 0, 1]
The 0.25- and 0.75-quantiles of the first column are 2 and 4, so 1 is raised to 2 and the outlier 100 is lowered to 4.