20 - Statistics: Estimating, Comparing, Deciding

Machine Learning Engineer course · Module 02 · Mathematics and statistics · Lesson 5

Goal. Move from probability to statistics: from the observed data go back to the unknown parameters that generated them, quantifying how far one can trust them. At the end you will be able to summarise a sample without being fooled, explain (with the proof) why the sample variance divides by n − 1, estimate a parameter with maximum likelihood, build confidence intervals (also with the bootstrap), read a hypothesis test correctly, and derive from scratch the least-squares regression line.

Myth to dispel. "The mean of my sample is the true mean." It is only the best estimate you have, and it changes from one sample to another. Statistics does not eliminate uncertainty: it measures it.

Introduction

Probability goes from the model to the data: "if the coin is fair, what will 100 tosses look like?". Statistics takes the opposite path: "I saw 62 heads in 100 tosses, what can I say about the coin?". The basic terms, explained here at first use:

• Population: the set of everything you want to talk about (all the customers, all the possible measurements).
• Sample: the part you actually observed, x₁, …, xₙ. It must be chosen at random (and independently) for the results to hold.
• Parameter: a number that describes the population (the true mean μ); it is fixed and unknown.
• Statistic: a number computed from the sample (the sample mean x̄); it is random, because it depends on the sample drawn.
• Estimator: the recipe (the formula) with which an estimate of a parameter is obtained from a sample.

The five items sit on two planes: parameter and population are fixed and unknown, statistic and sample are random because they change at every draw. All inferential statistics lives in the passage from one plane to the other, and every formula in this lesson is a way of saying how far the statistic can stray from the parameter.

The analogy: the soup. To know whether the soup is salty you do not drink it all: you taste a spoonful. The spoonful is the sample, the pot the population. It works if you have stirred well (random sample) and if the spoonful is big enough. And you know that two different spoonfuls will give slightly different answers: how different, statistics tells.

The path starts from descriptive statistics, goes through estimators, bias and the n − 1 correction, maximum likelihood, confidence intervals and the bootstrap, hypothesis tests and least-squares regression, and ends with the difference between correlation and causation. It rests on the central limit theorem and on conditional probability.

Descriptive statistics

Centre: the mean x̄ = (1/n) Σ xᵢ (the centre of mass), the median (the value in the middle after sorting), the mode (the most frequent). Dispersion: variance s² and standard deviation s; the interquartile range (IQR, the distance between the first and the third quartile, that is between the value below which 25% of the data lie and the one below which 75% lie) which ignores the tails.

The difference between mean and median is the quickest test to understand the shape: if the mean is much larger than the median, there is a right tail. Let us show it with a synthetic example of monthly incomes (log-normal distribution, long-tailed), to which we add a single extreme value.

import numpy as np
import pandas as pd

rng = np.random.default_rng(10)
redditi = rng.lognormal(mean=7.4, sigma=0.5, size=1000)          # typical values of a few thousand
con_estremo = np.append(redditi, 5_000_000.0)                       # a single outlier

def riassunto(v):
    q1, q3 = np.percentile(v, [25, 75])
    return {"mean": v.mean(), "median": np.median(v), "std.dev": v.std(ddof=1), "IQR": q3 - q1}

tab = pd.DataFrame({"without outlier": riassunto(redditi), "with one outlier": riassunto(con_estremo)}).round(0)
print(tab)
print("skewness without outlier:", round(float(pd.Series(redditi).skew()), 2))

Output obtained:

         without outlier  with one outlier
mean              1786.0            6779.0
median            1581.0            1582.0
std.dev            923.0          157981.0
IQR               1082.0            1082.0
skewness without outlier: 1.48

What the code does. It generates 1000 log-normal incomes (the exponential of a normal variable, hence always positive and with a right tail), then adds one of five million. The function riassunto computes mean, median, standard deviation (with ddof=1, the correction we explain shortly) and IQR. The comparison of the two columns shows who holds up: the mean and the standard deviation explode because of a single observation, median and IQR barely move. They are robust indicators. Which to use depends on the question: the total spend is known from the mean, the typical case from the median.

Estimators and the n − 1 correction

Bias, variance, squared error

An estimator θ̂ of a parameter θ is evaluated with three criteria:

• Bias: E[θ̂] − θ. An unbiased estimator errs "to the right and to the left" equally, on average.
• Variance: how much θ̂ oscillates from one sample to another.
• Mean squared error: MSE = E[(θ̂ − θ)²] = bias² + variance (proved by adding and subtracting E[θ̂] inside the square: the cross product vanishes).

The sample mean is unbiased (E[x̄] = μ, by the linearity of the expected value) and has variance σ²/n. Its standard deviation, σ/√n, is called the standard error: it says by how much the sample mean oscillates.

Why the sample variance divides by n − 1

The natural candidate to estimate σ² is (1/n) Σ (xᵢ − x̄)². But it is biased downwards. The intuition: x̄ is the "tailor-made" centre of this sample, it is the point that minimises the squared deviations; the data are therefore always closer to x̄ than to the true μ. So the deviations around x̄ underestimate those around μ.

Proof. We add and subtract μ: xᵢ − x̄ = (xᵢ − μ) − (x̄ − μ). Squaring and summing over i (the cross product is −2(x̄−μ)·Σ(xᵢ−μ) = −2n(x̄−μ)², because Σ(xᵢ−μ) = n(x̄−μ)):

Σ (xᵢ − x̄)² = Σ (xᵢ − μ)² − n (x̄ − μ)².

Take the expected value: E[Σ(xᵢ−μ)²] = nσ² (each term has mean σ²) and E[n(x̄−μ)²] = n · Var(x̄) = n · σ²/n = σ². Hence

E[Σ (xᵢ − x̄)²] = nσ² − σ² = (n − 1)σ².

Dividing by n − 1 (and not by n) gives an unbiased estimator: s² = (1/(n−1)) Σ (xᵢ − x̄)². The difference 1 is the "degree of freedom" consumed to estimate x̄ from the same data.

rng = np.random.default_rng(1)
sigma2_vera = 4.0                       # true variance of the population: N(0, 4)
n, ripetizioni = 5, 200_000
campioni = rng.normal(0, np.sqrt(sigma2_vera), size=(ripetizioni, n))

stima_n = campioni.var(axis=1, ddof=0)              # divides by n
stima_n1 = campioni.var(axis=1, ddof=1)             # divides by n - 1
print("true variance:", sigma2_vera)
print("mean of the estimates, divisor n    :", round(float(stima_n.mean()), 3), "| theory (n-1)/n * σ² =", (n - 1) / n * sigma2_vera)
print("mean of the estimates, divisor n - 1:", round(float(stima_n1.mean()), 3))

print("MSE divisor n    :", round(float(((stima_n - sigma2_vera) ** 2).mean()), 3))
print("MSE divisor n - 1:", round(float(((stima_n1 - sigma2_vera) ** 2).mean()), 3))

Output obtained:

true variance: 4.0
mean of the estimates, divisor n    : 3.19 | theory (n-1)/n * σ² = 3.2
mean of the estimates, divisor n - 1: 3.988
MSE divisor n    : 5.764
MSE divisor n - 1: 7.983

What the code does. It repeats 200,000 times the experiment "draw 5 values from a normal with true variance 4 and estimate the variance": each row of campioni is a sample, var(axis=1, ...) computes the variance of each row with the two divisors (in NumPy ddof = delta degrees of freedom, the number subtracted from n). The mean of the estimates with divisor n falls close to 3.2 = (4/5)·4, confirming the computed bias; with n − 1 it falls close to 4. The last two lines report a detail not to forget: the unbiased estimator is not the one with the smaller MSE. There is often a compromise between bias and variance: it will return forcefully throughout machine learning (it is the compromise between under- and overfitting).

Maximum likelihood

The analogy: the detective

You see certain footprints on the floor and must guess the shoe size of the suspect. For each candidate shoe size you can ask yourself: "if the culprit wore that size, how likely would it be to leave footprints like these?" You choose the size that makes the observed footprints most plausible. This is the principle of maximum likelihood (Fisher, 1922).

Formally. Given the independent sample x₁, …, xₙ with density (or mass) f(x; θ), the likelihood is, as a function of θ with the data fixed,

L(θ) = Π f(xᵢ; θ),

the product because the observations are independent. One maximises L, or better its logarithm ℓ(θ) = Σ ln f(xᵢ; θ) (the logarithm turns the product into a sum, does not change where the maximum is and prevents a product of many small numbers from becoming 0 in floating point: lesson 4 of module 01). The maximum of a smooth function is found where its slope (the derivative, which we will define properly in the calculus lesson) is zero.

Example 1: the coin (Bernoulli)

You observe k heads in n tosses. The log-likelihood is ℓ(p) = k ln p + (n − k) ln(1 − p). Zero slope:

k/p − (n − k)/(1 − p) = 0 ⇒ k(1 − p) = (n − k)p ⇒ p̂ = k/n.

The estimate is the observed frequency: the method justifies the intuition.

Example 2: the normal

For N(μ, σ²) the log-likelihood is ℓ = −(n/2) ln(2πσ²) − Σ(xᵢ−μ)²/(2σ²). With respect to μ, setting the slope to zero gives Σ(xᵢ − μ) = 0, hence μ̂ = x̄. With respect to σ², it gives σ̂² = (1/n) Σ(xᵢ − x̄)²: the estimator with divisor n, the biased one from before. Maximum likelihood does not guarantee an unbiased estimator, but guarantees good properties for large n (consistency).

from scipy.optimize import minimize_scalar

k, n = 62, 100
griglia = np.linspace(0.01, 0.99, 981)
loglik = k * np.log(griglia) + (n - k) * np.log(1 - griglia)
print("maximum on the grid at p =", round(float(griglia[np.argmax(loglik)]), 3), "| formula k/n =", k / n)

res = minimize_scalar(lambda p: -(k * np.log(p) + (n - k) * np.log(1 - p)), bounds=(0.001, 0.999), method="bounded")
print("numerical optimiser:", round(float(res.x), 4))

x = np.random.default_rng(2).normal(5.0, 2.0, size=300)
print("μ̂ (mean)            :", round(float(x.mean()), 4))
print("σ̂² with divisor n   :", round(float(x.var(ddof=0)), 4), "| with n-1:", round(float(x.var(ddof=1)), 4))

Output obtained:

maximum on the grid at p = 0.62 | formula k/n = 0.62
numerical optimiser: 0.62
μ̂ (mean)            : 4.893
σ̂² with divisor n   : 4.1309 | with n-1: 4.1447

What the code does. It computes the log-likelihood of the coin on a grid of values of p and looks for its maximum: it falls on 0.62 = 62/100, like the formula. Then the same problem is solved by a numerical optimiser (one minimises −ℓ, because optimisers look for minima): it is how one works when a closed formula does not exist, for example for logistic regression. The last part shows the estimates of a normal sample: the mean coincides with μ̂, and the two variances differ by a factor n/(n−1) (almost nothing at n = 300: the problem is serious for small samples).

Confidence intervals

A point estimate (x̄ = 4.87) does not say how much it is worth. A 95% confidence interval is an interval [a, b] built with a procedure such that, repeating it on many samples, 95% of the intervals obtained contain the true parameter.

Mind the reading. The 95% describes the procedure, not the single interval: once computed, [a, b] contains μ or it does not (in the drawing, 19 out of 20 do). Saying "μ has a 95% probability of lying in [a, b]" is, in the frequentist setting of this course, an incorrectness of language.

Construction for the mean

By the CLT (lesson 4), x̄ is approximately normal with standard deviation σ/√n: so P(|x̄ − μ| ≤ 1.96 σ/√n) ≈ 0.95, and the interval is x̄ ± 1.96 · σ/√n. But σ is almost always unknown: it is replaced by s. This substitution adds uncertainty (especially with small n), and William Gosset, signing as "Student" (1908), showed that the quantity (x̄ − μ)/(s/√n) follows a distribution with heavier tails than the normal, the Student t with n − 1 degrees of freedom. The correct interval is x̄ ± t · s/√n, with t the quantile of the t.

from scipy import stats

rng = np.random.default_rng(3)
n, prove, mu = 10, 40_000, 0.0
campioni = rng.normal(mu, 1.0, size=(prove, n))
medie = campioni.mean(axis=1)
s = campioni.std(axis=1, ddof=1)
se = s / np.sqrt(n)

z = 1.96
t = stats.t.ppf(0.975, df=n - 1)
copertura_z = (np.abs(medie - mu) <= z * se).mean()
copertura_t = (np.abs(medie - mu) <= t * se).mean()
print(f"t quantile with {n - 1} degrees of freedom: {t:.3f} (against 1.960 for the normal)")
print(f"coverage with 1.96 (wrong for small n): {copertura_z:.3f}")
print(f"coverage with Student's t              : {copertura_t:.3f}")

Output obtained:

t quantile with 9 degrees of freedom: 2.262 (against 1.960 for the normal)
coverage with 1.96 (wrong for small n): 0.919
coverage with Student's t              : 0.951

What the code does. It repeats 40,000 times the experiment "sample of 10 observations from a normal with μ = 0", computes for each the mean and the standard error s/√n and checks in how many cases the interval contains μ: the coverage. With the normal multiplier 1.96 it covers less than the promised 95% (about 92%): the interval is too narrow because it ignores the uncertainty of s. With the quantile of the t (2.262 for 9 degrees of freedom) the coverage returns to 95%. For large samples the two quantiles coincide, for small ones they do not.

The bootstrap: when there is no formula

How do you build an interval for the median, or for a ratio, or for a model metric? For many statistics there is no simple formula. The bootstrap (Efron, 1979) offers a universal recipe:

1. You have a sample of n observations. Treat it as if it were the population.
2. Draw with replacement n observations from it: a bootstrap sample.
3. Compute on it the statistic you care about.
4. Repeat thousands of times: you obtain a distribution of the statistic.
5. The 95% percentile interval lies between the 2.5th and the 97.5th percentile of that distribution.

The idea is that the sample is the best description of the population that you have: resampling it simulates "what would have happened with other samples". The order of the five steps is mandatory, because each uses the result of the previous one, and step 2 is the one that makes the difference: by drawing with replacement every bootstrap sample has the same size as the original but a different composition, and it is the variability among these compositions that estimates the uncertainty.

rng = np.random.default_rng(4)
campione = rng.lognormal(7.4, 0.5, size=60)                  # 60 synthetic incomes

mediane = np.array([np.median(rng.choice(campione, size=campione.size, replace=True)) for _ in range(10_000)])
lo, hi = np.percentile(mediane, [2.5, 97.5])
print(f"sample median: {np.median(campione):.0f}")
print(f"95 % bootstrap interval for the median: [{lo:.0f}, {hi:.0f}]")

# fraction of distinct observations in a bootstrap sample
n = campione.size
distinti = np.mean([len(np.unique(rng.integers(0, n, n))) / n for _ in range(5_000)])
print(f"mean fraction of distinct observations: {distinti:.4f} | 1 - 1/e = {1 - np.exp(-1):.4f}")

Output obtained:

sample median: 1515
95 % bootstrap interval for the median: [1361, 1789]
mean fraction of distinct observations: 0.6359 | 1 - 1/e = 0.6321

What the code does. rng.choice(campione, size=n, replace=True) is the bootstrap sample (draw with replacement); a list comprehension produces 10,000 of them and computes their median. The 2.5 and 97.5 percentiles delimit the interval. The last part verifies an elegant fact: in a bootstrap sample the fraction of distinct observations is about 63.2%. Why: the probability that a given observation is never drawn in n draws is (1 − 1/n)ⁿ, which for large n tends to e⁻¹ ≈ 0.368 (a notable limit that we will see in the calculus lesson); hence about 1 − 0.368 = 63.2% of distinct observations are present, the others are repeated. This share is that of the "out-of-bag" observations used by random forests to evaluate themselves.

Trade-off. The bootstrap is powerful but not magic: with very small samples, or for statistics sensitive to extremes (the maximum), it can work badly; and it does not correct a non-representative sample.

Hypothesis tests and p-value

You already met the idea in module 00 (lesson 6, A/B test). We repeat it with a different question: do two groups have the same mean?

• Null hypothesis H₀: "no difference" (the true means are equal).
• Test statistic: a number that measures how far the data stray from H₀. For two means: t = (x̄₁ − x̄₂) / √(s₁²/n₁ + s₂²/n₂) (Welch's test).
• p-value: the probability, if H₀ were true, of observing a statistic at least as extreme as the one obtained. It is not the probability that H₀ is true.

The three items are ordered from the simplest to the most misunderstood. The null hypothesis is a statement that can be rejected but not proved, the test statistic is a number computed from the data, and the p-value is the probability of that number (or of a more extreme one) computed assuming H₀ true: it does not answer the question "how likely is it that I am right".

rng = np.random.default_rng(5)
a = rng.normal(100, 15, size=50)                     # control group
b = rng.normal(106, 15, size=50)                     # treated group: true difference of 6 points

t_stat, p = stats.ttest_ind(b, a, equal_var=False)   # Welch's test
diff = b.mean() - a.mean()
se_diff = np.sqrt(a.var(ddof=1) / 50 + b.var(ddof=1) / 50)
d_cohen = diff / np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
print(f"observed difference: {diff:.2f} | standard error: {se_diff:.2f}")
print(f"t statistic = {t_stat:.3f} | p-value = {p:.4f}")
print(f"standardised effect (Cohen's d) = {d_cohen:.2f}")
print(f"95 % confidence interval of the difference: [{diff - 1.98 * se_diff:.2f}, {diff + 1.98 * se_diff:.2f}]")

Output obtained:

observed difference: 9.02 | standard error: 2.66
t statistic = 3.393 | p-value = 0.0010
standardised effect (Cohen's d) = 0.68
95 % confidence interval of the difference: [3.76, 14.28]

What the code does. It builds two groups with a true difference of 6 points (known because we inserted it), computes Welch's test (equal_var=False does not assume equal variances) and reports, besides the p-value, the difference with its standard error, the standardised effect d (the difference in units of standard deviation, independent of scale) and a confidence interval (with the multiplier 1.98 approximated for about 98 degrees of freedom). Note that the observed difference (about 9) is not exactly the true one (6): it is the variability of the sample, and the confidence interval, indeed, contains 6. This is the correct way to report a result: the p-value says whether the difference is distinguishable from noise, d and the interval say how large it is and how uncertain.

What a p-value really does, by simulation

rng = np.random.default_rng(6)
p_values = np.array([stats.ttest_ind(rng.normal(0, 1, 30), rng.normal(0, 1, 30)).pvalue for _ in range(5_000)])
print("H0 true, 5000 tests: fraction with p < 0.05 =", round(float((p_values < 0.05).mean()), 4))
print("deciles of the p-values (uniform if H0 is true):", np.round(np.percentile(p_values, np.arange(10, 100, 10)), 2))

Output obtained:

H0 true, 5000 tests: fraction with p < 0.05 = 0.0498
deciles of the p-values (uniform if H0 is true): [0.11 0.2  0.3  0.4  0.5  0.59 0.7  0.8  0.9 ]

What the code does. It repeats 5000 times a test on two groups drawn from the same population (H₀ true by construction). The result shows the most important property of p-values: if H₀ is true, they are uniformly distributed between 0 and 1, so 5% lie below 0.05 by pure chance. These are the false positives. If you run 20 tests on data without any effect, you expect about one "significant" (1 − 0.95²⁰ ≈ 64% of having at least one). The American Statistical Association's statement on p-values (2016) reminds us that the p-value measures neither the probability that the hypothesis is true nor the importance of the effect, and that conclusions should not be based on a single threshold.

Least-squares linear regression

The analogy: the rubber band

You have a cloud of points (xᵢ, yᵢ) and want the line ŷ = a + b·x that describes them best. Imagine connecting each point to the line with a vertical spring: the line settles where the overall tension is minimal. The tension of a spring grows with the square of the stretch, so the equilibrium line minimises the sum of the squares of the residuals rᵢ = yᵢ − (a + b·xᵢ). It is the method of least squares (Legendre 1805, Gauss 1809).

Derivation

We want to minimise S(a, b) = Σ (yᵢ − a − b·xᵢ)². A function of two variables has its minimum where the partial derivatives (the slope along a holding b still, and vice versa; we define them in the calculus lesson) are both zero.

With respect to a: −2 Σ (yᵢ − a − b·xᵢ) = 0, that is n·ȳ − n·a − b·n·x̄ = 0, hence a = ȳ − b·x̄: the line always passes through the point of the means (x̄, ȳ).

With respect to b: −2 Σ xᵢ (yᵢ − a − b·xᵢ) = 0. Substituting a = ȳ − b x̄ and rearranging: Σ xᵢ (yᵢ − ȳ) = b · Σ xᵢ (xᵢ − x̄). Since Σ x̄ (yᵢ − ȳ) = 0 and Σ x̄ (xᵢ − x̄) = 0, one can replace xᵢ with (xᵢ − x̄), and we arrive at:

b = Σ (xᵢ − x̄)(yᵢ − ȳ) / Σ (xᵢ − x̄)² = Cov(x, y) / Var(x).

The slope is the ratio between the covariance and the variance of x (the sample one, with the same divisor in numerator and denominator). One also defines the coefficient of determination R² = 1 − Σ rᵢ² / Σ (yᵢ − ȳ)²: the share of the variability of y that the line explains; for simple regression it coincides with the square of the correlation.

rng = np.random.default_rng(8)
x = rng.uniform(0, 10, size=80)
y = 2.0 + 0.5 * x + rng.normal(0, 1.0, size=80)         # true line: a = 2, b = 0.5, noise σ = 1

b = ((x - x.mean()) * (y - y.mean())).sum() / ((x - x.mean()) ** 2).sum()
a = y.mean() - b * x.mean()
residui = y - (a + b * x)
r2 = 1 - (residui ** 2).sum() / ((y - y.mean()) ** 2).sum()
print(f"a = {a:.3f}, b = {b:.3f}, R² = {r2:.3f}")
print("b = cov/var:", round(float(np.cov(x, y)[0, 1] / np.var(x, ddof=1)), 3), "| R² = corr²:", round(float(np.corrcoef(x, y)[0, 1] ** 2), 3))
print("np.polyfit check:", np.round(np.polyfit(x, y, 1), 3), "| scipy linregress:", round(float(stats.linregress(x, y).slope), 3))
print("sum of the residuals (≈ 0):", round(float(residui.sum()), 10))

# uncertainty on the slope with the bootstrap on the (x, y) pairs
pendenze = []
for _ in range(5_000):
    idx = rng.integers(0, x.size, x.size)
    xb, yb = x[idx], y[idx]
    pendenze.append(((xb - xb.mean()) * (yb - yb.mean())).sum() / ((xb - xb.mean()) ** 2).sum())
print("95 % bootstrap interval of the slope:", np.round(np.percentile(pendenze, [2.5, 97.5]), 3))

Output obtained:

a = 1.849, b = 0.524, R² = 0.664
b = cov/var: 0.524 | R² = corr²: 0.664
np.polyfit check: [0.524 1.849] | scipy linregress: 0.524
sum of the residuals (≈ 0): 0.0
95 % bootstrap interval of the slope: [0.45  0.604]

What the code does. It generates 80 points from a known line (a = 2, b = 0.5) plus noise, and applies the derived formulas. The estimates are close to the true values but not identical (it is the uncertainty of the sample). Then it verifies the same slope with np.cov/np.var (covariance over variance), R² as the square of the correlation, np.polyfit and scipy.stats.linregress. The sum of the residuals is zero (a consequence of the derivative with respect to a). The last block resamples the pairs (x, y) with replacement and recomputes the slope, obtaining an interval for the slope: if the procedure is correct, it contains it in about 95% of experiments (in this one it contains the true slope, 0.5; in one experiment out of twenty it would not).

Trade-off. Least squares is very sensitive to outliers (the square amplifies large deviations); it assumes a linear relationship, residuals with constant and independent dispersion. Always: plot the residuals (visualisation lesson: Anscombe's quartet has the same line and R² for four very different situations).

Correlation and causation

Two quantities can be correlated because one causes the other, because one causes the other the other way round, by chance, or because a third influences both (a confounder). We simulate the third case: temperature increases both the ice creams sold and the people who go swimming (and therefore, unfortunately, certain accidents).

rng = np.random.default_rng(9)
temp = rng.normal(25, 5, size=5_000)
gelati = 10 + 2.0 * temp + rng.normal(0, 4, 5_000)
bagnanti = -5 + 1.5 * temp + rng.normal(0, 4, 5_000)

print("corr(ice creams, swimmers) =", round(float(np.corrcoef(gelati, bagnanti)[0, 1]), 3))

def residuo(y, x):
    pend, interc = np.polyfit(x, y, 1)
    return y - (interc + pend * x)

print("corr after removing the effect of temperature =", round(float(np.corrcoef(residuo(gelati, temp), residuo(bagnanti, temp))[0, 1]), 3))

Output obtained:

corr(ice creams, swimmers) = 0.822
corr after removing the effect of temperature = -0.007

What the code does. It builds two variables that depend only on temperature (plus independent noise): there is no direct link between them, but the correlation is high. Then it computes the residuals of each with respect to temperature (what remains after removing the linear effect of temperature) and correlates the residuals: the correlation almost disappears. This is the logic of "controlling for a variable" and the reason why randomised experiments (lesson 6 of module 00) are the most solid method for talking about causes.

Typical mistakes (symptom, cause, remedy)

• Summarising with the mean alone. Remedy: median, IQR and a histogram.
• Divisor n for the variance of a sample. Remedy: ddof=1.
• Reading the confidence interval as a probability on the parameter. Remedy: "procedure that covers 95% of the time".
• Using 1.96 with small samples. Remedy: quantiles of the t.
• "p < 0.05 so the effect is large/true". Remedy: report effect and interval; know that with H₀ true 5% of tests are "significant".
• Multiple tests without correction. Remedy: a declared primary metric, or Bonferroni's correction (α/m).
• Bootstrap on dependent data (time series, rows of the same subject): resampling single rows breaks the structure. Remedy: block bootstrap or by group.
• Extrapolating the line outside the range of the data. Remedy: declare the domain of validity.
• Confusing correlation and causation. Remedy: confounders, experiments.

The mistakes fall into three families: those of description (the mean alone, the wrong divisor), those of reading (the 95%, the p-value, 1.96 on few data) and those of method (multiple tests, bootstrap on dependent data, extrapolation, causation). For the last two families no line of code protects you: you must declare before the calculation what you want to conclude and under which assumptions.

How to work in practice

Describe first, with a table of robust summaries and a chart. Always report estimate, uncertainty and sample size, and prefer intervals to p-values alone, declaring the test and the threshold in advance. Simulate to verify coverages and false positives, and document the assumptions (independence, approximate normality, variances), checking them on the data.

The vocabulary boils down to a few lines. Population and sample are the whole and the observed part; parameter, statistic and estimator are the unknown true number, the number computed from the data and the formula that produces it. The bias is the systematic mean error of an estimator, the standard error its standard deviation (for the mean, σ/√n) and the degrees of freedom are the values free to vary after the constraints imposed by the estimates. The likelihood is the probability of the data seen as a function of the parameter; the confidence interval comes from a procedure that covers the parameter with a declared frequency; the residual is the gap between observed and predicted; the confounder is a third variable that influences both quantities.

Mini-project

1. Simulate 10,000 samples from a Poisson with λ = 4 and n = 15. For each compute a 95% confidence interval for λ with x̄ ± 1.96·√(x̄/n) and measure the coverage. Change n to 5, 15, 100: how does it change?
2. With the bootstrap, estimate a 95% interval for the difference between the medians of the two groups of Welch's test and compare it with the one for the means.
3. Build an example in which the R² of the line is high but the residuals show a curvature, and plot them.

The first point is the check of the "coverage" criterion: the Poisson interval uses a normal approximation, and expecting it to work worse for small n is the part to measure, not to assume. The second reuses the bootstrap code with a different statistic, the third shows why a high R² is not enough to declare a relationship linear.

Challenge

• Easy. Compute with the formula and with np.polyfit the line for the points (1,2), (2,3), (3,5), (4,4) and verify that the sum of the residuals is zero.
• Medium. Prove that MSE = bias² + variance by adding and subtracting E[θ̂].
• Hard. Prove that for simple regression R² = r² (the square of the correlation coefficient).
• Expert. Derive the maximum likelihood estimator of λ for a Poisson (ℓ(λ) = Σ(xᵢ ln λ − λ − ln xᵢ!)) and prove that it is x̄. Is it unbiased? What is its standard error?
• Enterprise. Take a metric from one of your dashboards (conversion, response time). Write estimate, 95% bootstrap interval, and one sentence on what you cannot conclude with that sample.

The first four challenges are proofs and calculations that you can check with the code of the lesson; the fifth is the most formative because it forces you to write the concluding sentence with its limits. A statistical result is not complete until it also says what the sample does not allow you to claim.

What you can do now

You can summarise a sample with robust indicators, explain and verify by simulation why the sample variance divides by n − 1, estimate a parameter with maximum likelihood, build a confidence interval with Student's t or with the bootstrap and measure its coverage. You can read a p-value for what it is, derive the least-squares line and recognise a confounder. The next step is mathematical analysis: the derivative that here we used "by eye" to find the maxima of the likelihood and the minima of least squares; regression with many variables returns in Linear regression.

The same series also includes: 15 - scikit-learn: Estimators, Pipelines and Validation, 17 - Combinatorics: Counting Without Listing, 24 - Linear Algebra 2: Eigenvalues, SVD and PCA, 25 - Matrix Calculus and Automatic Differentiation: From the Jacobian to an Autograd Written From Scratch, 28 - Bayesian Statistics: Priors, Posteriors, Conjugates and Why Regularisation Is a Belief, 29 - Sampling and Monte Carlo Methods: Estimating What You Cannot Compute, 74 - Reinforcement Learning: The Agent That Learns by Trial and Error (and the Exploration-Exploitation Dilemma) (which refer back to this lesson).