9. Basic Statistics

Level: Beginner You'll use this for: dashboards, SLOs, experiments, and "is this normal" judgment calls.

Mean, median, mode

• Mean: sum ÷ count. Sensitive to outliers.
• Median: the middle value when sorted. Robust to outliers.
• Mode: the most frequent value.

Why this matters for latency dashboards: if your API has 99 requests at 50ms and 1 request at 50,000ms (a stuck connection), the mean jumps to ~549ms — misleadingly bad. The median stays at 50ms — correctly reflecting "almost everyone had a fine experience." This is exactly why SRE dashboards report percentiles, not averages.

Percentiles (p50, p95, p99)

The value below which X% of observations fall. p99 = 500ms means 99% of requests were faster than 500ms (and 1% were slower — that 1% is often where your angriest users live).

import numpy as np
latencies = [45, 48, 50, 52, 51, 49, 5000, 47, 53, 46]  # ms, one outlier
print("mean:", np.mean(latencies))     # 544.1  <- misleading
print("median:", np.median(latencies)) # 49.0   <- representative
print("p95:", np.percentile(latencies, 95))
print("p99:", np.percentile(latencies, 99))

Variance and standard deviation

Variance measures spread: average of squared distances from the mean. Standard deviation (σ, the square root of variance) puts that spread back into the original units, which is why alerting thresholds are usually phrased as "mean ± k·σ" ("3-sigma alert").

variance = (1/n) · Σ(xᵢ - mean)²
std_dev  = sqrt(variance)

data = [10, 12, 9, 11, 60]   # one clear outlier
print(np.mean(data), np.std(data))

Computing this by hand, step by step:

data = [10, 12, 9, 11, 60]

1. Mean:  (10+12+9+11+60) / 5  =  102 / 5  =  20.4

2. Deviation from mean, each point (xi - mean):
     10 - 20.4 = -10.4
     12 - 20.4 =  -8.4
      9 - 20.4 = -11.4
     11 - 20.4 =  -9.4
     60 - 20.4 =  39.6

3. Square each deviation (removes negative signs, penalizes big gaps more):
     (-10.4)² = 108.16
     (-8.4)²  =  70.56
     (-11.4)² = 129.96
     (-9.4)²  =  88.36
     (39.6)²  = 1568.16

4. Average the squared deviations (this average IS the variance):
     (108.16+70.56+129.96+88.36+1568.16) / 5 = 1965.2 / 5 = 393.04

5. Standard deviation = sqrt(variance):
     sqrt(393.04) ≈ 19.82

Notice how step 3 makes the outlier (60) dominate the whole result — its squared deviation (1568.16) is more than 10× any other point's, which is exactly why one wild data point drags both variance and standard deviation up so much.

The normal (Gaussian) distribution and the 68-95-99.7 rule

For normally distributed data: ~68% of values fall within 1σ of the mean, ~95% within 2σ, ~99.7% within 3σ. This is why "3-sigma" is a common default anomaly threshold — a point beyond 3σ has < 0.3% chance of occurring "naturally," making it a reasonable (if crude) alert trigger.

Sample vs. population: why some formulas divide by n - 1

The variance above divides by n. That is correct when your data is the entire population. Usually it's only a sample (the last 1,000 requests, not every request ever), and then you should divide by n - 1:

population variance:  σ² = (1/N)     · Σᵢ (xᵢ - μ)²
sample variance:      s² = (1/(n-1)) · Σᵢ (xᵢ - x̄)²

Why n - 1 (Bessel's correction)? The deviations are measured from x̄, which was computed from the same data. That makes the data look closer to its own centre than it is to the true μ, so dividing by n systematically underestimates the spread. Another way to see it: once you know x̄, the deviations must sum to zero, so only n - 1 of them are free.

Example: the same data [10, 12, 9, 11, 60] gives 1965.2 / 5 = 393.04 with n, but 1965.2 / 4 = 491.3 (s ≈ 22.17) with n - 1. NumPy's np.std divides by n by default (ddof=0) and pandas' .std() divides by n - 1 (ddof=1), so the "same" calculation can disagree between two notebooks. Pass ddof explicitly.

z-scores: putting anything on a common ruler

z = (x - μ) / σ

Example: a host normally sends μ = 200 MB/hour of egress traffic with σ = 40 MB. This hour it sent 380 MB:

z = (380 - 200) / 40 = 4.5

Under a normal model, a 4.5σ event happens about 3 times in a million hours. That's worth paging for, since it might be data exfiltration. Dividing by σ is also why ML pipelines standardize features before training. Packet sizes (hundreds) and port numbers (thousands) end up on the same scale, so gradient descent (Ch. 21) doesn't spend all its effort on whichever feature has the biggest raw numbers.

Covariance and correlation: do two things move together?

Cov(X, Y) = (1/(n-1)) · Σᵢ (xᵢ - x̄)(yᵢ - ȳ)

r = Cov(X, Y) / (s_X · s_Y)                 (Pearson correlation, always in [-1, 1])

Example (worked by hand): five load-test runs, with x = payload size (KB) and y = p50 latency (ms):

x = [1, 2, 3, 4, 5]     x̄ = 3
y = [2, 4, 5, 4, 5]     ȳ = 4

xᵢ - x̄:            -2  -1   0   1   2
yᵢ - ȳ:            -2   0   1   0   1
product:             4   0   0   0   2     sum = 6

Cov = 6 / 4 = 1.5
s_X² = (4+1+0+1+4)/4 = 2.5      s_Y² = (4+0+1+0+1)/4 = 1.5
r = 1.5 / sqrt(2.5 · 1.5) = 1.5 / 1.936 ≈ 0.77

That's a strong positive relationship. The best-fit line's slope is Cov/s_X² = 1.5/2.5 = 0.6 ms per KB, and the intercept is ȳ - 0.6·x̄ = 2.2 ms. That's linear regression in two lines. Ch. 10 derives the same answer from matrices, which scales to thousands of features.

import numpy as np
x = np.array([1, 2, 3, 4, 5]); y = np.array([2, 4, 5, 4, 5])
print(np.cov(x, y)[0, 1], np.corrcoef(x, y)[0, 1])   # 1.5  0.7746
print(np.polyfit(x, y, 1))                           # [0.6 2.2]  slope, intercept

Remember: r = 0 does not mean "unrelated". y = x² on x ∈ [-1, 1] has r = 0 but is perfectly determined by x. Correlation only measures straight-line relationships. Mutual information (Ch. 23) catches the rest.

Exponentially weighted moving average (EWMA): statistics with a memory budget

A plain average needs every past sample. An EWMA keeps a single number and nudges it toward each new observation:

S_new = (1 - α) · S_old + α · x_new

Unroll it and you'll see sample k steps old gets weight α(1-α)ᵏ, which is geometric decay (Ch. 6). Its "memory" is roughly 1/α samples.

Networking in the wild: TCP's retransmission timer (RFC 6298). Every TCP connection in your kernel runs two EWMAs, one on round-trip time and one on its variability, to decide how long to wait before retransmitting:

RTTVAR ← (1 - β) · RTTVAR + β · |SRTT - R|       β = 1/4
SRTT   ← (1 - α) · SRTT   + α · R                α = 1/8
RTO    =  SRTT + 4 · RTTVAR

Worked by hand. The first sample R = 100 ms initializes SRTT = 100, RTTVAR = R/2 = 50, so RTO = 100 + 4·50 = 300 ms. Then:

R = 120 ms:  RTTVAR = 0.75·50   + 0.25·|100   - 120| = 37.5   + 5      = 42.5
             SRTT   = 0.875·100 + 0.125·120          = 102.5
             RTO    = 102.5 + 4·42.5                 = 272.5 ms

R = 300 ms:  RTTVAR = 0.75·42.5  + 0.25·|102.5 - 300| = 31.875 + 49.375 = 81.25
             SRTT   = 0.875·102.5 + 0.125·300         = 127.19
             RTO    = 127.19 + 4·81.25                = 452.19 ms

One slow sample moves the mean by only 25 ms, but the timeout by 180 ms, because the variance term reacts to the surprise. (Real stacks also clamp RTO to a minimum: 1 s in the RFC, 200 ms on Linux.)

AI in the wild: the Adam optimizer (Ch. 21) is two EWMAs, one over gradients (β₁ = 0.9, memory ≈ 10 steps) and one over squared gradients (β₂ = 0.999, memory ≈ 1,000 steps). BatchNorm's "running mean" and "running var" are EWMAs with momentum 0.1. Model-weight EMA, used to stabilize diffusion models and LLM checkpoints, is the same formula applied to every parameter.

Real-life engineering ties

• SWE: SLO definitions (p99 latency < 300ms), capacity planning off percentile traffic, not average traffic.
• AI: every model evaluation metric is a statistic over a test set; a single accuracy number without variance/confidence interval is incomplete (Ch. 22).
• Network: jitter is literally the standard deviation of packet latency; SLA definitions for ISPs use percentile packet loss/latency.
• Security: baselining "normal" user/system behavior statistically is the foundation of anomaly-based intrusion detection — and its Achilles' heel (see "base rate fallacy," Ch. 16 and 26).

More real-life examples:

• Google's Site Reliability Engineering practice formalizes "error budgets" directly from percentile/SLO statistics — a team's freedom to ship risky changes is a function of how much of their p99-based budget remains this quarter.
• Monitoring systems like Prometheus and Datadog store latency as histograms specifically so they can compute p50/p95/p99 after the fact, rather than only an average — a direct product decision driven by the mean-vs-median problem shown above.
• Credit card fraud detection scores transactions against a customer's historical spending distribution and flags statistical outliers (large deviations from the customer's typical mean/variance) for review.
• Manufacturing and DevOps alike use "three-sigma" control charts (originally from Six Sigma manufacturing QA) to decide when a metric — defect rate or deploy failure rate — has drifted outside normal statistical variation and warrants investigation.

Common mistake

Reporting only the mean/average of a skewed distribution (like latency, which is almost always right-skewed with a long tail) — always pair it with a median or percentile, and look at the shape, not just one number.