Statistics turns data you happened to collect into claims about a world you did not observe; probability is the machinery that makes those claims auditable. This post is a single self-contained pass over the foundation — the axioms, random variables, the distribution zoo, the limit theorems, estimation, hypothesis testing (including the two definitions everyone botches), regression, and the exact point where classical statistics becomes machine learning. Written for engineers who already write code and now need to reason about uncertainty without hand-waving. Every number below has been checked and is reproducible with numpy and scipy.stats.

Learning Outcome

Apply the probability axioms, conditional probability, independence, the law of total probability and Bayes' theorem to a real diagnostic problem, and see the base-rate fallacy numerically.
Work with random variables through PMF / PDF / CDF, and compute expectation, variance, covariance and correlation from their algebraic properties.
Pick a distribution from the standard zoo by its generating mechanism, and state the law of large numbers and the central limit theorem precisely — including what the CLT does not claim.
Derive a maximum-likelihood estimator from scratch, and say correctly what a 95% confidence interval means.
Run the standard tests, state a p-value correctly, correct for multiple comparisons, and prefer effect sizes to bare significance.
Fit and diagnose linear and logistic regressions, and describe the bias–variance trade-off separating models built to explain from models built to predict.

Chapter 1 — Probability foundations

Terminology

Sample space (Ω) — the set of all possible outcomes. One die: Ω = {1,2,3,4,5,6}.
Event — any subset of Ω. "The roll is even" is {2,4,6}.
Probability measure (P) — assigns a number in [0,1] to each event, obeying the axioms.
Mutually exclusive (disjoint) — A ∩ B = ∅, the events cannot both occur. Exhaustive — a collection whose union is all of Ω.
Independent — knowing one event occurred tells you nothing about the other.

The axioms (Kolmogorov, 1933)

1
2
3
4
5
6
7 A1  Non-negativity       P(A) ≥ 0 for every event A
A2  Normalisation        P(Ω) = 1
A3  Countable additivity for pairwise disjoint A₁, A₂, ...
                         P(A₁ ∪ A₂ ∪ ...) = P(A₁) + P(A₂) + ...
Consequences:  P(Aᶜ) = 1 − P(A)    P(∅) = 0    A ⊆ B ⇒ P(A) ≤ P(B)
               P(A ∪ B) = P(A) + P(B) − P(A ∩ B)

Watch out: the addition rule loses its correction term only when the events are disjoint. P(A ∪ B) = P(A) + P(B) is a special case, not the general rule.

Counting

When every outcome in a finite Ω is equally likely, P(A) = |A| / |Ω| and the problem is pure counting. The only question is whether order matters: P(n,r) = n!/(n−r)! for ordered selections, C(n,r) = n!/(r!(n−r)!) for unordered, and n₁·n₂·...·n_k for k independent stages.

Three distinct medals among 10 runners (ordered) is P(10,3) = 10·9·8 = 720; a three-person committee from the same 10 (unordered) is C(10,3) = 720/3! = 120.
Five-card poker hands: C(52,5) = 2,598,960. A flush (any suit, straight flushes included) is 4·C(13,5)/C(52,5) = 5148/2,598,960 ≈ 0.00198 — about 0.198%.
Birthday problem: P(match) = 1 − (365/365)(364/365)...((365−n+1)/365) gives 0.5073 at n = 23 and 0.9704 at n = 50. 23 people beats a coin flip because they generate C(23,2) = 253 pairs, not 23 opportunities.

Conditional probability and independence

1
2
3
4 P(A | B) = P(A ∩ B) / P(B)                       defined when P(B) > 0
P(A ∩ B) = P(A | B)·P(B) = P(B | A)·P(A)
A ⊥ B  ⟺  P(A ∩ B) = P(A)·P(B)  ⟺  P(A | B) = P(A)

Roll two fair dice. A = "sum is 8" (5/36: 2-6, 3-5, 4-4, 5-3, 6-2), B = "first die is 3" (6/36), A ∩ B = {3-5} (1/36). So P(A|B) = (1/36)/(6/36) = 1/6 ≠ 5/36 — not independent.

Common pitfall: mutually exclusive and independent are near-opposites, not synonyms. Disjoint events with positive probability satisfy P(A|B) = 0 ≠ P(A), so they are maximally dependent.

Law of total probability and Bayes' theorem

If B₁, ..., B_k partition Ω then P(A) = Σᵢ P(A | Bᵢ)·P(Bᵢ). Two factories: F1 makes 70% of units at a 2% defect rate, F2 makes 30% at 5%, so P(defect) = 0.70·0.02 + 0.30·0.05 = 0.029.

1
2
3
4              P(B | A)·P(A)
P(A | B) = ───────────────────    with P(B) = P(B|A)·P(A) + P(B|Aᶜ)·P(Aᶜ)
                  P(B)

Read it as an update: prior P(A), likelihood P(B|A), posterior P(A|B). Which factory made a defective part? P(F1|defect) = 0.014/0.029 ≈ 0.483 and P(F2|defect) = 0.015/0.029 ≈ 0.517 — F2 makes under a third of the parts but is the more likely source of any given defect.

The medical test and the base-rate fallacy

A disease affects 0.5% of the population. A test has 99% sensitivity (P(+|D) = 0.99) and 95% specificity (so a 5% false-positive rate). You test positive. What is P(D|+)? Natural frequencies over 100,000 people make it obvious:

1
2
3
4
5
6
7
8 100,000 people
├──    500 have the disease (0.5%)   →   495 test positive  ← true positives
└── 99,500 do not                    → 4,975 test positive  ← false positives
Total positives = 5,470        P(D | +) = 495 / 5,470 ≈ 0.0905   →  9.05%
 
Same in Bayes form:
P(D | +) = (0.99·0.005) / (0.99·0.005 + 0.05·0.995) = 0.00495 / 0.0547 ≈ 0.0905

1
2
3 sens, spec, prev = 0.99, 0.95, 0.005
print(sens * prev / (sens * prev + (1 - spec) * (1 - prev)))   # 0.09049360146252285

A "99% accurate" test leaves you with a 9% chance of being ill, because the false positives come from a pool 199 times larger than the true positives. Ignoring P(D) when computing P(D|+) is the base-rate fallacy. Retest with an independent second test and the posterior becomes the new prior: of the 5,470, 495·0.99 = 490.05 diseased and 4,975·0.05 = 248.75 healthy test positive again, so P(D|++) = 490.05/738.80 ≈ 0.663 — 66.3%. That is why screening programmes always confirm.

Watch out: that second calculation assumes the two tests' errors are conditionally independent given disease status. Two tests using the same biochemical mechanism fail on the same patients, and the true posterior is lower than 66%.

Chapter 2 — Random variables

Terminology

Random variable (X) — a function from the sample space to the reals. A rule, not a storage slot.
PMF p(x) = P(X = x) — discrete only. PDF f(x) — continuous only; a density, not a probability. CDF F(x) = P(X ≤ x) — both; non-decreasing, right-continuous, 0 → 1.
Support — the values with non-zero probability or density. i.i.d. — independent and identically distributed, the default assumption for a sample.

1
2
3
4
5
6 Discrete                           Continuous
p(x) ≥ 0,  Σ p(x) = 1              f(x) ≥ 0,  ∫ f(x) dx = 1
P(a ≤ X ≤ b) = Σ_{a≤x≤b} p(x)      P(a ≤ X ≤ b) = ∫ₐᵇ f(x) dx
F(x) = Σ_{t≤x} p(t)                F(x) = ∫_{−∞}^{x} f(t) dt,  f(x) = F′(x)
P(X = a) = p(a)                    P(X = a) = 0  always  ⇒  P(X<a) = P(X≤a)

Common pitfall: because P(X = a) = 0 in the continuous case, P(X < a) and P(X ≤ a) agree there — but for a discrete variable they differ by p(a). Off-by-one errors in binomial and Poisson tails almost always trace back to this.

Expectation, variance, and their algebra

1
2
3
4
5
6
7
8
9
10
11 E[X] = Σ x·p(x)  or  ∫ x·f(x) dx        E[g(X)] = Σ g(x)p(x)  or  ∫ g(x)f(x) dx
Var(X) = E[(X − μ)²] = E[X²] − (E[X])²  SD(X) = √Var(X)   (units of X)
E[a] = a   E[aX + b] = a·E[X] + b   Var(a) = 0   Var(aX + b) = a²·Var(X)
E[X + Y] = E[X] + E[Y]      ← linearity holds even for DEPENDENT X, Y
E[XY] = E[X]·E[Y]           ← independence required
Var(X ± Y) = Var(X) + Var(Y) ± 2·Cov(X, Y)
 
One fair die:  E[X] = 3.5,  E[X²] = 91/6 ≈ 15.1667
               Var(X) = 91/6 − 3.5² = 35/12 ≈ 2.9167,  SD ≈ 1.7078
Two dice:      E[sum] = 7,  Var(sum) = 35/6 ≈ 5.833,   SD ≈ 2.415

For independent variables the covariance term vanishes and variances add in both cases — subtracting random variables does not subtract uncertainty.

Covariance and correlation

1
2
3 Cov(X,Y) = E[(X − μₓ)(Y − μᵧ)] = E[XY] − E[X]·E[Y]      ρ = Cov(X,Y)/(σₓ·σᵧ)
Sample:  r = Σ(xᵢ − x̄)(yᵢ − ȳ) / √( Σ(xᵢ − x̄)² · Σ(yᵢ − ȳ)² ),  −1 ≤ ρ, r ≤ 1

Covariance carries units (kg·cm, dollars·hours) and is hard to read; correlation is the unitless rescaling, hitting ±1 exactly when Y is an affine function of X.

Correlation measures linear association only. Let X take −2, −1, 1, 2 with equal probability and Y = X² — a perfect deterministic relationship. Then E[X] = 0 and E[XY] = E[X³] = (−8 − 1 + 1 + 8)/4 = 0, so Cov(X,Y) = 0 and ρ = 0. Zero correlation, total dependence. (np.corrcoef([-2,-1,1,2], [4,1,1,4])[0,1] returns 0 up to floating point; the continuous twin is X ~ Uniform(−1,1), Y = X².)

Watch out: ρ = 0 means "no linear trend", never "independent". Independence implies zero correlation; the converse is false except for jointly normal variables. Plot the scatter before trusting r. And correlation is still not causation — X → Y, Y → X, a common cause Z, and selection effects in how the data was gathered all produce the same r.

Chapter 3 — The distribution zoo

Every named distribution is a generating mechanism plus its consequences. Choose one by asking what process produced the number.

Distribution	Parameters	Mean	Variance	Mechanism / typical use
Bernoulli(p)	`p ∈ [0,1]`	`p`	`p(1−p)`	one trial, success or failure
Binomial(n,p)	`n`, `p`	`np`	`np(1−p)`	successes in `n` independent trials, fixed `p`
Geometric(p)	`p`	`1/p`	`(1−p)/p²`	trials until first success (support 1, 2, 3, ...)
Poisson(λ)	`λ > 0`	`λ`	`λ`	count of rare events in a fixed window
Uniform(a,b)	`a < b`	`(a+b)/2`	`(b−a)²/12`	no reason to prefer any value in a range
Exponential(λ)	`λ > 0`	`1/λ`	`1/λ²`	waiting time between Poisson events
Normal(μ,σ²)	`μ`, `σ²`	`μ`	`σ²`	sums of many small independent effects
Log-normal(μ,σ²)	`μ`, `σ²`	`exp(μ + σ²/2)`	`(exp(σ²)−1)·exp(2μ+σ²)`	products of effects: incomes, latencies, file sizes
Student t(ν)	`ν > 0`	`0` for `ν>1`	`ν/(ν−2)` for `ν>2`	a mean when `σ` is estimated from the data
Chi-square(k)	`k`	`k`	`2k`	sums of squared standard normals: count tables, variances
F(d₁,d₂)	`d₁`, `d₂`	`d₂/(d₂−2)` for `d₂>2`	long closed form	ratio of two variances: ANOVA, model comparison

Bernoulli, binomial, geometric

P(X = k) = C(n,k)·pᵏ·(1−p)ⁿ⁻ᵏ, valid when trials are independent and p constant. Ten independent requests each failing with probability 0.3:

1
2
3
4 P(X = 3) = C(10,3)·0.3³·0.7⁷ = 120 · 0.027 · 0.0823543 = 0.26683
P(X ≤ 3) = 0.02825 + 0.12106 + 0.23347 + 0.26683 = 0.64961
E[X] = 3     Var(X) = 10·0.3·0.7 = 2.1     SD ≈ 1.449

Geometric: P(X = k) = (1−p)^(k−1)·p. With p = 0.2 the expected number of attempts is 5 and P(X > 10) = 0.8¹⁰ ≈ 0.1074. Like the exponential it is memoryless — past failures do not improve your prospects.

Poisson and the binomial limit

P(X = k) = e^(−λ)·λᵏ/k! — rare events at a constant average rate λ in a fixed window. A support line averaging 3 calls/hour: P(X = 5) = e⁻³·243/120 = 0.10082, P(X = 0) = e⁻³ = 0.04979, P(X ≥ 1) = 0.95021.

If n → ∞ and p → 0 with np = λ fixed, Binomial(n,p) → Poisson(λ):

1
2
3
4 n = 1000, p = 0.002, λ = 2   Binomial   Poisson
              P(X = 0)       0.13506    0.13534
              P(X ≤ 1)       0.40573    0.40601

Three decimals of agreement. Rules of thumb: good when n ≥ 100 and p ≤ 0.01, usually acceptable from n ≥ 20 with p ≤ 0.05. The error is of order p, so it is p — not n — that sets the quality.

Trick: use Poisson when you know a rate but not a trial count — packets dropped per second, defects per square metre, 500-errors per minute. If you cannot say what n is, the binomial was never the model.

Uniform, exponential, normal

Uniform(a,b) has constant density 1/(b−a); its practical roles are the output of a random number generator and the input to inverse-transform sampling (F⁻¹(U) has CDF F). Exponential is the waiting time between Poisson events:

1
2
3
4
5
6
7
8 f(t) = λe^(−λt)    F(t) = 1 − e^(−λt)    P(T > t) = e^(−λt)
E[T] = 1/λ         Var(T) = 1/λ²
 
λ = 0.1 per minute (mean gap 10 min):
  P(T ≤ 5) = 1 − e^(−0.5) = 0.39347      P(T > 15) = e^(−1.5) = 0.22313
Memorylessness:  P(T > s+t | T > s) = P(T > t)
  P(T > 25 | T > 10) = e^(−2.5)/e^(−1.0) = e^(−1.5) = 0.22313 = P(T > 15)

The ten minutes already spent waiting are worth nothing. The exponential is the only continuous distribution with this property (geometric is its discrete twin).

Watch out: memorylessness is usually false for anything that ages. Disks, mechanical parts and user sessions have hazard rates that change over time — reach for a Weibull or a log-normal when the failure rate is not constant.

Normal — f(x) = (1/(σ√(2π)))·exp(−(x−μ)²/(2σ²)), standardised by Z = (X − μ)/σ. The empirical rule is 68.27% / 95.45% / 99.73% within 1, 2, 3 σ; the exact 95% multiplier is 1.96, not 2. Sums of independent normals stay normal: N(μ₁,σ₁²) + N(μ₂,σ₂²) = N(μ₁+μ₂, σ₁²+σ₂²). Heights as Normal(170, 8²) cm: P(X > 182) has z = 1.5 so 1 − Φ(1.5) = 0.06681; P(162 < X < 178) has z = ±1 so Φ(1) − Φ(−1) = 0.68269.

1
2
3
4
5
6
7
8 from scipy import stats
print(stats.binom.pmf(3, 10, 0.3), stats.binom.cdf(3, 10, 0.3))
# 0.26682793199999977 0.6496107184
print(stats.poisson.pmf(5, 3), stats.expon.sf(15, scale=10))
# 0.10081881344492448 0.22313016014842982
print(1 - stats.norm.cdf(182, 170, 8), stats.norm.ppf(0.975))
# 0.06680720126885807 1.959963984540054

The four support distributions

Log-normal — log X is normal. Arises from multiplicative accumulation, so it is right-skewed with a long tail. Response times, incomes and city sizes fit far better here than under a normal; the mean is pulled well above the median, so report both.
Student t(ν) — the sampling distribution of (x̄ − μ)/(s/√n) when σ is estimated. Heavier tails than normal, converging to it as ν → ∞ (already close by ν ≈ 30). Used for t-tests and small-sample intervals.
Chi-square(k) — a sum of k squared independent standard normals, used for goodness-of-fit and independence tests on counts and for inference about a variance. F(d₁,d₂) — a ratio of two independent chi-squares each over its df, used for ANOVA and for comparing nested regression models.

Chapter 4 — Sampling and the limit theorems

Terminology

Population — every unit you want to claim something about; described by parameters μ, σ, p.
Sample — what you measured; described by statistics x̄, s, p̂. A statistic is a function of random data, so it is itself a random variable.
Sampling distribution — the distribution of a statistic across all possible samples of size n. The central object of inference, and not the distribution of your data. Standard error (SE) is its standard deviation.

1
2
3
4
5 Population   Sample                        E[x̄] = μ     Var(x̄) = σ²/n
    μ        x̄ = (1/n) Σ xᵢ                SE(x̄) = σ/√n (use s/√n if σ unknown)
    σ²       s² = (1/(n−1)) Σ (xᵢ − x̄)²    σ = 15:  n=25 → 3.00    n=100 → 1.50
    p        p̂ = k/n                                n=400 → 0.75   n=1600 → 0.375

Quadrupling the sample halves the error. Precision is expensive and gets more expensive the more of it you already have — this single fact should govern every "should we collect more data?" conversation. For a proportion, SE(p̂) = √(p(1−p)/n), largest at p = 0.5.

Watch out: all of this assumes independent sampling. Clustered data — many measurements per user, repeated requests from one client — has an effective sample size far below n. Treating it as n gives standard errors that are too small and intervals that are too narrow. This is the most common way real analyses go wrong.

Law of large numbers

x̄ₙ → μ in probability (weak LLN) and almost surely (strong LLN). That is all it says, and it needs a finite mean — for a Cauchy distribution the sample mean never settles down.

Common pitfall: the gambler's fallacy is the LLN misapplied. After ten heads it does not predict compensating tails; it predicts the ten heads become negligible against a growing n. The imbalance is diluted, never corrected.

The central limit theorem

1
2
3
4
5
6
7 Let X₁, ..., Xₙ be i.i.d. with finite mean μ and finite variance σ². Then
 
     x̄ₙ − μ
    ────────  →  Normal(0, 1)  in distribution as n → ∞
     σ / √n
equivalently, x̄ₙ is approximately Normal(μ, σ²/n) for large n.

What it does say: the sampling distribution of the mean (or of any sum) becomes approximately normal regardless of the population shape, provided the variance is finite. What it does not say:

It does not say your data becomes normal. Ten thousand exponential measurements are still exponentially distributed; only the mean is approximately normal.
It does not license a normality assumption about individual observations — which is what regression and t-test assumptions are partly about.
It does not apply with infinite variance: the mean of n Cauchy draws is Cauchy for every n.
n ≥ 30 is a rule of thumb, not a theorem. Mild skew converges by n ≈ 15; strong skew may need hundreds. Bootstrap to check rather than reciting 30.
It says nothing about the tails at finite n. The centre converges fastest and extreme quantiles slowest — which matters if you care about p99 latency.

1
2
3
4
5
6
7
8 import numpy as np
rng = np.random.default_rng(0)
pop = rng.exponential(scale=1.0, size=1_000_000)       # heavily right-skewed
for n in (2, 5, 30, 200):
    means = rng.choice(pop, size=(20_000, n)).mean(axis=1)
    print(n, round(means.std(), 4), round(1 / np.sqrt(n), 4))   # sd vs σ/√n
# the histogram is visibly skewed at n=2, near-symmetric at n=30, normal at n=200

Worked example. Population Uniform(0,1): μ = 0.5, σ = √(1/12) ≈ 0.288675. With n = 100, SE = 0.0288675, so P(x̄ > 0.55) ≈ 1 − Φ(0.05/0.0288675) = 1 − Φ(1.7321) = 0.0416 — about 4.2%.

Chapter 5 — Estimation

Terminology

Estimator θ̂ — a rule (function of the sample) for guessing a parameter. Estimate — the number it produces on one sample.
Bias — E[θ̂] − θ. Zero bias means right on average across repeated samples.
Consistency — θ̂ₙ → θ as n → ∞. Efficiency — smallest variance among unbiased estimators; the Cramér–Rao bound is the floor.
Likelihood — L(θ) = P(data ; θ), a function of θ with the data fixed. It is not a probability distribution over θ.

Since MSE(θ̂) = E[(θ̂ − θ)²] = Var(θ̂) + Bias(θ̂)², bias and variance are different failure modes and a slightly biased estimator with much lower variance can beat an unbiased one — that observation is the seed of regularisation in Chapter 8.

Why sample variance divides by n − 1: using x̄ in place of the unknown μ makes Σ(xᵢ − x̄)² systematically too small, because x̄ is precisely the value minimising that sum for this sample. Bessel's correction compensates exactly, giving E[s²] = σ².

Watch out: numpy defaults to ddof=0 (biased, the MLE) and pandas to ddof=1, so the same column yields two different standard deviations depending on which library you asked — np.var([2,4,4,4,5,5,7,9]) is 4.0 while ddof=1 gives 4.5714. Set ddof explicitly. Note also that s is a biased estimator of σ even though s² is unbiased: unbiasedness does not survive nonlinear transformation.

Method of moments and maximum likelihood

Method of moments sets sample moments equal to population moments and solves. For Exponential(λ) the population mean is 1/λ, so 1/λ = x̄ gives λ̂ = 1/x̄. MLE instead picks the parameter making the observed data most probable — worked all the way through:

1
2
3
4
5
6
7
8
9
10 1. Likelihood        L(λ) = Πᵢ λ e^(−λxᵢ) = λⁿ · exp(−λ Σᵢ xᵢ)
2. Log-likelihood    ℓ(λ) = n·ln(λ) − λ·Σᵢ xᵢ    (same maximiser, sums not products)
3. Differentiate     dℓ/dλ = n/λ − Σᵢ xᵢ = 0
4. Solve             λ̂ = n / Σᵢ xᵢ = 1 / x̄
5. Confirm maximum   d²ℓ/dλ² = −n/λ² < 0 for all λ > 0  →  strictly concave ✓
Bernoulli, n trials with k successes:
   ℓ(p) = k·ln(p) + (n − k)·ln(1 − p)
   dℓ/dp = k/p − (n − k)/(1 − p) = 0
     ⇒ k(1 − p) = (n − k)p  ⇒  k − kp = np − kp  ⇒  p̂ = k/n

Here MLE and method of moments coincide; frequently they do not. The Bernoulli MLE being the sample proportion is a good sign the machinery does something sensible. MLEs are consistent, asymptotically normal, asymptotically efficient and invariant (the MLE of g(θ) is g(θ̂)), but not generally unbiased — the MLE of σ² divides by n.

1
2
3
4
5
6
7 import numpy as np
from scipy import optimize
data = np.random.default_rng(1).exponential(scale=1 / 2.5, size=500)   # true λ = 2.5
neg_ll = lambda lam: -np.sum(np.log(lam) - lam * data)
print(1 / data.mean(),                                                 # closed-form MLE
      optimize.minimize_scalar(neg_ll, bounds=(1e-6, 50), method="bounded").x)

Confidence intervals

1
2
3
4
5
6 σ known      x̄ ± z₍₁₋α/2₎ · σ/√n        Proportion   p̂ ± z₍₁₋α/2₎ · √(p̂(1−p̂)/n)
σ unknown    x̄ ± t₍₁₋α/2, n−1₎ · s/√n
n = 36, x̄ = 52, σ = 6 known:  SE = 1,    CI₉₅ = 52 ± 1.96·1   = (50.04, 53.96)
n = 25, x̄ = 52, s = 6:         SE = 1.2,  t₀.₉₇₅,₂₄ = 2.0639
                               CI₉₅ = 52 ± 2.4767 = (49.52, 54.48)

The t-interval is wider; that width is the price of estimating σ from the same 25 numbers. What "95% confidence" actually means — a statement about the procedure, not this interval:

If you repeated the whole experiment many times — new sample, same method — 95% of the intervals constructed this way would contain the true parameter. The 95% is a long-run coverage property of the recipe.

Once (49.52, 54.48) is computed, μ is either inside it or not; there is no randomness left to attach a probability to. Under the frequentist framework μ is a fixed unknown constant. Misreadings to avoid:

❌ "There is a 95% probability μ lies in (49.52, 54.48)." That is the Bayesian credible interval statement and it requires a prior — a legitimate thing to want, just not what this gives you.
❌ "95% of the data lies in this interval." That is a prediction or tolerance interval, and it is far wider.
❌ "95% of future sample means will fall in this interval." No — different interval, lower coverage. And ❌ "a wider interval means the estimate is wrong": a wider interval means less information, honestly reported.

1
2
3
4
5
6
7
8
9 import numpy as np
from scipy import stats
rng = np.random.default_rng(7); mu, sigma, n, hits = 100.0, 15.0, 25, 0
for _ in range(100_000):
    s = rng.normal(mu, sigma, n)
    lo, hi = stats.t.interval(0.95, n - 1, s.mean(), s.std(ddof=1) / np.sqrt(n))
    hits += lo <= mu <= hi
print(hits / 100_000)      # ≈ 0.95 — coverage is a property of the procedure

Chapter 6 — Hypothesis testing

Terminology

Null hypothesis (H₀) — the "nothing interesting happening" model, specific enough to compute probabilities under. Alternative (H₁) — what you conclude if H₀ is rejected.
Test statistic — a number from the sample whose distribution under H₀ is known. Significance level (α) — the Type I error rate chosen before seeing the data.
Type I error — rejecting a true H₀ (false positive, rate α). Type II — failing to reject a false H₀ (false negative, rate β). Power — 1 − β.
Effect size — how big the difference is, in units someone cares about.

The p-value, stated correctly

A p-value is the probability of obtaining a test statistic at least as extreme as the one observed, computed under the assumption that H₀ and all the model's assumptions are true.

Every clause carries weight. What it is not: ❌ the probability H₀ is true — that is P(H₀ | data), while a p-value is P(data this extreme | H₀), the same confused conditional as the medical test in Chapter 1; ❌ the probability the result is "due to chance"; ❌ 1 − P(H₁ true); ❌ a measure of effect size, since at n = 10,000,000 a trivial difference gives p < 0.001; ❌ a statement about whether the finding will replicate; ❌ something with a bright line at 0.05, since p = 0.049 and p = 0.051 are the same evidence.

A large p-value is also not evidence for H₀ — absence of evidence is not evidence of absence, especially at small n. "We failed to reject H₀" is the honest phrasing.

Mechanics. State H₀ and H₁ before seeing the data; choose α and a sample size giving adequate power; compute the statistic; compute the p-value from its null distribution; reject if p ≤ α; and always report the effect size with its confidence interval. A two-tailed test asks "is there a difference?", a one-tailed test asks "is it bigger?" and puts all of α in one tail — easier to reach significance in that direction, impossible in the other. Choosing the tail after seeing which way the data went doubles the real Type I rate.

Worked tests

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 z-test — σ = 15 known, H₀: μ = 100, n = 25, x̄ = 106
  SE = 15/√25 = 3        z = (106 − 100)/3 = 2.00
  two-tailed p = 2·(1 − Φ(2)) = 0.0455    CI₉₅ = 106 ± 1.96·3 = (100.12, 111.88)
  → rejects at α = 0.05, barely; the CI excludes 100 and is more informative
 
One-sample t — n = 20, x̄ = 5.4, s = 1.2, H₀: μ = 5
  SE = 1.2/√20 = 0.26833     t = 0.4/0.26833 = 1.4907,  df = 19
  two-tailed p = 0.1525      → do not reject; the 8% difference is not
                               resolvable with 20 observations
 
Two-sample t — A: n=12, x̄=7.8, s=1.4    B: n=12, x̄=6.5, s=1.6
  s²ₚ = (11·1.4² + 11·1.6²)/22 = 2.26      sₚ = 1.5033
  SE = sₚ·√(1/12 + 1/12) = 0.61373         t = 1.3/0.61373 = 2.118,  df = 22
  two-tailed p = 0.0457                    Cohen's d = 1.3/1.5033 = 0.865 (large)

Trick: default to Welch's t-test (equal_var=False) for two independent samples. It does not assume equal variances, costs almost nothing in power when they are equal, and avoids a pre-test that inflates the error rate.

Chi-square goodness-of-fit. A die rolled 60 times, observed [8, 9, 12, 11, 10, 10], expected 10 each: χ² = (4+1+4+1+0+0)/10 = 1.00, df = 5, p = 0.963 — no evidence of unfairness.

Chi-square independence. A 2×2 table with Eᵢⱼ = (row·col)/grand = 25 in every cell:

1
2
3
4
5             Success Failure Total     χ² = 4·(5²/25) = 4.00, df = 1, p = 0.0455
Treatment      30      20      50      Yates' correction:
Control        20      30      50        χ² = 4·(4.5²/25) = 3.24, p = 0.0719
Total          50      50     100

The correction flips the verdict — which is exactly why a bare "significant / not significant" is a poor summary.

1
2
3
4
5
6
7
8 import numpy as np
from scipy import stats
print(stats.ttest_ind_from_stats(7.8, 1.4, 12, 6.5, 1.6, 12, equal_var=True))
print(stats.chisquare([8, 9, 12, 11, 10, 10]))          # 1.0, 0.9625657732472964
t = np.array([[30, 20], [20, 30]])
print(stats.chi2_contingency(t, correction=False)[:2])  # (4.0, 0.045500263...)
print(stats.chi2_contingency(t, correction=True)[:2])   # (3.24, 0.071860638...)

One-way ANOVA decomposes variability as SS_total = SS_between + SS_within and forms F = MS_between/MS_within = (SS_between/(k−1)) / (SS_within/(N−k)). A significant F says "not all means are equal" — not which. Follow with pairwise comparisons and a correction (Tukey's HSD is built for it).

Common pitfall: running all C(k,2) pairwise t-tests instead of an ANOVA. With k = 5 that is 10 tests and roughly a 40% family-wise false positive rate.

Assumptions and nonparametric alternatives

Test	Question	Key assumptions	Nonparametric alternative
One-sample z	mean vs known value, `σ` known	independence; normal or large `n`	sign test
One-sample t	mean vs known value, `σ` unknown	independence; approx. normal	Wilcoxon signed-rank
Two-sample t (pooled)	two independent means	independence; normality; equal variances	Mann–Whitney U
Welch's t	two independent means	independence; normality; unequal variances OK	Mann–Whitney U
Paired t	mean of within-pair differences	differences approx. normal; pairs independent	Wilcoxon signed-rank
Chi-square GoF / independence	counts vs expected; two categoricals related	independent counts; expected cells ≥ 5	exact multinomial; Fisher's exact
One-way ANOVA	3+ group means all equal	independence; normality; equal variances	Kruskal–Wallis
Pearson correlation	linear association	bivariate normal, no strong outliers	Spearman rank correlation

Paired vs two-sample is the choice that most often costs power. If each subject is measured twice (before/after, A/B on the same user), analyse the n differences with a one-sample t-test — pairing removes between-subject variability, usually the largest noise source.

Power, sample size and multiple comparisons

1
2
3
4
5              2·(z₍₁₋α/2₎ + z₍₁₋β₎)²         with α = 0.05 and 80% power:
  n_group ≈ ───────────────────────         2·(1.95996 + 0.84162)² ≈ 15.70
                      d²                    so n ≈ 16/d² per group
  d = 0.2 (small) → n ≈ 393    d = 0.5 (medium) → n ≈ 63    d = 0.8 → n ≈ 25

Underpowered studies are doubly bad: they usually miss real effects, and when they do reach significance the estimated effect is inflated, because only the luckiest samples cross the threshold. That inflation is one mechanical driver of failed replications.

Run m independent tests at α = 0.05 and P(at least one false positive) = 1 − 0.95^m: 5.0% at m = 1, 22.6% at 5, 40.1% at 10, 64.2% at 20, 92.3% at 50. Bonferroni tests each hypothesis at α/m, controlling the family-wise error rate — simple, always valid, conservative. Benjamini–Hochberg controls the false discovery rate: sort p-values ascending, find the largest k with p₍ₖ₎ ≤ (k/m)·q, reject 1 through k.

1
2
3
4
5
6
7 m = 10, q = 0.05   k:       1      2      3      4      5      6
                   p₍ₖ₎:    0.001  0.008  0.012  0.021  0.033  0.049
                   (k/m)·q: 0.005  0.010  0.015  0.020  0.025  0.030
                   reject:  yes    yes    yes ←  no     no     no
                                        largest k satisfying the condition
BH rejects the first 3.  Bonferroni (0.05/10 = 0.005) rejects only the first 1.

1
2
3
4
5 from statsmodels.stats.multitest import multipletests   # pip install statsmodels
p = [0.001, 0.008, 0.012, 0.021, 0.033, 0.049, 0.13, 0.24, 0.55, 0.71]
print(multipletests(p, alpha=0.05, method="fdr_bh")[0])       # first 3 True
print(multipletests(p, alpha=0.05, method="bonferroni")[0])   # first 1 True

Use Bonferroni when a single false positive is expensive (a confirmatory clinical endpoint); use FDR when screening many candidates with a tolerable proportion of duds (genomics, feature screening, dashboards with dozens of metrics).

p-hacking and the replication crisis

Every test above is valid only if the analysis was fixed before the data was seen. The moment analytic decisions respond to the data, the reported p-value's null distribution is no longer the one you looked up — and this mostly happens without dishonest intent. Optional stopping: peeking and stopping when p < 0.05; with enough peeks, an A/B test on two identical variants crosses the line with probability approaching 1. Outcome switching: measuring ten metrics, reporting the one that worked. Subgroup mining: "significant among iOS users aged 25–34". Flexible exclusions: choosing which outliers to drop after seeing the effect on p. Model shopping: trying covariates and transformations until something is significant. The garden of forking paths: the same inflation without explicit re-testing, because the analysis would have differed had the data differed. HARKing: hypothesising after results are known and presenting exploratory work as confirmatory.

The replication crisis is the consequence. The Open Science Collaboration's 2015 Estimating the Reproducibility of Psychological Science re-ran 100 studies: 97% of originals reported significant results, 36% of replications did, with replication effect sizes averaging roughly half the originals. Comparable gaps appeared in cancer biology and experimental economics. What helps: pre-registration (hypotheses, outcomes, exclusions and analysis plan written down first — keep exploratory work, just label it); report effect sizes with confidence intervals, since "d = 0.31, 95% CI (0.05, 0.57)" carries size, direction and precision while "p = 0.03" carries almost nothing; power the study beforehand, because post-hoc power computed from the observed effect is a deterministic function of the p-value and adds no information; publish nulls, since the file-drawer problem makes the literature a biased sample; and share data and code, because reproducibility is the precondition for replicability.

Watch out: if you must peek at a running A/B test, use a method designed for it — group sequential designs with alpha spending, or always-valid sequential inference. Repeatedly applying a fixed-sample t-test is not one of those methods.

Chapter 7 — Regression

Terminology

Response y — what you model; predictor x — what you model it from. Residual — eᵢ = yᵢ − ŷᵢ, the part the model failed to explain (not the same as the unobservable error term).
Homoscedasticity — residual variance constant across fitted values; the opposite is heteroscedasticity.
Multicollinearity — predictors strongly linearly related to each other. Leverage — how far a point's predictor values sit from the centre; high-leverage points can dominate a fit.

Simple linear regression

1
2
3
4
5
6
7
8
9
10
11
12
13
14 Model    yᵢ = β₀ + β₁xᵢ + εᵢ,     εᵢ ~ i.i.d. Normal(0, σ²)
 
         Σ(xᵢ − x̄)(yᵢ − ȳ)     Sxy           sᵧ
  β̂₁ = ─────────────────── =  ───  =  r · ──        β̂₀ = ȳ − β̂₁·x̄
              Σ(xᵢ − x̄)²        Sxx           sₓ
 
x:  1   2   3   4   5   (x̄ = 3)        y:  2   4   5   4   5   (ȳ = 4)
Sxy = (−2)(−2) + (−1)(0) + (0)(1) + (1)(0) + (2)(1) = 6
Sxx = 4+1+0+1+4 = 10        Syy = 4+0+1+0+1 = 6
β̂₁ = 0.6      β̂₀ = 4 − 0.6·3 = 2.2         ŷ = 2.2 + 0.6x
fitted      2.8   3.4   4.0   4.6   5.2
residuals  −0.8   0.6   1.0  −0.6  −0.2
SSE = 2.40   SST = 6.00   R² = 1 − 2.4/6 = 0.60   r = 6/√60 = 0.7746, r² = 0.60 ✓

The fitted line always passes through (x̄, ȳ); stats.linregress on those five points returns slope 0.6, intercept 2.2, rvalue**2 = 0.6. Interpretation: β̂₁ = 0.6 means a one-unit increase in x is associated with a 0.6-unit increase in the mean of y — "causes" only if x was randomised. β̂₀ = 2.2 is the prediction at x = 0, meaningless whenever x = 0 is outside the data or physically impossible; centring x at its mean makes the intercept the prediction at average x.

R² and its limits. R² is the fraction of the variance of y explained by the model — useful and routinely over-read. It never decreases when you add a predictor, even a column of pure noise, so compare models with adjusted R² or out-of-sample error. A high R² does not mean the functional form is right: Anscombe's quartet is four datasets with identical x̄ = 9, ȳ = 7.50, r = 0.816, the same line ŷ = 3.00 + 0.50x and the same R² = 0.667 — yet one is linear, one a clean parabola, one a perfect line plus an outlier, and one a vertical stack plus one influential point. Conversely a low R² does not mean the model is useless: with irreducible noise (human behaviour, financial returns) an R² of 0.05 can carry a real, well-estimated, economically significant slope. And R² says nothing about causality, bias, or whether the model holds on new data.

Trick: always plot the data and the residuals. Summary statistics exist to compress, and what breaks a regression is exactly what compression discards.

Multiple regression and its assumptions

In y = β₀ + β₁x₁ + ... + β_p x_p + ε, each βⱼ is the expected change in y per unit of xⱼ holding the other predictors in the model fixed. That conditional clause is why a coefficient can flip sign when a variable is added — the classic Simpson's paradox situation.

Assumption	Meaning	Diagnostic	Failure symptom	Fix
Linearity	mean of `y` linear in predictors	residuals vs fitted	systematic curve / U shape	transform, polynomial or spline terms
Independence	errors uncorrelated across observations	residuals vs order; ACF; Durbin–Watson	runs of same sign; time pattern	time-series or mixed model; cluster-robust SEs
Homoscedasticity	constant error variance	residuals vs fitted; scale–location	funnel / megaphone shape	log-transform `y`; WLS; robust (HC) SEs
Normality of residuals	errors approximately normal	Q-Q plot of residuals	heavy tails, strong skew	transform; robust regression; lean on the CLT at large `n`

Independence matters most — violating it corrupts every standard error and p-value in the output, and no amount of n rescues it. Normality matters least — it is a condition on the residuals, never on y or the predictors, and moderate departures are survivable at decent n.

1
2
3
4
5
6
7
8
9
10
11
12 import statsmodels.api as sm
import matplotlib.pyplot as plt
 
model = sm.OLS(df["price"], sm.add_constant(df[["area", "rooms", "age"]])).fit()
print(model.summary())                   # coefficients, SEs, t, p, R², adj-R², F
print(model.get_influence().cooks_distance[0].max())     # influential points
 
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
ax[0].scatter(model.fittedvalues, model.resid)     # linearity + homoscedasticity
sm.qqplot(model.resid, line="45", fit=True, ax=ax[1])              # normality
ax[2].plot(model.resid.values, marker=".")                         # independence

Multicollinearity and categorical predictors

Strongly correlated predictors leave the model unable to attribute credit between them. Predictions stay fine; individual coefficients become unstable, with huge standard errors and signs that flip on refit. VIFⱼ = 1/(1 − R²ⱼ), where R²ⱼ comes from regressing xⱼ on all other predictors: VIF = 1 is uncorrelated, > 5 is worth investigating, > 10 is serious and the coefficients are not individually interpretable. Remedies: drop a redundant predictor, combine them into an index or principal component, collect more data, or use ridge regression (Chapter 8), designed for exactly this.

A categorical variable with k levels enters as k − 1 dummy columns; the omitted level is the reference and every coefficient is a contrast against it.

1
2
3 colour ∈ {red, green, blue}, reference = red:  y = β₀ + β₁·is_green + β₂·is_blue
   β₀ = mean(red)      β₁ = green − red      β₂ = blue − red

Including all k dummies alongside an intercept makes the design matrix singular — the dummy variable trap. pd.get_dummies(..., drop_first=True) and the C() operator in statsmodels formulas handle it.

Watch out: encoding an unordered category as integers (red=1, green=2, blue=3) tells the model that blue is three times red and green sits halfway between. Ordinal encoding is valid only for genuinely ordered levels with roughly equal spacing.

Logistic regression

With binary y, linear regression predicts probabilities outside [0,1] and has structurally heteroscedastic errors. Model the log-odds instead, fitted by maximum likelihood using the Bernoulli likelihood from Chapter 5:

1
2
3
4
5
6
7 odds = p/(1 − p)      logit(p) = ln(p/(1 − p)) = β₀ + β₁x₁ + ... + β_p x_p
inverse (logistic / sigmoid):   p = 1 / (1 + e^(−(β₀ + Σ βⱼxⱼ)))
logit(p) = −4 + 0.05·score
  score =  80:  logit = 0.0  →  p = 1/(1 + e⁰)  = 0.500
  score = 100:  logit = 1.0  →  p = 1/(1 + e⁻¹) = 0.731
  score = 120:  logit = 2.0  →  p = 1/(1 + e⁻²) = 0.881

βⱼ is the change in log-odds per unit of xⱼ; exponentiate for an odds ratio, so β₁ = 0.7 → OR = e^0.7 = 2.01 means each extra unit multiplies the odds by about 2.01, holding the others fixed (OR = 1 no effect, > 1 increases, < 1 decreases). Note that the same 20-point score change above moves the probability by 23 points then 15: the effect on the probability is nonlinear even though the effect on the log-odds is constant. In statsmodels, np.exp(sm.Logit(y, sm.add_constant(X)).fit().params) gives the odds ratios and np.exp(fit.conf_int()) their intervals.

Common pitfall: an odds ratio is not a risk ratio. They are close only when the outcome is rare (p < 0.1); for common outcomes the OR sits much further from 1 than the relative risk, so reading "twice as likely" from OR = 2 overstates the effect. Report absolute probabilities at representative predictor values too. Evaluate with the confusion matrix, precision/recall, ROC-AUC and calibration — and note that accuracy on an imbalanced problem is the base-rate fallacy in disguise: always predicting "no disease" is 99.5% accurate on the Chapter 1 population and useless.

Chapter 8 — Where statistics becomes machine learning

Everything above is mostly inference: estimate a parameter, quantify uncertainty, decide whether an effect is real. Machine learning uses the same mathematics to optimise a different target — out-of-sample prediction error. The methods overlap almost entirely; the objectives do not.

Terminology

Training error — error on the data used to fit; optimistically biased, always. Generalisation error — expected error on new data from the same distribution; the quantity you actually care about.
Overfitting — the model learned noise specific to the training set (train error low, generalisation error high). Underfitting — too rigid to capture the signal (both high). Hyperparameter — a knob not learned from the training data. Regularisation — a penalty on complexity, trading a little bias for a lot of variance.

1
2
3
4
5 E[(y − f̂(x))²]  =  Bias(f̂(x))²  +  Var(f̂(x))       +  σ²
                   ────────────     ────────────      ──────────────────
                   wrong shape      unstable across   irreducible noise —
                                    samples           no model removes it

High bias / low variance is a linear model on curved data: wrong the same way every time. Low bias / high variance is a degree-15 polynomial through 20 points: fits each training set beautifully and differently. The best model is rarely the most flexible; it is the one minimising the sum.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 import numpy as np
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
 
rng = np.random.default_rng(0)
X = rng.uniform(-3, 3, size=(60, 1)); y = np.sin(X).ravel() + rng.normal(0, 0.3, 60)
for degree in (1, 3, 9, 15):
    m = make_pipeline(PolynomialFeatures(degree), LinearRegression())
    cv = -cross_val_score(m, X, y, scoring="neg_mean_squared_error", cv=5).mean()
    m.fit(X, y)
    print(f"deg {degree:2d}  train {np.mean((y - m.predict(X))**2):.4f}  CV {cv:.4f}")
# training MSE falls monotonically with degree; CV MSE bottoms out mid-range
# then blows up — that turning point is overfitting

Train / validation / test

Train fits the parameters, validation chooses hyperparameters and model families, and test is touched exactly once to estimate generalisation error (60/20/20, or 80/20 with cross-validation inside the 80). k-fold cross-validation replaces a single validation split: partition into k folds, fit on k − 1, score the held-out one, rotate, average. k = 5 or 10 is standard, uses the data far more efficiently, and gives a variance estimate alongside the mean.

Watch out — leakage. Fit scalers, imputers and feature selectors inside the CV loop, on the training fold only; standardising the full dataset before splitting leaks test information and produces a score that will not survive production. sklearn.pipeline.Pipeline makes the correct thing the easy thing. Two more traps: split time-ordered data by time (TimeSeriesSplit), never randomly, or the model sees the future; and split grouped data by group (GroupKFold), or the same user lands on both sides.

Regularisation as a statistical idea

Ridge and lasso are not ML inventions bolted onto regression — they are MSE = Var + Bias² from Chapter 5 turned into an algorithm: accept a little bias, cut a lot of variance, lower the total.

1
2
3
4
5 OLS       minimise  Σ(yᵢ − ŷᵢ)²                              (no penalty)
Ridge     minimise  Σ(yᵢ − ŷᵢ)² + λ·Σ βⱼ²                    (L2)
Lasso     minimise  Σ(yᵢ − ŷᵢ)² + λ·Σ |βⱼ|                   (L1)
Elastic   minimise  Σ(yᵢ − ŷᵢ)² + λ(α·Σ|βⱼ| + (1−α)·Σ βⱼ²)

	Ridge (L2)	Lasso (L1)
Coefficients	shrunk toward 0, never exactly 0	some set exactly to 0
Feature selection	no	yes, automatically
Correlated predictors	keeps all, splits the weight	picks one arbitrarily, drops the rest
Bayesian reading	MAP with a Gaussian prior on `β`	MAP with a Laplace prior on `β`
Best when	many small effects; multicollinearity	few predictors matter; `p > n`

λ is chosen by cross-validation, never by looking at the test set. Both penalties require standardised predictors — otherwise a variable in millimetres is penalised a thousand times more than the same variable in metres — and the intercept is conventionally left unpenalised.

1
2
3
4
5
6
7
8
9
10 import numpy as np
from sklearn.linear_model import RidgeCV, LassoCV
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
 
ridge = make_pipeline(StandardScaler(), RidgeCV(alphas=np.logspace(-3, 3, 50), cv=5))
lasso = make_pipeline(StandardScaler(), LassoCV(cv=5, random_state=0, max_iter=10_000))
ridge.fit(X_tr, y_tr); lasso.fit(X_tr, y_tr)
print("lasso zeroed:", int(np.sum(lasso[-1].coef_ == 0)), "of", lasso[-1].coef_.size)

Common pitfall: the p-values printed for a regularised model are not valid in the usual sense — the penalty biases the coefficients on purpose. Selecting variables with lasso and then reporting ordinary OLS inference on the survivors badly understates the uncertainty. Honest p-values are a job for an inferential model on a separate dataset.

Inference vs prediction

	Inferential model	Predictive model
Question	Why? Does `x` affect `y`, by how much?	What next? What is `y` for this case?
Success criterion	unbiased coefficients, correct SEs	low error on unseen data
Complexity / bias	simple and interpretable; avoid bias	whatever generalises; trade bias for variance
Correlated predictors	a serious problem	mostly harmless
Discipline	pre-registration, fixed analysis plan	held-out test set touched once
Validation	assumption diagnostics, residual plots	cross-validated out-of-sample score
Typical output	"effect 0.6 units, 95% CI (0.31, 0.89)"	"AUC 0.87 on held-out data"
Risk	a wrong causal claim	silent accuracy decay under drift

The failure mode is using one where you needed the other. A gradient-boosted model with excellent AUC says nothing about whether raising the price causes churn — it will happily exploit any correlate, including consequences of churn. Conversely, a beautifully specified regression with every assumption satisfied may generalise poorly, because unbiasedness is not low prediction error. Decide which you are building before fitting anything: it changes the split, the model class, the metric, and what you may conclude.

Where the data comes from

None of this happens on clean arrays — the real pipeline is extract, aggregate, reshape, then model.

The aggregation stage is usually a database job. The MongoDB course's aggregation basics post covers $group, $bucket and the accumulators that compute group means, counts and variances server-side — where the group summaries feeding an ANOVA or a chi-square table should come from.
Time-ordered measurements bring their own independence problems; the time-series at extreme scale post covers bucketing and windowing, and the statistical point is that consecutive observations are autocorrelated, so the effective n is below the row count.
Vector search and RAG are cosine similarity plus nearest-neighbour retrieval — geometry built on the covariance ideas of Chapter 2; see Atlas Search and vector search and the production RAG build, and note that evaluating a retrieval system is a hypothesis-testing problem carrying every multiple-comparison hazard from Chapter 6. The Python side starts with data structures — dict, list, comprehensions — which is what pandas is built on and what you reach for when reshaping records before they hit numpy.

Practice Exercises

Exercise 1 — Base rates

A screening test has 98% sensitivity and 90% specificity. Compute P(disease | positive) for prevalence 0.1%, 1% and 10%, then find the prevalence at which a positive result first becomes more likely than not.

Solution

1
2
3
4
5
6
7
8
9 sens, spec = 0.98, 0.90
posterior = lambda prev: sens * prev / (sens * prev + (1 - spec) * (1 - prev))
print([round(posterior(p), 4) for p in (0.001, 0.01, 0.10)])
# [0.0097, 0.0901, 0.5213]
 
# posterior > 0.5 ⟺ odds > 1 ⟺ prev/(1−prev) · sens/(1−spec) > 1
print(1 / (1 + sens / (1 - spec)))    # 0.0926 → below ~9.3% prevalence a
                                      # positive still leaves disease unlikely

Exercise 2 — Variance algebra

X and Y are independent with E[X] = 3, Var(X) = 4, E[Y] = 5, Var(Y) = 9. Compute E[2X − 3Y + 1] and Var(2X − 3Y + 1). Then redo the variance assuming Cov(X, Y) = −2.

Solution

1
2
3
4
5 E[2X − 3Y + 1] = 2·3 − 3·5 + 1 = −8
Independent:   Var = 2²·4 + (−3)²·9 = 16 + 81 = 97
Cov = −2:      Var(aX+bY) = a²Var(X) + b²Var(Y) + 2ab·Cov(X,Y)
                          = 97 + 2·(2)·(−3)·(−2) = 121   (the +1 never matters)

Exercise 3 — Choose the distribution

Name the distribution and parameters for each, justifying from the generating mechanism: (1) heads in 50 fair coin flips; (2) typos on a page averaging 1.4 per page; (3) minutes to the next request at 12 requests/hour, memoryless; (4) interviews until the first acceptance at 15% per candidate; (5) end-to-end service latency formed as a product of independent multiplicative slowdowns.

Solution

1
2
3
4
5
6
7 1. Binomial(n = 50, p = 0.5)   fixed trials, constant p    mean 25, var 12.5
2. Poisson(λ = 1.4)            rate in a window, no n      mean = var = 1.4
3. Exponential(λ = 0.2 /min)   memoryless waiting time     mean 5, var 25
4. Geometric(p = 0.15)         trials until first success  mean ≈ 6.67
5. Log-normal                  products → normal on the log scale; right-skewed,
                               so report median and p95, not just the mean

Exercise 4 — CLT by simulation

Draw 20,000 samples of size n ∈ {1, 2, 5, 30} from Uniform(0,1). Histogram the sample means and compare their standard deviation to σ/√n with σ = √(1/12). At which n does the histogram stop looking triangular? Then repeat with rng.standard_cauchy and explain why the spread never shrinks.

Exercise 5 — Derive an MLE

Derive the MLE of λ for n i.i.d. observations from Poisson(λ), including the second-derivative check.

Solution

1
2
3
4
5 L(λ)  = Πᵢ e^(−λ)·λ^(xᵢ)/xᵢ!      ℓ(λ) = −nλ + (Σᵢ xᵢ)·ln(λ) − Σᵢ ln(xᵢ!)
dℓ/dλ = −n + (Σᵢ xᵢ)/λ = 0    ⇒    λ̂ = (Σᵢ xᵢ)/n = x̄
d²ℓ/dλ² = −(Σᵢ xᵢ)/λ² < 0     ⇒    maximum ✓
Σ ln(xᵢ!) is constant in λ and drops out of the derivative entirely.

Exercise 6 — Say what the interval means

A study reports "mean improvement 4.2 points, 95% CI (0.3, 8.1), p = 0.037." Which statements are correct? (1) There is a 95% probability the true mean lies in (0.3, 8.1). (2) Repeated many times, about 95% of intervals built this way would contain the true mean. (3) There is a 3.7% probability the null is true. (4) Data would be at least this extreme in 3.7% of studies if the true effect were exactly zero. (5) It is significant at α = 0.05, so the effect is important. (6) The interval is wide relative to the effect, so the study is consistent with anything from negligible to large.

Solution

1
2
3
4
5
6
7
8
9 1. ✗  Bayesian credible-interval phrasing; a CI describes the procedure.
2. ✓  The correct long-run coverage statement.
3. ✗  A p-value is P(data this extreme | H₀), not P(H₀ | data).
4. ✓  The correct p-value statement.
5. ✗  Significance is detectability, not importance — importance is decided by
      the effect size in domain units.
6. ✓  0.3 to 8.1 is nearly a 27-fold range: "positive effect, magnitude poorly
      determined, needs a larger study".

Exercise 7 — Fit and diagnose

Simulate y = 3 + 2x + ε with ε ~ Normal(0,1), n = 100. Fit with statsmodels and confirm the coefficients recover the truth within their standard errors. Then break one assumption at a time and note which diagnostic reacts: multiply the errors by x (heteroscedasticity → funnel in residuals-vs-fitted and wrong SEs; refit with cov_type="HC3"); add a 0.5x² term (nonlinearity → U-curve in residuals-vs-fitted); draw ε from a t with 2 df (heavy tails → Q-Q plot bending at both ends); generate ε as an AR(1) series with φ = 0.9 (dependence → Durbin–Watson far below 2 and visible runs in residuals-vs-order).

Cheat-sheet

What	How
Union / complement / conditional	`P(A ∪ B) = P(A) + P(B) − P(A ∩ B)`; `P(Aᶜ) = 1 − P(A)`; `P(A given B) = P(A ∩ B)/P(B)`; independent ⟺ `P(A ∩ B) = P(A)P(B)`
Bayes / counting	`P(A given B) = P(B given A)·P(A)/P(B)` — always check the base rate; `P(n,r) = n!/(n−r)!` ordered, `C(n,r) = n!/(r!(n−r)!)` unordered
Variance toolkit	`Var(X) = E[X²] − (E[X])²`; `Var(aX+b) = a²Var(X)`; `Var(X ± Y) = Var(X) + Var(Y)` if independent
Correlation	`ρ = Cov(X,Y)/(σₓσᵧ)` in `[−1,1]` — linear association only
Discrete zoo	Binomial `np`, `np(1−p)`; Poisson `mean = var = λ`; Geometric `1/p`; binomial → Poisson when `n` large, `p ≤ 0.01`, `λ = np`
Continuous zoo	Exponential `1/λ`, `P(T>t) = e^(−λt)`, memoryless; Normal 68/95/99.7, exact 95% is `1.96`
Standard error / CLT	`SE = σ/√n`, so quadruple `n` to halve it; the CLT normalises the sampling distribution of the mean, not your data
Sample variance / MLE	divide by `n − 1` (`numpy` needs `ddof=1`); MLE = write `L`, take `log`, differentiate, solve, check `ℓ″ < 0`
CI	`x̄ ± t₍₁₋α/2, n−1₎·s/√n`; 95% of intervals from this procedure cover `θ` — not `P(θ ∈ this one)`
p-value	`P(statistic at least this extreme given H₀ true)` — never `P(H₀ given data)`
Errors and power	Type I `α`, Type II `β`, power `1 − β`; `n ≈ 16/d²` per group for 80% power
Test picker	two independent means → Welch's t; paired → t on the differences; counts → chi-square; 3+ groups → ANOVA then Tukey
Multiple tests	Bonferroni `α/m` (family-wise) or Benjamini–Hochberg (FDR)
Regression	`β̂₁ = Sxy/Sxx = r·(sᵧ/sₓ)`; line passes through `(x̄, ȳ)`; `R² = 1 − SSE/SST`
Diagnostics	residuals vs fitted, Q-Q, residuals vs order; `VIF > 10` serious; `exp(β)` is an odds ratio, not a risk ratio
ML side	`MSE = Bias² + Var + noise`; ridge shrinks, lasso zeroes; fit transforms inside the CV fold

Where to go deeper

Bayesian statistics — Chapter 1's theorem generalises to full posteriors over parameters, which is where the "95% probability the parameter is in this interval" statement legitimately lives. Gelman et al., Bayesian Data Analysis; implement with PyMC or Stan.
Causal inference — regression coefficients are associations; turning them into causal claims needs potential outcomes, DAGs, backdoor adjustment, instrumental variables, difference-in-differences, regression discontinuity.
Experimental design and time series — allocation, blocking, randomisation and stratification (directly applicable to A/B testing, sequential designs and CUPED-style variance reduction); then ARIMA, state-space models and cross-validation that respects time ordering, for the case where Chapter 7's independence assumption fails by construction.
Generalised linear models — logistic regression is one member of a family, alongside Poisson regression for counts, gamma for skewed positives and ordinal models for ranked outcomes. Hierarchical / mixed-effects models — the principled answer to clustered data; partial pooling is regularisation derived from the data structure. Resampling — the bootstrap gives intervals for statistics with no tractable sampling distribution, and permutation tests give exact p-values without distributional assumptions.
Statistical learning — An Introduction to Statistical Learning is the natural sequel to Chapter 8, with The Elements of Statistical Learning as the deeper version.

References

Kolmogorov, A. N. — Foundations of the Theory of Probability (1933; Chelsea translation). Source of the Chapter 1 axioms.
Casella, G. & Berger, R. L. — Statistical Inference, 2nd ed. (the standard graduate reference for estimation and testing); Wasserman, L. — All of Statistics (compressed and rigorous, written for computer scientists).
James, G., Witten, D., Hastie, T. & Tibshirani, R. — An Introduction to Statistical Learning (free PDF).
Hastie, T., Tibshirani, R. & Friedman, J. — The Elements of Statistical Learning (free PDF); Hernán, M. A. & Robins, J. M. — Causal Inference: What If (free PDF).
Wasserstein, R. L. & Lazar, N. A. — The ASA Statement on p-Values, The American Statistician, 2016.
Open Science Collaboration — Estimating the Reproducibility of Psychological Science, Science 349 (2015); Benjamini, Y. & Hochberg, Y. — Controlling the False Discovery Rate, JRSS-B 57 (1995); Anscombe, F. J. — Graphs in Statistical Analysis, The American Statistician 27 (1973).
SciPy stats reference · statsmodels · scikit-learn user guide · NumPy random sampling

Planned future extensions:

A dedicated Bayesian inference post — conjugate priors, MCMC and posterior predictive checks, using the Chapter 1 medical-test example as the running thread.
A causal inference companion covering DAGs, confounding, and why "control for everything" is a bad strategy.
A practical A/B testing post: sequential designs, variance reduction, guardrail metrics, and what to do when someone peeks.
A time-series post extending Chapter 7 to autocorrelated data, pairing with the MongoDB time-series material, plus worked notebooks reproducing every numeric example here as runnable cells.