The Statistics You Actually Need — Python Primer

You do not need a statistics degree to do useful data work. You do need about six ideas, and you need them properly rather than as remembered definitions.

Everything here is demonstrated by simulation. If you can generate a situation a thousand times and count what happens, you do not need the formula, and you will believe the answer in a way that a proof does not achieve.

The Only Idea That Matters

Your data is a sample. A different sample would have given a different answer.

Everything else in this chapter follows from that one sentence.

You measure the average order value in March and get £47.20. That number is not the truth about your customers. It is what this particular March produced. Another March, with the same underlying reality, would have given £46.80 or £48.10.

The whole point of statistics is quantifying how much that number would have wobbled, so you know whether a difference you observed is signal or shuffle.

Sampling Variation, Demonstrated

Rather than defining it, generate it.

import numpy as np
rng = np.random.default_rng(0)

# A world where the true mean is exactly 50
truth = 50

# Take 1,000 samples of 30 observations each
sample_means = [rng.normal(truth, 15, size=30).mean() for _ in range(1000)]

print(f"true mean:          {truth}")
print(f"average of samples: {np.mean(sample_means):.2f}")
print(f"spread of samples:  {np.std(sample_means):.2f}")
print(f"range seen:         {min(sample_means):.1f} to {max(sample_means):.1f}")

The samples centre on the truth, and they scatter around it substantially. Some individual samples are five or more away from 50, and each of those, taken alone, would have looked like a solid finding.

Now increase the sample size and watch the scatter shrink:

for n in [10, 30, 100, 1000]:
    means = [rng.normal(truth, 15, size=n).mean() for _ in range(1000)]
    print(f"n={n:>5}: spread of sample means = {np.std(means):.2f}")

The spread falls with the square root of the sample size. Four times the data halves the wobble, not quarters it. This is why the jump from 100 to 400 observations helps a lot and the jump from 10,000 to 40,000 helps far less.

The Bootstrap

Here is a technique that replaces a great deal of formula-based statistics with something you can reason about directly.

You have one sample. You cannot collect more. But you can resample your own data with replacement, which simulates the variation you would have seen.

data = rng.normal(50, 15, size=200)          # your actual sample

boot_means = [rng.choice(data, size=len(data), replace=True).mean()
              for _ in range(10_000)]

lo, hi = np.percentile(boot_means, [2.5, 97.5])
print(f"estimate: {data.mean():.2f}")
print(f"95% interval: {lo:.2f} to {hi:.2f}")

That interval is a confidence interval, and you built it from scratch without a formula.

The bootstrap works for almost any statistic, which is its real advantage. The median, the 90th percentile, a correlation, a ratio, the difference between two groups. All the same code with one line changed:

boot_medians = [np.median(rng.choice(data, len(data), replace=True))
                for _ in range(10_000)]

Try finding the standard formula for the confidence interval of a 90th percentile. The bootstrap does not care.

What a Confidence Interval Means

This is worth getting right, because almost everyone states it incorrectly.

A 95% confidence interval does not mean there is a 95% chance the true value is inside it. The true value is a fixed number; it is either in there or it is not.

It means: if you repeated this whole procedure many times, 95% of the intervals you construct would contain the truth.

The confidence is a property of the method, not of any particular interval. Demonstrate it:

truth = 50
contained = 0

for _ in range(1000):
    sample = rng.normal(truth, 15, size=30)
    boot = [rng.choice(sample, 30, replace=True).mean() for _ in range(500)]
    lo, hi = np.percentile(boot, [2.5, 97.5])
    contained += (lo <= truth <= hi)

print(f"{contained/10:.1f}% of intervals contained the truth")

You get roughly 95%. That is what the number means.

The practical value: report an interval, not a point. “£47.20” invites false precision. “£47.20, plausibly between £44 and £50” tells the reader how much to trust it, and stops a meeting arguing about a difference of thirty pence.

Hypothesis Tests and p-values

A p-value answers one narrow question: if there were genuinely no effect, how often would I see a difference at least this large by chance alone?

Simulate it rather than looking it up:

group_a = rng.normal(50, 15, 100)
group_b = rng.normal(52, 15, 100)
observed = group_b.mean() - group_a.mean()

# Shuffle the labels: what differences appear when there is NO real effect?
pooled = np.concatenate([group_a, group_b])
null = []
for _ in range(10_000):
    shuffled = rng.permutation(pooled)
    null.append(shuffled[100:].mean() - shuffled[:100].mean())

p = np.mean(np.abs(null) >= abs(observed))
print(f"observed difference: {observed:.2f}")
print(f"p-value: {p:.4f}")

That is a permutation test, and it is the clearest way to understand what a p-value is. You destroyed the group structure by shuffling, measured how big a difference chance produces, and counted how often chance beat your actual result.

What a p-value is not

Four things it is routinely mistaken for:

Not the probability the null hypothesis is true. It assumes the null and asks about the data.
Not the probability your finding is a fluke.
Not a measure of effect size. With a million rows, a difference of 0.01 gets p < 0.001 and means nothing commercially.
Not a threshold at which something becomes true. 0.049 and 0.051 are the same evidence.

The practical framing: a small p-value says “chance alone struggles to produce this”. It says nothing about whether the effect is large enough to care about, which is a business question the arithmetic cannot answer.

Effect Size

Always report the size of the difference alongside any test.

diff = group_b.mean() - group_a.mean()
pooled_sd = np.sqrt((group_a.var(ddof=1) + group_b.var(ddof=1)) / 2)
cohens_d = diff / pooled_sd

print(f"difference: {diff:.2f} units")
print(f"standardised: {cohens_d:.2f}")

Roughly, 0.2 is small, 0.5 is moderate, 0.8 is large. But the units matter more than the standardised number: a difference of £3 in average order value has a meaning your business already understands.

Multiple Comparisons

Test twenty things at the 5% level and you expect one false positive. That is not bad luck, it is arithmetic.

# Twenty groups, all drawn from the SAME distribution. No real effects.
false_positives = 0
for _ in range(1000):
    groups = [rng.normal(50, 15, 30) for _ in range(20)]
    baseline = groups[0]
    for g in groups[1:]:
        pooled = np.concatenate([baseline, g])
        null = [np.abs(rng.permutation(pooled)[30:].mean()
                     - rng.permutation(pooled)[:30].mean()) for _ in range(200)]
        if np.mean(null >= abs(g.mean() - baseline.mean())) < 0.05:
            false_positives += 1
            break

print(f"{false_positives/10:.0f}% of runs found at least one 'significant' "
      f"result where none exists")

You get roughly 60%. Nineteen comparisons, each with a one-in-twenty chance of a false alarm, and in most runs at least one fires.

This is what “we sliced the data and found that customers in the north-west aged 35 to 44 respond differently” usually is. If you searched, say you searched, and adjust:

from statsmodels.stats.multitest import multipletests
reject, adjusted, _, _ = multipletests(p_values, alpha=0.05, method='fdr_bh')

Or better, decide what you are testing before you look.

Correlation and Confounding

df[['age', 'income', 'revenue']].corr()

Correlation measures linear association only, on a scale from minus one to one. It misses curves entirely. Always plot the scatter; a correlation of zero can hide a perfect parabola.

The deeper trap is confounding, and Simpson’s paradox is its extreme form. A relationship can reverse when you account for a group.

# Two departments with different acceptance rates and different applicant mixes
dept_a = {'men': (400, 0.60), 'women': (100, 0.65)}   # (applicants, rate)
dept_b = {'men': (100, 0.20), 'women': (400, 0.25)}

for label, sex in [('men', 'men'), ('women', 'women')]:
    total = dept_a[sex][0] + dept_b[sex][0]
    admitted = dept_a[sex][0]*dept_a[sex][1] + dept_b[sex][0]*dept_b[sex][1]
    print(f"{label}: overall rate {admitted/total:.1%}")

Women have a higher acceptance rate in both departments and a lower rate overall, because more of them applied to the harder one. Both statements are arithmetically true. Which one you report is a choice, and it should be a considered one.

The lesson: before believing an aggregate, check whether it holds within the obvious subgroups.

Distributions

The normal distribution is less universal than introductions imply. Real business data is usually skewed: revenue, session length, income, order value all have a long right tail.

df.revenue.skew()        # 0 is symmetric; above 1 is strongly right-skewed

For skewed data, the median describes the typical case better than the mean, because a handful of large values drag the mean away from where most of the data sits.

print(f"mean:   {df.revenue.mean():.2f}")
print(f"median: {df.revenue.median():.2f}")

When those differ a lot, report both and say which one you are using. Someone will assume the other.

The reason the normal distribution appears so often is the central limit theorem: averages of many observations tend toward normal even when the individual values are not. That is about the sampling distribution of the mean, not about your raw data, and conflating the two causes real errors.

What Carries Forward

Your number would have been different with different data. Quantify how different.

Bootstrap it. Resampling gives you an interval for almost any statistic with five lines and no formula lookup.

Report intervals, not points.

A p-value is not an effect size. With enough rows, trivial differences become significant.

If you tested twenty things, say so. One will look real.

Check aggregates against subgroups before believing them.

Next: evaluation, which is the same scepticism applied to models.