#stats - statistics: descriptive, distributions, and inference

The stats library is the whole statistical toolkit: descriptive statistics over a list of int or list of float, the probability distributions those build on, and the classical inference layer (regression, confidence intervals, hypothesis tests). It is pure-value and dependency-free (Go stdlib only, the distribution CDFs reusing math's special functions), so both binaries build it. There is deliberately no separate prob library - distributions live with the descriptive stats and inference they are used with, as in scipy.stats or R's base dnorm / pnorm / qnorm.

use stats;
def xs as list of int init [2, 4, 4, 4, 5, 5, 7, 9];
stats.mean($xs);              # 5.0
stats.median($xs);            # 4.5
stats.stddev($xs);            # 2.0
stats.percentile($xs, 90);    # 7.6

#Functions

Each function accepts a list of int or a list of float (an int element promotes to float for the real-valued reductions). A non-numeric element, or a non-list argument, is a positioned error.

#Kind-preserving vs. real-valued

The reductions that are inherently real-valued - mean, median, variance, stddev, percentile, correlation - always return float. The selections that pick an element - min, max, mode - and sum preserve the input kind: over a list of int they return an int, over a list of float a float.

stats.min([2, 4, 9]);        # 2   (int)
stats.min([2.0, 4.0]);       # 2.0 (float)
stats.sum([1, 2, 3]);        # 6   (int, overflow-checked)
stats.sum([1.5, 2.5]);       # 4.0 (float)

#Population vs. sample (σ vs s)

variance / stddev are population statistics (σ², σ - divide by n) - the list is treated as the whole data set, matching NumPy's default. sampleVariance / sampleStddev are the sample (s², s - Bessel-corrected, ÷ (n-1)) forms for inferential use, and need at least two elements:

population   σ² = Σ(x - μ)² / n           σ = √σ²
sample       s² = Σ(x - x̄)² / (n - 1)     s = √s²

where μ (population) / x̄ (sample) is the mean. The sample form s is what some sources call the empirical standard deviation; Jennifer names the two forms explicitly (stddev vs sampleStddev) rather than by that ambiguous label. covariance follows the population convention (σ_xy, ÷ n) to match variance; sampleCovariance is the sample form (s_xy, ÷ (n-1)). The normalized correlation divides out the standard deviations, so its n factor cancels and the population/sample distinction does not matter there. The companion linalg library does linear algebra over the same value types; further ML primitives atop the two are on the horizon.

#Percentile

stats.percentile(xs, p) sorts the data and interpolates linearly between the two nearest ranks (NumPy's default method): p = 0 is the minimum, p = 100 the maximum, and p = 50 equals median. p may be an int or float; a value outside [0, 100], or an empty list, is an error.

#Correlation

stats.correlation(xs, ys) is the Pearson correlation coefficient in [-1, 1]: 1.0 for a perfect positive linear relationship, -1.0 for a perfect negative one, 0.0 for none. The two lists must be the same length with at least two elements, and neither may have zero variance (a constant list) - all of which are catchable errors. stats.covariance(xs, ys) is the raw, units-aware version (the numerator, ÷ n).

#Quartiles, IQR, and z-scores

stats.quartiles(xs) returns [Q1, Q2, Q3] ([p25, p50, p75]) as a list of float; stats.iqr(xs) is p75 - p25, useful for spread and outlier detection. stats.range(xs) is the full spread max - min (kind-preserving). stats.zscore(xs) standardizes every element to (x - mean) / stddev (population), returning a list of float - a constant list has no z-score and errors.

#Geometric and harmonic means

stats.geometricMean(xs) (the n-th root of the product, computed via logs so a long list does not overflow) suits growth rates and ratios; stats.harmonicMean(xs) (n / sum of reciprocals) suits averaging rates. Both require every element to be strictly positive. stats.weightedMean(xs, weights) averages xs weighted by weights (same length; the weights must not sum to zero).

#Shape

stats.skewness(xs) is the population skewness (the standardized third moment): 0 for a symmetric distribution, positive when the right tail is longer. stats.kurtosis(xs) is the excess kurtosis (the standardized fourth moment minus 3), so a normal distribution is 0 and heavier-than-normal tails are positive. Both need at least two elements and a non-constant list. stats.mad(xs) (median absolute deviation) and stats.modes(xs) (every modal value, versus mode's single one) round out the robust / multimodal summaries.

#describe

stats.describe(xs) computes a whole summary in one pass and returns a stats.Summary struct - count (int) plus min, q1, median, mean, q3, max, and stddev (population, all float):

def s as stats.Summary init stats.describe($data);
io.printf("median %v, IQR %v\n", $s.median, $s.q3 - $s.q1);

#Distributions

The probability distributions inference is built on, with flat R-style names. Every function returns a float. The continuous CDFs reduce to math's regularized incomplete gamma / beta; the quantiles (inverse CDFs) invert them numerically. Sampling draws on math's shared random source (so math.randSeed makes it reproducible).

A quantile probability must be in the open interval (0, 1) (the boundaries map to infinity); an out-of-domain parameter is a catchable error, never a NaN.

use stats;
stats.normalCdf(1.96, 0, 1);        # 0.975...
stats.normalQuantile(0.975, 0, 1);  # 1.959...
stats.tQuantile(0.975, 10);         # 2.228... (t critical value)
stats.binomialPmf(5, 10, 0.5);      # 0.2461

#Inference

Regression, confidence intervals, and hypothesis tests. Results come back as small structs whose fields are:

def struct Regression { n as int, slope as float, intercept as float, r as float, r2 as float, stdErr as float, pValue as float };
def struct Interval { lower as float, upper as float };
def struct Test { statistic as float, df1 as float, df2 as float, pValue as float };

Regression is a simple-OLS fit: slope / intercept, the correlation r and its square r2, the slope's standard error stdErr, its two-sided pValue, and the point count n. Interval is a [lower, upper] confidence interval. Test carries the test statistic, its degrees of freedom df1 / df2 (df2 is 0 for a single-df test like the t or chi-square; both are set for an F / ANOVA), and the pValue.

A degenerate input (zero variance, a singular design, a non-positive expected count, a negative observed count, or magnitudes that overflow the computation) is a catchable error.

A few boundary and numerical notes:

• Proportion intervals at the extremes. For proportionCi at successes == 0 or successes == n, the "wald" and "wilson" formulas collapse to a degenerate point interval ([0, 0] / [1, 1]) - the well-known boundary defect of those methods. Use "clopper-pearson" (the exact interval) when a boundary count matters; it yields a proper one-sided bound there.
• Regression conditioning. linearRegression / multipleRegression solve the normal equations by Gaussian elimination (the same approach as linalg.solve). That is fast and exact for well-conditioned data, but numerically weak for near-collinear predictors or extreme scales; the coefficients can lose precision before the design is flagged singular. Center / scale the predictors when they span very different magnitudes.
• Cost. multipleRegression is O(n*k²) to form the system and O(k³) to solve it, for k predictors - fine for the modest tabular data a tree-walker handles, not a many-thousand-column design.

use stats;
use io;

def r as stats.Regression init stats.linearRegression([1, 2, 3, 4, 5], [2.1, 3.9, 6.2, 7.8, 10.1]);
io.printf("y = %.2fx + %.2f (r2 %.3f, p %.4f)\n", $r.slope, $r.intercept, $r.r2, $r.pValue);

def ci as stats.Interval init stats.proportionCi(8, 10, 0.95, "clopper-pearson");
io.printf("95%% CI: [%.3f, %.3f]\n", $ci.lower, $ci.upper);

def t as stats.Test init stats.tTest([5.1, 4.9, 5.2, 4.8, 5.0], 4.5);
io.printf("t = %.3f, p = %.4f\n", $t.statistic, $t.pValue);

#Strictness

Following math's stance, a non-finite or mathematically undefined result is a catchable error rather than a NaN or Inf: an empty list (for every reduction except sum), a percentile p outside [0, 100], mismatched or too-short correlation inputs, a zero-variance correlation, and an int sum that overflows int64 all raise a positioned error you can try / catch. The same applies when an input's magnitudes overflow the computation - a list of float with values near the float64 ceiling (~1.8e308) can push an intermediate sum or power to +/-Inf (and thence NaN); every real-valued reduction rejects a non-finite result instead of returning it, so a NaN never escapes into the type system.