Descriptive Statistics and EDA

You are handed a users table with a lifetime_value column and asked for one number that "describes the typical customer." You compute the mean: $1,284. You compute the median: $92. Both are correct, and quoting the wrong one will make your product roadmap either celebratory or panicked. This is the first question of every EDA session and the most common way summaries mislead: center and spread are separate claims, and skewed data splits them apart.

The mean is the arithmetic average, (1/n)·Σxᵢ. Every value pulls it — one $900k whale moves the mean of 100 users by $9,000. The median is the 50th percentile: sort the values, take the middle. It moves only when values cross the middle rank, so it resists outliers by construction. That resistance is the whole reason both numbers exist:

For ML, this choice is not cosmetic. Linear regression minimizes squared error, so it is a mean-seeking model: a log transform on a heavy-tailed target literally asks the model to fit the geometric mean instead of the arithmetic mean. When your mean and median disagree by more than ~20–30% of the median, treat the feature as skewed and plan a transform or a tree-based model. You will feel this again in Key Distributions for ML, where the mean/variance pair uniquely parameterizes the normal but tells you far less about skewed families.

Work the numbers by hand once so the formulas stop being abstract. Five customers paid [40, 42, 45, 50, 500]: mean = (40+42+45+50+500)/5 = 677/5 = 135.4; median = 45 (the middle of the sorted list). One paying customer moved the mean by $90 over the median. Trimmed means exist for this: drop the top and bottom k% and average the rest — scipy.stats.trim_mean(x, 0.2) would return 46.5 here.

import numpy as np
from scipy.stats import trim_mean

paid = np.array([40, 42, 45, 50, 500])
print(paid.mean(), np.median(paid))        # 135.4  45.0
print(trim_mean(paid, 0.2))                # 46.5 — one whale removed
print(paid.mean() < 0.8 * np.median(paid) or
      paid.mean() > 1.2 * np.median(paid)) # True → flag as skewed

An interview-ready framing: "I report the median for skewed targets, the mean only after checking the distribution, and I log-transform targets where the ratio of mean to median exceeds roughly 1.5." That single sentence signals you have looked at real data, not just textbook columns.

Record the population and sampling window beside every summary. A median from one region or season may not describe a different deployment cohort, even when its calculation is correct.

Compare summaries across meaningful groups and include sample counts. A median based on a handful of observations is less stable than one supported by a broad cohort.

Two recommendation models are evaluated on the same test set. Model A has mean absolute error $3.10 with a standard deviation of $0.40. Model B has the same $3.10 mean with a standard deviation of $4.80. Same headline, different products: Model A is boringly consistent; Model B occasionally prices an item $20 off. Ship decisions ride on spread, and the two numbers that describe it are variance and standard deviation.

Variance is the average squared deviation from the mean: σ² = (1/n)·Σ(xᵢ − μ)². The population form divides by n; the sample estimate s² divides by n − 1 (Bessel's correction) because dividing by n systematically underestimates the population variance when you estimate it from a sample — the deviations from a sample mean are smaller than the deviations from the true mean. NumPy's default np.var(x) divides by n; use np.var(x, ddof=1) for the unbiased estimate, which is what most statistical tests expect.

The squared units are variance's interpretability problem: if your feature is in dollars, the variance is in dollars squared. Standard deviation fixes that by taking the square root: σ = √(variance), back in the original units. Variance is the algebra-friendly quantity (it adds cleanly across independent sources); standard deviation is the reportable one. Both are maximized by distant outliers — that is their weakness.

What makes standard deviation intuitive is its behavior under the normal distribution: the empirical rule puts ~68% of the data within ±1σ of the mean, ~95% within ±2σ, ~99.7% within ±3σ. So an error of $15.50 with σ = $0.40 is a 38σ event — impossible under normality, which is exactly how you detect that the error distribution has a tail your mean is hiding. For non-normal data, Chebyshev's inequality still guarantees at least 1 − 1/k² of the data lies within k standard deviations (≥75% within 2σ, ≥89% within 3σ) — weaker, but it never assumes a shape.

For skewed features, the interquartile range (IQR = Q3 − Q1) is the robust alternative: it ignores everything outside the middle 50% and therefore cannot be blown up by outliers. The rule of thumb for boxplots and outlier fences follows directly: any point beyond Q1 − 1.5·IQR or Q3 + 1.5·IQR is flagged as a potential outlier.

import numpy as np

errors = np.array([2.9, 3.0, 3.1, 3.1, 3.2, 3.4, 15.5])
var_n  = np.var(errors)            # 12.06 — biased (divides by n)
var_1  = np.var(errors, ddof=1)    # 14.07 — unbiased estimate
sd     = np.std(errors, ddof=1)    # 3.75 — original units
q1, q3 = np.percentile(errors, [25, 75])
iqr = q3 - q1
lo, hi = q1 - 1.5 * iqr, q3 + 1.5 * iqr
print(var_n, sd, lo, hi)
print(errors[(errors < lo) | (errors > hi)])  # [15.5] — the tail you'd miss by the mean

Interview framing: "Standard deviation measures typical distance from the mean and is only fully meaningful when the distribution is roughly symmetric; for skewed features I report IQR because the variance is dominated by the tail." Add that variance is additive across independent variables (σ²_total = σ²_A + σ²_B) and you have shown the algebraic reason models optimize squared error in Simple Linear Regression.

A colleague fits a linear regression on purchase_amount and reports an R² of 0.11 with residual plots that fan out to the right. Nothing is wrong with the code. The distribution of the target is heavily right-skewed (skewness ≈ 2.4) with a fat tail (excess kurtosis ≈ 11), and the model is spending its capacity chasing a handful of extreme orders. Shape statistics exist precisely to catch this before the fit, not after.

Skewness is the standardized third moment: the average cubed deviation from the mean, scaled by σ³ so it is unitless. Positive skew means a long right tail: mean > median > mode. Negative skew mirrors it. The practical read: positive skew plus squared-error loss means large residuals dominate the gradient, so a few points steer the fit. The standard remedy is a log transform — log1p(x) handles zeros and turns multiplicative tails into additive ones — after which the residuals usually tighten and the mean/median gap closes.

Kurtosis is the standardized fourth moment; most libraries report excess kurtosis (subtracting 3 so a normal distribution reads 0). Positive excess kurtosis means more mass in the tails and a sharper peak than a normal — the technical name for "extreme values happen more often than your model assumes." Excess kurtosis of 11 says tail events are common enough that a squared-error model will meet them constantly. One caveat worth saying aloud in interviews: kurtosis is extremely sensitive to a single outlier, so always pair the number with a histogram or QQ plot rather than acting on the statistic alone.

The EDA loop that catches all three problems — center, spread, shape — in ten minutes:

import pandas as pd
import numpy as np
from scipy.stats import skew, kurtosis

rng = np.random.default_rng(7)
orders = pd.DataFrame({
    "amount": np.round(np.exp(rng.normal(3.2, 0.9, 5000)), 2),
})

x = orders["amount"]
print(x.describe()[["mean", "50%", "std"]])   # mean well above median
print(skew(x), kurtosis(x))                      # ~2+, ~8+ → fat right tail

log_x = np.log1p(x)
print(skew(log_x), kurtosis(log_x))              # near 0 — tamed

q1, q3 = x.quantile([0.25, 0.75])
fences = (q1 - 1.5 * (q3 - q1), q3 + 1.5 * (q3 - q1))
print((x < fences[0]).sum(), (x > fences[1]).sum())  # count of flagged outliers

Notice the loop is a loop: measure the shape, transform, re-measure. skew(log_x) dropping toward 0 is the evidence the transform worked — not the fact that you applied one. This habit of verifying rather than assuming is exactly what Missing Values and Outliers formalizes when it decides between dropping, imputing, and flagging. And once your features look like a shape you can name, Key Distributions for ML gives you the parameter tables to model that shape directly.