Statistics for Data Science Part 1 — Descriptive Statistics & Distributions

You run one number over your data and paste it into the report: the average response time is 32 ms. It is arithmetically correct. It is also a lie that everyone in the room will believe.

Because your median response time is 20 ms, and 81% of your requests are faster than that “average.” The 32 comes from a thin tail of slow requests — timeouts, cold starts, a garbage-collection pause — that a single number folds in and hides. You reported the mean; your reader heard “typical”; those are not the same thing, and nobody in the meeting knows the difference.

import numpy as np
latency = ...                    # 10,000 real response times, mostly ~20ms, a few ~500ms
print(f"mean   = {latency.mean():.1f} ms")     # => mean   = 32.1 ms
print(f"median = {np.median(latency):.1f} ms") # => median = 20.2 ms
print(f"faster than the mean: {(latency < latency.mean()).mean():.0%}")  # => 81%

Statistics is not a branch of mathematics you were spared in school and can keep avoiding. It is the set of tools that stop you from fooling yourself with your own data — and for a programmer, every one of those tools is three lines of Python. This lesson teaches the concept, computes it in numpy/scipy/pandas, and then — the part textbooks skip — shows you when the number lies and what to reach for instead. Every statistic and every simulation below was executed on CPython 3.12 with a seeded random generator, so if you run the code you will get the numbers on the page.

---

Why this matters

Here is the uncomfortable truth about descriptive statistics: each summary number throws away almost all of your data on purpose, and the whole skill is knowing which one throws away the part that doesn’t matter. A mean crushes ten thousand numbers into one. That is the point — you cannot email someone ten thousand numbers — but the crushing is lossy, and different summaries lose different things. The mean loses the shape. The median loses the tail. The standard deviation assumes a shape that your data may not have.

A total beginner reaches for the mean for everything, because it is the “average” and school taught exactly one average. Then reality hands them skewed data — response times, incomes, file sizes, session lengths, anything with a floor of zero and no ceiling — and the mean quietly stops describing anything typical. The people it describes don’t exist: on the salary data later in this lesson, the “average” salary is one that nobody on the team earns and nobody is near.

The second thing nobody warns beginners about is that the tools disagree with each other by default. numpy.std() and pandas.std() give you different standard deviations for the same data — not a bug, a deliberate convention difference (ddof) that we will pin down exactly — and if you copy a threshold computed with one into code that uses the other, your alerts fire at the wrong level and you never notice. This lesson is as much about those silent disagreements as about the formulas.

And the third thing — the idea that makes the whole field work — is the Central Limit Theorem: no matter how weird and skewed your raw data is, the averages of samples drawn from it are normal. That single fact is why you can put error bars on a number, run an A/B test, or trust a poll of 1,000 people to predict millions. We will not assert it; we will draw twenty thousand samples from a lopsided distribution and watch the bell curve appear.

You need three libraries. Set up a virtual environment and install them — never pip install into the system Python:

python3 -m venv .venv
source .venv/bin/activate          # Windows: .venv\Scripts\activate
pip install numpy scipy pandas
python -c "import numpy, scipy, pandas; print(numpy.__version__, scipy.__version__, pandas.__version__)"
# => 2.5.1 1.18.0 3.0.3   (any recent versions are fine)

This lesson is the descriptive half. It builds directly on vectorised arrays — if .mean(axis=1) and broadcasting aren’t yet reflex, read NumPy: Arrays, Broadcasting & Vectorization first — and it sets up Part 2 (inference: confidence intervals, hypothesis tests, A/B testing), which is entirely built on the Central Limit Theorem you will prove at the end.

---

Central tendency: mean, median, mode

Central tendency is the question “where is the middle?” — and it has three different answers because “middle” means three different things.

The mean (arithmetic average) is the balance point: add everything up, divide by the count. Formally, for values x₁ … xₙ, the mean x̄ = (Σ xᵢ) / n. It uses every value, which is its strength and its weakness — every value includes the crazy ones.

The median is the middle rank: sort the data and take the value in the middle (or the average of the two middle values if n is even). Exactly half the data sits below it. It doesn’t care what the extreme values are, only that they’re extreme — moving the largest value from 100 to 100,000,000 doesn’t budge the median by a hair.

The mode is the most frequent value. It’s the only one of the three that works on categorical data (“which browser is most common?”) and the only one that can be undefined (all values unique) or multiple (a tie). On continuous data it’s nearly useless without binning — more on that trap later.

Here is all three, three ways — the stdlib statistics module, numpy, and pandas — on a deliberately tiny, hand-checkable dataset: the annual salaries (in ₹ lakh) of an eight-person engineering team, plus the founder.

import statistics as stats
import numpy as np
import pandas as pd

salaries = [6.5, 7.0, 7.2, 7.8, 8.0, 8.5, 9.0, 9.5, 48.0]   # 8 engineers + 1 founder

# --- stdlib statistics: pure Python, no dependency, scalar in / scalar out
print(stats.mean(salaries))     # => 12.38888888888889
print(stats.median(salaries))   # => 8.0
print(stats.mode(salaries))     # => 6.5   (first value seen; all are unique here)

# --- numpy: the workhorse for arrays
print(np.mean(salaries))        # => 12.38888888888889
print(np.median(salaries))      # => 8.0

# --- pandas: for a labelled Series/DataFrame column
s = pd.Series(salaries)
print(s.mean(), s.median())     # => 12.38888888888889 8.0
print(s.mode().tolist())        # => [6.5]   (Series.mode returns ALL modes)

The mean is 12.39. The median is 8.0. That is not a rounding difference — the mean is 55% higher than the median, and the reason is sitting at the end of the list. One founder earning 48 drags the average up past every single engineer. Count it: eight of the nine people earn below the “average” salary. The mean here describes a salary that literally nobody earns.

Three details that bite beginners, all visible above:

• statistics.mode used to raise on ties; since Python 3.8 it returns the first mode and never raises. pandas.Series.mode() instead returns every tied value as a Series — pd.Series([1,1,2,2,3]).mode().tolist() is [1, 2]. If your code assumes mode is one number, pandas will surprise you.
• numpy has no mode function. For a mode on discrete data use collections.Counter(data).most_common(1) or scipy.stats.mode.
• statistics functions take a plain iterable and return a scalar; numpy/pandas want an array/Series. Mixing them is fine, but statistics is pure Python and slow on large data — use it for teaching and tiny data, numpy for anything real.

---

When the mean lies: skew and outliers

The salary example was small enough to see the lie by eye. Real data isn’t. Let’s build a realistic dataset — 10,000 API response times — with a seeded generator so your numbers match mine exactly, then watch the mean and median diverge at scale.

Real latencies are right-skewed: there’s a hard floor (nothing is faster than instant), a dense body of normal requests, and a long thin tail of slow ones (timeouts, retries, cold starts). The log-normal distribution models that body perfectly; we bolt on a 2% tail of genuinely slow requests to be honest about production.

import numpy as np
import pandas as pd
from scipy import stats as sps

rng = np.random.default_rng(42)              # seed -> reproducible
base = rng.lognormal(mean=3.0, sigma=0.5, size=9800)   # ~20ms body, right-skewed
slow = rng.uniform(200, 800, size=200)                 # 2% slow requests (GC, timeouts)
latency = np.concatenate([base, slow])                 # 10,000 response times in ms

print(f"n      = {latency.size}")                # => n      = 10000
print(f"min    = {latency.min():.3f} ms")        # => min    = 2.238 ms
print(f"max    = {latency.max():.3f} ms")        # => max    = 799.680 ms
print(f"mean   = {latency.mean():.3f} ms")       # => mean   = 32.077 ms
print(f"median = {np.median(latency):.3f} ms")   # => median = 20.196 ms
print(f"mean / median = {latency.mean()/np.median(latency):.3f}")   # => 1.588
print(f"faster than mean = {(latency < latency.mean()).mean():.1%}")  # => 81.2%

The mean (32.08 ms) is 59% larger than the median (20.20 ms), and 81.2% of requests are faster than the mean. If you put “average latency: 32 ms” on a dashboard, you have overstated the typical experience by more than half, because 200 slow requests out of 10,000 — two percent of the data — moved the headline number. The median shrugged them off.

This is the single most important habit in descriptive statistics: the gap between the mean and the median is a skew detector you get for free.

When do you use which? The rule of thumb that will not let you down:

That last row is the professional’s move: for latency you rarely report the mean or the median — you report p95 and p99, because your users feel the slow tail, not the balance point. We’ll compute those next.

---

Dispersion: range, variance, standard deviation (and the ddof trap)

“The middle is 20 ms” tells you nothing about whether the data is all clustered at 20 or scattered from 2 to 800. Dispersion measures the spread — how far, typically, values sit from the centre.

The crude measure is the range (max − min). It uses exactly two values and is therefore hostage to both outliers; on our latency data the range is 799.680 − 2.238 = 797.442 ms, which tells you the tail exists but nothing about the body.

The real measures are variance and its square root, standard deviation. Variance is the average squared distance from the mean: σ² = Σ(xᵢ − x̄)² / n. We square the deviations so that above and below don’t cancel, and so that big deviations count disproportionately. The standard deviation σ = √variance un-squares it back into the original units (milliseconds, not milliseconds²), which is why we report it. Roughly: the standard deviation is the typical distance of a value from the mean.

print(f"range = {latency.max() - latency.min():.3f} ms")   # => range = 797.442 ms
print(f"variance = {latency.var():.3f}")                    # => variance = 5105.090
print(f"std dev  = {latency.std():.4f} ms")                 # => std dev  = 71.4499 ms

The ddof trap — why numpy and pandas disagree

Now the gotcha that silently corrupts real analyses. Compute the standard deviation of the same data four ways:

print(np.std(latency))              # => 71.44991343128834   (numpy default)
print(latency.std())                # => 71.44991343128834   (same thing)
print(pd.Series(latency).std())     # => 71.4534861949194    (pandas default — DIFFERENT!)
print(stats.stdev(latency))         # => 71.4534861949194    (statistics — matches pandas)
print(stats.pstdev(latency))        # => 71.44991343128834   (statistics — matches numpy)

numpy says 71.4499. pandas says 71.4535. Same numbers, different answer. This is not a bug and it is not floating-point noise — it is a deliberate convention difference called ddof (delta degrees of freedom), and it is one of the most common silent errors in data code.

The divisor in the variance formula is the switch:

• Population standard deviation divides by n (ddof=0). Use it when your data is the entire population — you have every value there is.
• Sample standard deviation divides by n − 1 (ddof=1). Use it when your data is a sample drawn from a larger population and you want to estimate that population’s spread. Dividing by n − 1 instead of n is Bessel’s correction, and it exists because a sample hugs its own mean a little too tightly — using n would systematically underestimate the true spread. The n − 1 nudges the estimate up to compensate.

The libraries pick opposite defaults:

You control it explicitly on either side — always do, in code that matters:

np.std(latency, ddof=0)              # => 71.44991...  population
np.std(latency, ddof=1)              # => 71.45348...  sample
pd.Series(latency).std(ddof=0)       # => 71.44991...  force population to match numpy

On 10,000 points the two answers differ by 0.005% — invisible, which is exactly why it’s dangerous: it’s never wrong enough to notice on big data, so you learn to ignore it, and then it burns you on small data. Watch what happens with n = 5:

small = [12.0, 15.0, 14.0, 10.0, 19.0]
print(np.std(small))            # => 3.0331501776206203   (ddof=0, ÷5)
print(np.std(small, ddof=1))    # => 3.391164991562634    (ddof=1, ÷4)
print(pd.Series(small).std())   # => 3.391164991562634    (pandas default = ddof=1)
print(np.std(small, ddof=1) / np.std(small))   # => 1.118033988749895
print(np.sqrt(5 / 4))                          # => 1.118033988749895

At n = 5 the sample std is 11.8% larger than the population std — and that ratio is exactly √(n / (n−1)) = √(5/4). If you compute a control-limit threshold in numpy (population) and your monitoring library recomputes it in pandas (sample), your thresholds are 12% apart on small windows. Nobody will ever find that bug by reading the code, because both lines say .std().

The rule: in real analysis you almost always have a sample, so ddof=1 is usually right — but the discipline that saves you is to pass ddof explicitly every time rather than trust a default you might misremember.

That last column previews the next section: when data is skewed or has outliers, the standard deviation is inflated by the same tail that inflates the mean. The IQR — built from quartiles — is the robust measure of spread, exactly as the median is the robust measure of centre.

---

Quantiles: quartiles, percentiles, IQR and the 5-number summary

A percentile is a threshold with a rank meaning: the 90th percentile is the value below which 90% of the data falls. Quartiles are just the 25th, 50th, and 75th percentiles — they cut the sorted data into four equal-sized quarters. The 50th percentile is the median you already know.

q1, q2, q3 = np.percentile(latency, [25, 50, 75])
print(f"Q1 (25th) = {q1:.3f} ms")     # => Q1 (25th) = 14.425 ms
print(f"Q2 (50th) = {q2:.3f} ms")     # => Q2 (50th) = 20.196 ms   (== median)
print(f"Q3 (75th) = {q3:.3f} ms")     # => Q3 (75th) = 28.648 ms

p90, p95, p99 = np.percentile(latency, [90, 95, 99])
print(f"p90 = {p90:.3f}   p95 = {p95:.3f}   p99 = {p99:.3f}")
# => p90 = 40.655   p95 = 51.347   p99 = 499.707

Look at the jump from p95 to p99: 51 ms to 500 ms. The slowest 1% of requests are ten times slower than the 95th percentile. That is the story the mean buried and the median ignored — and it’s why SRE dashboards live on p95/p99, not averages.

The interquartile range (IQR) is Q3 − Q1 — the width of the middle 50% of the data. It’s the robust cousin of the standard deviation: the outliers that inflate std to 71 ms live outside the quartiles, so the IQR doesn’t see them.

iqr = q3 - q1
print(f"IQR = {iqr:.3f} ms")          # => IQR = 14.223 ms
# std = 71.45 ms, but the middle-50% spread is only 14.22 ms — the std is nearly all TAIL

The 5-number summary and the boxplot

Five numbers describe any distribution’s location and spread without assuming a shape: min, Q1, median, Q3, max. This is exactly what a boxplot draws — the summary is the picture:

The box holds the middle 50%; the line in it is the median; the whiskers reach to the last point within 1.5 × IQR of the box, and anything beyond is drawn as an individual outlier dot. Which brings us to the rules for calling something an outlier.

Outlier rules: 1.5×IQR and the z-score

An outlier is a value far enough from the rest to be suspicious. “Far enough” needs a rule, and there are two standard ones.

The 1.5×IQR rule (Tukey’s, what boxplots use) sets fences at Q1 − 1.5×IQR and Q3 + 1.5×IQR; anything outside is an outlier. It’s built from quartiles, so it’s robust — the outliers don’t move the fences that catch them.

lo_fence = q1 - 1.5 * iqr
hi_fence = q3 + 1.5 * iqr
print(f"fences: [{lo_fence:.3f}, {hi_fence:.3f}]")   # => fences: [-6.908, 49.982]
outliers_iqr = latency[(latency < lo_fence) | (latency > hi_fence)]
print(f"IQR outliers: {outliers_iqr.size} ({outliers_iqr.size/latency.size:.2%})")
# => IQR outliers: 533 (5.33%)

The z-score rule measures each value’s distance from the mean in standard deviations: z = (x − mean) / std. A common cutoff is |z| > 3. But notice the flaw for skewed data — it’s built from the mean and std, the two measures the outliers inflate, so the outliers partly hide themselves:

z = (latency - latency.mean()) / latency.std(ddof=0)
print(f"max z = {z.max():.2f}")                        # => max z = 10.74
outliers_z = latency[np.abs(z) > 3]
print(f"z>3 outliers: {outliers_z.size}")              # => z>3 outliers: 177
# IQR flagged 533; z-score flagged only 177 — the inflated std raised the bar

The two rules disagree (533 vs 177) precisely because the data is skewed. The z-score’s own std is bloated by the tail, so its 3σ fence sits further out and catches fewer points. For skewed data, prefer the IQR rule — or the robust modified z-score built from the median and MAD (median absolute deviation), which flagged 417:

⚠️ “Outlier” does not mean “delete it.” Those 533 slow requests are real — they’re your production tail, the thing users complain about. Flagging is for investigation, not automatic removal. Delete an outlier only when you can name why it’s invalid (a sensor glitch, a test row, a −1 sentinel). Silently dropping real data to make a chart prettier is how analyses lie.

---

Shape: skewness and kurtosis

Central tendency locates the data, dispersion sizes it, and shape describes its silhouette — two numbers that quantify what a histogram would show you.

Skewness measures asymmetry — which way the tail leans:

• Skew ≈ 0 → symmetric (both tails equal).
• Skew > 0 → right-skewed / positive skew (long tail to the high side; mean > median). Latency, income.
• Skew < 0 → left-skewed / negative skew (long tail to the low side; mean < median). Exam scores near a ceiling.

Kurtosis measures tailedness — how much of the variance comes from rare extreme values (fat tails and a sharp peak vs. thin tails). scipy reports excess kurtosis, which subtracts 3 so that a normal distribution scores 0; positive means fatter tails than normal (more outliers), negative means thinner.

from scipy import stats as sps
print(f"skew          = {sps.skew(latency):.3f}")             # => skew          = 7.754
print(f"excess kurt   = {sps.kurtosis(latency):.3f}")         # => excess kurt   = 63.807
print(f"Pearson kurt  = {sps.kurtosis(latency, fisher=False):.3f}")  # => 66.807 (excess + 3)

# a symmetric normal sample for contrast:
normal_sample = rng.normal(50, 10, size=10000)
print(f"normal skew   = {sps.skew(normal_sample):.3f}")       # => normal skew   = -0.015
print(f"normal kurt   = {sps.kurtosis(normal_sample):.3f}")   # => normal kurt   = -0.038

The latency data has skew 7.75 (violently right-skewed — that tail) and excess kurtosis 63.8 (enormously fat-tailed — those slow requests are far, far out). The normal sample sits near 0 on both, as it must. That contrast is the interpretation: the further from zero, the less you should trust any tool that assumes normality (which includes the z-score rule, the empirical rule, and most of Part 2 applied to raw data).

Two gotchas:

• scipy’s default kurtosis is excess (Fisher, normal = 0); pandas’ .kurt() is also excess but bias-corrected — they’ll differ slightly (sps.kurtosis(latency) = 63.807 vs pd.Series(latency).kurt() = 63.839). If you expected “normal = 3,” pass fisher=False.
• Both need enough data to mean anything. Skew and kurtosis on n = 8 are noise; treat them as descriptive only past a few hundred points.

---

Probability distributions with scipy.stats

So far we’ve described data we already have. A probability distribution is a model for data we might get — a mathematical recipe that says which values are likely and which are rare. Naming the distribution your data follows lets you compute probabilities, simulate more of it, and (in Part 2) test hypotheses about it.

scipy.stats gives every distribution the same interface, which is the thing worth memorising once:

One subtlety that trips everyone: for continuous distributions, .pdf(x) is a density, not a probability — it can even exceed 1. The probability of any exact value is zero; only intervals have probability, which you get from .cdf. For discrete distributions, .pmf(k) genuinely is P(X = k).

The five you must know

Normal (Gaussian) — the bell curve. Models anything that’s the sum of many small independent effects: measurement error, heights, aggregate quantities, and — crucially — sample means (that’s the CLT, next section). Two parameters: the mean μ (loc, where it’s centred) and the standard deviation σ (scale, how wide).

mu, sigma = 100, 15                              # e.g. IQ scale
print(sps.norm.pdf(100, mu, sigma))              # => 0.02659615...  density at the peak
print(sps.norm.cdf(115, mu, sigma))              # => 0.84134475...  P(X ≤ 115) = 84.1%
print(sps.norm.cdf(130, mu, sigma))              # => 0.97724987...  P(X ≤ 130) = 97.7%
print(sps.norm.sf(130, mu, sigma))               # => 0.02275013...  P(X > 130)  = 2.3%
print(sps.norm.ppf(0.975, mu, sigma))            # => 129.39945...   the 97.5th percentile
print(sps.norm.ppf(0.5, mu, sigma))              # => 100.0          the median = mu

Those two constants every statistician has memorised come straight from ppf on the standard normal (μ=0, σ=1):

print(sps.norm.ppf(0.975))    # => 1.959963984540054   (~1.96 — the 95% CI multiplier)
print(sps.norm.ppf(0.95))     # => 1.6448536269514722  (~1.645 — one-sided 95%)

Frozen vs. unfrozen — the API gotcha. You can pass the parameters on every call (unfrozen), or bind them once into a frozen distribution object and call methods on it. Frozen is cleaner and less error-prone when you reuse the same distribution:

# unfrozen — parameters every time (easy to fat-finger)
sps.norm.cdf(115, 100, 15)          # => 0.8413447...
# frozen — bind once, reuse
rv = sps.norm(loc=100, scale=15)    # a frozen distribution
rv.cdf(115)                         # => 0.8413447...   same answer, no repeated params
rv.ppf(0.975)                       # => 129.39945...
rv.rvs(size=3, random_state=rng)    # draw from it

Uniform — every value in a range equally likely. Models “pick a random point in [a, b],” jitter, naive random sampling. Params: loc = start, scale = width.

print(sps.uniform.pdf(5, 0, 10))    # => 0.1    flat: 1/width everywhere in range
print(sps.uniform.cdf(2.5, 0, 10))  # => 0.25   a quarter of the way along
print(sps.uniform.mean(0, 10))      # => 5.0

Binomial — the count of successes in n independent yes/no trials, each with success probability p. Models “how many of 10 servers fail if each fails independently with p=5%?”, conversion counts, defect counts.

print(sps.binom.pmf(0, 10, 0.05))   # => 0.59873...  P(exactly 0 of 10 fail) = 59.9%
print(sps.binom.pmf(2, 10, 0.05))   # => 0.07463...  P(exactly 2 fail)       = 7.5%
print(sps.binom.sf(1, 10, 0.05))    # => 0.08613...  P(2 or more fail)       = 8.6%
print(sps.binom.mean(10, 0.05))     # => 0.5   = n·p    (expect half a failure)

Poisson — the count of events in a fixed interval when they happen at a constant average rate λ, independently. Models requests per second, arrivals per minute, defects per wafer. Its defining quirk: mean = variance = λ.

print(sps.poisson.pmf(0, 3))        # => 0.04978...  P(0 requests in a second) at λ=3
print(sps.poisson.pmf(3, 3))        # => 0.22404...  P(exactly 3) — the most likely count
print(sps.poisson.cdf(5, 3))        # => 0.91608...  P(≤ 5 requests) = 91.6%
print(sps.poisson.sf(5, 3))         # => 0.08391...  P(> 5) = 8.4%  (your capacity headroom)

Exponential — the waiting time between Poisson events. Where Poisson counts arrivals, exponential measures the gaps. Models time-between-requests, time-to-failure, inter-arrival times. Param: scale = the mean wait (= 1/λ). Its quirk: memoryless — how long you’ve already waited tells you nothing about how much longer.

print(sps.expon.cdf(20, scale=20))  # => 0.63212...  P(next request within 20ms) if mean gap=20ms
print(sps.expon.sf(50, scale=20))   # => 0.08208...  P(gap longer than 50ms) = 8.2%
print(sps.expon.mean(scale=20))     # => 20.0    mean = scale
print(sps.expon.ppf(0.5, scale=20)) # => 13.862... the MEDIAN wait < the MEAN wait (right-skew!)

Sampling: rng vs .rvs

There are two idioms for drawing random samples, and mixing up their seeding is a classic reproducibility bug.

# 1) numpy Generator methods — for the common distributions, fastest
rng = np.random.default_rng(42)
rng.normal(100, 15, size=5)         # normal draws
rng.poisson(3, size=5)              # poisson draws

# 2) scipy .rvs — for ANY scipy distribution; pass random_state to seed it
sps.poisson.rvs(3, size=5, random_state=rng)     # seed with the SAME rng for reproducibility
sps.expon.rvs(scale=20, size=5, random_state=np.random.default_rng(42))

⚠️ Seed with np.random.default_rng() (the modern Generator), not the legacy np.random.seed(). They are different generators producing different streams — np.random.seed(0) then np.random.randn(3) gives [1.7641, 0.4002, 0.9787], while np.random.default_rng(0).standard_normal(3) gives [0.1257, -0.1321, 0.6404]. The legacy global functions (np.random.rand, np.random.randn) still work but are discouraged: they share hidden global state, so one library reseeding breaks another’s reproducibility. Create an explicit rng and pass it around.

The empirical rule (68-95-99.7), verified

For a normal distribution, a fixed fraction of the data falls within each band of standard deviations from the mean — the empirical rule: ~68% within ±1σ, ~95% within ±2σ, ~99.7% within ±3σ. Don’t take it on faith; simulate a million normal values and count:

rng2 = np.random.default_rng(0)
sample = rng2.normal(100, 15, size=1_000_000)
m, sd = sample.mean(), sample.std()
for k in (1, 2, 3):
    within = np.mean(np.abs(sample - m) <= k * sd)
    print(f"within {k}σ: {within:.3%}")
# => within 1σ: 68.263%
# => within 2σ: 95.451%
# => within 3σ: 99.729%

Bang on the theoretical 68.27 / 95.45 / 99.73. But here’s the warning the rule comes with — it is only true for normal data. Run the identical count on our skewed latency:

for k in (1, 2, 3):
    within = np.mean(np.abs(latency - latency.mean()) <= k * latency.std())
    print(f"latency within {k}σ: {within:.2%}")
# => latency within 1σ: 97.96%   (not 68%!)
# => latency within 2σ: 98.00%
# => latency within 3σ: 98.23%

97.96% within 1σ instead of 68%. The empirical rule is wildly wrong here, because the tail inflates σ so much that ±1σ already swallows almost everything, while a ±2σ and ±3σ band barely add anything (they’d need to reach negative latency on the low side). Applying normal-distribution rules to non-normal data is one of the most common ways beginners produce confident nonsense.

---

Standardisation and z-scores

You can’t compare a value on one scale with a value on another — a “78” on a test scored out of 100 versus a “150” on a test scored out of 200 — until you put them on a common footing. Standardisation does that by converting each value to a z-score: how many standard deviations it sits above (or below) its own mean.

z = (x − μ) / σ. A z-score of 0 is exactly average; +1 is one std above; −2 is two std below. It’s unitless, so it compares across any scales.

# Candidate scored 78 on Test A (mean 70, sd 8) and 150 on Test B (mean 120, sd 20).
# Which is the better RELATIVE performance?
za = (78 - 70) / 8       # => 1.0   one std above A's mean
zb = (150 - 120) / 20    # => 1.5   one-and-a-half std above B's mean
print(za, zb)            # => 1.0 1.5

# Turn z-scores into percentiles via the standard normal cdf:
print(f"Test A: {sps.norm.cdf(za):.1%}")   # => Test A: 84.1%
print(f"Test B: {sps.norm.cdf(zb):.1%}")   # => Test B: 93.3%

The raw 150 looks worse than 78 as a fraction (0.75 vs 0.78), but standardised it’s the stronger result: the candidate is at the 93rd percentile on B versus the 84th on A. Standardisation is how you compare apples to oranges honestly — and it’s the required input for many ML models (which assume comparable feature scales) and for the outlier z-score rule you saw earlier.

To standardise a whole array, scipy.stats.zscore does it in one call (and confirms the defining property — standardised data has mean 0 and std 1):

arr = np.array([50, 60, 70, 80, 90])
print(sps.zscore(arr))          # => [-1.4142 -0.7071  0.  0.7071  1.4142]
print(sps.zscore(arr).mean())   # => 0.0            (always, by construction)
print(sps.zscore(arr).std())    # => 0.9999999999999999   (≈ 1, ddof=0 default)

A z-score maps straight to a percentile through the standard normal cdf — a handful are worth knowing by heart, because they turn “how many σ out” into “how unusual”:

⚠️ sps.zscore defaults to ddof=0 (population). If you’re standardising a sample and want n−1, pass ddof=1 — the same trap as before, in a new place.

---

The Central Limit Theorem — the idea everything rests on

Everything in Part 2 — confidence intervals, t-tests, A/B testing, “the poll has a margin of error of ±3%” — depends on one theorem. It sounds too good to be true, and it is provably true:

Take any population with a finite variance — normal, skewed, uniform, bimodal, anything. Draw samples of size n and compute each sample’s mean. As n grows, the distribution of those sample means approaches a normal distribution — regardless of the population’s shape — centred on the true mean, with standard deviation σ / √n.

Read that twice. It does not say your data becomes normal. It says the averages of samples become normal, even when the data is violently non-normal. That’s why you can put error bars on an average without knowing the population’s shape — and it’s why the normal distribution is the single most important one in statistics.

We’ll prove it by brute force. Take a strongly right-skewed population (an exponential, skew ≈ 2 — nothing like a bell), draw 20,000 samples at each of several sizes n, and watch the skew of the sample means march toward zero and their spread shrink like σ/√n.

rng = np.random.default_rng(7)
pop = rng.exponential(scale=20, size=1_000_000)     # the POPULATION — heavily skewed
print(f"population: mean={pop.mean():.3f}  median={np.median(pop):.3f}  skew={sps.skew(pop):.3f}")
# => population: mean=19.993  median=13.860  skew=1.992   (nothing like normal)

for n in (2, 5, 30, 100):
    means = rng.choice(pop, size=(20000, n)).mean(axis=1)   # 20,000 sample means
    print(f"n={n:>3}: std={means.std():6.3f}  (σ/√n={pop.std()/np.sqrt(n):6.3f})  "
          f"skew={sps.skew(means):+.3f}")
# => n=  2: std=14.093  (σ/√n=14.118)  skew=+1.376
# => n=  5: std= 8.979  (σ/√n= 8.929)  skew=+0.876
# => n= 30: std= 3.641  (σ/√n= 3.645)  skew=+0.331
# => n=100: std= 1.992  (σ/√n= 1.997)  skew=+0.205

Watch both columns move. The skew of the sample means falls from the population’s 1.99 to 1.38 (n=2) to 0.88 (n=5) to 0.33 (n=30) to 0.21 (n=100) — the lopsided exponential is becoming a symmetric bell. And the spread tracks σ/√n almost exactly (3.641 measured vs 3.645 predicted at n=30). The population never changed; only the sample size did.

Two facts fall out that Part 2 lives on:

• Symmetry emerges around n = 30, which is why “n ≥ 30” is the folk rule for when the CLT has “kicked in.” (It’s a rule of thumb — heavier skew needs a larger n. Our exponential is still visibly skewed at 30; a milder population would be bell-shaped sooner.)
• The standard error is σ/√n. To halve your uncertainty you must quadruple your sample — the √n is why bigger studies have diminishing returns.

It works from any shape. A perfectly flat uniform population (skew ≈ 0 already) also yields normal means — the theorem doesn’t need the source to be skewed, it just doesn’t care:

uni = rng.uniform(0, 1, size=1_000_000)          # flat population, skew ≈ 0
um = rng.choice(uni, size=(20000, 30)).mean(axis=1)
print(f"uniform pop skew={sps.skew(uni):+.4f} -> sample-mean skew={sps.skew(um):+.4f}")
# => uniform pop skew=-0.0002 -> sample-mean skew=-0.0217   (still normal-shaped means)

The diagram traces this whole arc — a lopsided population where the mean and median disagree, described with the 5-number summary, resampled twenty thousand times, and the sample means collapsing into the normal bell that makes 68-95-99.7 apply again:

The badges mark the six ideas to carry out: the mean lies on skewed data because the tail drags it above the median (1); describe with the robust 5-number summary before you model (2); resampling with one line of vectorised numpy is the experiment (3); the spread of the means shrinks as σ/√n (4); the bell always wins whatever the source shape (5); and that is precisely why the whole of inferential statistics — Part 2 — works at all (6).

---

Describing a dataset end-to-end

In practice you don’t compute these one at a time — you call df.describe() and read the whole story at a glance. It returns count, mean, std (ddof=1 — sample, remember!), min, the quartiles, and max in one shot:

df = pd.DataFrame({"latency_ms": latency})
print(df.describe())

         latency_ms
count  10000.000000
mean      32.077267
std       71.453486
min        2.237687
25%       14.425496
50%       20.196152
75%       28.648053
max      799.679666

Now read it critically — the whole point of this lesson in one glance:

• mean (32.1) ≫ median/50% (20.2) → right-skewed. Don’t quote the mean as “typical.”
• std (71.5) ≫ IQR-implied spread (28.6 − 14.4 = 14.2) → the standard deviation is almost all tail; the body is tight.
• max (799.7) is far past Q3 (28.6) → a long upper tail; there are outliers worth investigating.
• min (2.2) is close to Q1 (14.4) → the lower side is well-behaved; the skew is entirely on the high end.

You can widen the percentiles to expose the tail that the default quartiles hide — essential for latency:

print(df.describe(percentiles=[.5, .9, .95, .99]).round(3))
#  50%     20.196
#  90%     40.655
#  95%     51.347
#  99%    499.707     <- the p99 is 10x the p95: THAT is the SLA story

A note on correlation (setup for Part 2)

Describing one variable is where we stop today; describing how two move together is correlation, the bridge to Part 2. A quick preview of the two you’ll meet, because the choice between them is another silent trap:

Pearson’s r measures linear correlation (−1 to +1). Spearman’s ρ measures monotonic correlation — it ranks the data first, so it catches any consistent up-or-down relationship even if it’s curved, and it’s robust to outliers.

rng = np.random.default_rng(3)
x = rng.uniform(1, 100, 500)
y = x**2 + rng.normal(0, 500, 500)          # a real relationship, but CURVED
print(f"Pearson  = {sps.pearsonr(x, y).statistic:.4f}")    # => Pearson  = 0.9537
print(f"Spearman = {sps.spearmanr(x, y).statistic:.4f}")   # => Spearman = 0.9708
# Perfect monotonic curve y = x**3 exposes the gap:
xm = np.arange(1, 101); ym = xm**3
print(sps.pearsonr(xm, ym).statistic, sps.spearmanr(xm, ym).statistic)
# => 0.9176 1.0   -- Spearman sees the PERFECT monotonic link; Pearson underrates the curve

And a single outlier can invent or destroy a Pearson correlation while Spearman barely flinches:

xa = np.arange(1., 11.); ya = xa.copy()                    # perfectly correlated
print(sps.pearsonr(xa, ya).statistic)                      # => 1.0
xa2, ya2 = np.append(xa, 20.), np.append(ya, -50.)         # add ONE wild point
print(sps.pearsonr(xa2, ya2).statistic)                    # => -0.7308  (flipped to NEGATIVE!)
print(sps.spearmanr(xa2, ya2).statistic)                   # => 0.5      (far less fooled)

⚠️ And the rule to tattoo on the inside of your eyelids: correlation is not causation. A high r between two variables means they move together — not that one causes the other. Ice-cream sales and drowning deaths correlate strongly; neither causes the other (summer causes both). Correlation is a hint to investigate, never a conclusion. Part 2 makes this rigorous with hypothesis testing; for now, just never write “X drives Y” because you saw an r of 0.9.

---

Hands-on lab

Build a realistic skewed dataset and run the full descriptive workflow on it — mean-vs-median, the ddof gap, outliers two ways, the empirical rule, and a CLT simulation. Everything is seeded, so your numbers will match the comments exactly.

Setup (once):

python3 -m venv .venv && source .venv/bin/activate    # Windows: .venv\Scripts\activate
pip install numpy scipy pandas

Put the following in stats_lab.py and run it with python stats_lab.py.

Step 1 — build the dataset.

import numpy as np, pandas as pd
from scipy import stats as sps

rng = np.random.default_rng(42)
base = rng.lognormal(mean=3.0, sigma=0.5, size=9800)   # normal traffic, right-skewed
slow = rng.uniform(200, 800, size=200)                 # 2% slow requests
latency = np.concatenate([base, slow])
s = pd.Series(latency, name="latency_ms")
print("n =", s.size)                                   # => n = 10000

What just happened: a reproducible 10,000-point latency sample with a heavy right tail. (Statistics like mean/median/std are order-independent, so we don’t bother shuffling.)

Step 2 — the mean lies; the median doesn’t.

print(f"mean   = {s.mean():.3f} ms")                   # => mean   = 32.077 ms
print(f"median = {s.median():.3f} ms")                 # => median = 20.196 ms
print(f"mean/median = {s.mean()/s.median():.3f}")      # => mean/median = 1.588
print(f"faster than mean = {(s < s.mean()).mean():.1%}")  # => faster than mean = 81.2%

What just happened: the mean is 59% above the median and 81% of requests beat it — proof, on your own data, that “average latency” overstates the typical experience.

Step 3 — the ddof gap.

print(f"population std (ddof=0) = {s.std(ddof=0):.4f}")   # => 71.4499
print(f"sample std     (ddof=1) = {s.std(ddof=1):.4f}")   # => 71.4535
print(f"numpy default (ddof=0)  = {np.std(latency):.4f}") # => 71.4499
print(f"pandas default (ddof=1) = {s.std():.4f}")         # => 71.4535

What just happened: numpy and pandas gave different std by default. On n=10,000 it’s a 0.005% gap — but you saw it hit 11.8% at n=5. Pass ddof explicitly.

Step 4 — 5-number summary and IQR outliers.

q1, med, q3 = s.quantile([.25, .5, .75])
iqr = q3 - q1
lo, hi = q1 - 1.5*iqr, q3 + 1.5*iqr
print(f"Q1={q1:.3f}  med={med:.3f}  Q3={q3:.3f}  IQR={iqr:.3f}")
# => Q1=14.425  med=20.196  Q3=28.648  IQR=14.223
print(f"fences [{lo:.3f}, {hi:.3f}]")                  # => fences [-6.908, 49.982]
iqr_out = s[(s < lo) | (s > hi)]
print(f"IQR outliers = {iqr_out.size} ({iqr_out.size/s.size:.2%})")  # => 533 (5.33%)

What just happened: the robust IQR flagged 533 slow requests using only the quartiles — the outliers can’t hide the fence that catches them.

Step 5 — z-score outliers (and why they disagree).

z = (s - s.mean()) / s.std(ddof=0)
print(f"z>3 outliers = {(z.abs() > 3).sum()}   max z = {z.max():.2f}")
# => z>3 outliers = 177   max z = 10.74

What just happened: the z-score flagged only 177 vs the IQR’s 533 — the inflated std raised the bar. On skewed data, trust the IQR.

Step 6 — shape, and the empirical rule failing.

print(f"skew = {sps.skew(latency):.3f}   excess kurtosis = {sps.kurtosis(latency):.3f}")
# => skew = 7.754   excess kurtosis = 63.807
for k in (1, 2, 3):
    band = ((latency >= s.mean()-k*s.std(ddof=0)) & (latency <= s.mean()+k*s.std(ddof=0))).mean()
    print(f"within {k}σ = {band:.2%}  (normal says {['68','95','99.7'][k-1]}%)")
# => within 1σ = 97.96%  (normal says 68%)   <- rule BROKEN by skew

What just happened: skew 7.75 confirms the heavy tail, and the empirical rule gives 97.96% instead of 68% — a loud reminder not to apply normal-distribution rules to non-normal data.

Step 7 — describe it all at once.

print(s.describe().round(3).to_string())
# count    10000.000 / mean 32.077 / std 71.453 / min 2.238
# 25% 14.425 / 50% 20.196 / 75% 28.648 / max 799.680

What just happened: one call, the whole descriptive story — and you can now read it critically (mean ≫ median → skew; std ≫ IQR → tail-driven spread).

Step 8 — watch the Central Limit Theorem happen.

rng2 = np.random.default_rng(7)
pop = rng2.exponential(scale=20, size=1_000_000)       # skewed population
print(f"population skew = {sps.skew(pop):.3f}")        # => population skew = 1.992
for n in (2, 5, 30):
    means = rng2.choice(pop, size=(20000, n)).mean(axis=1)
    print(f"n={n:>2}: sample-mean skew={sps.skew(means):+.3f}  "
          f"std={means.std():.3f} (σ/√n={pop.std()/np.sqrt(n):.3f})")
# => n= 2: sample-mean skew=+1.376  std=14.093 (σ/√n=14.118)
# => n= 5: sample-mean skew=+0.876  std= 8.979 (σ/√n= 8.929)
# => n=30: sample-mean skew=+0.331  std= 3.641 (σ/√n= 3.645)

What just happened: you drew 20,000 samples at each size and watched the skew of their means fall from 1.99 toward 0 while the spread tracked σ/√n — the theorem that makes Part 2 possible, proven on your laptop in eight lines.

⚠️ Optional visual: if you also pip install matplotlib, add import matplotlib.pyplot as plt; plt.hist(means, bins=60); plt.show() after the n=30 line to see the bell. Plotting is covered in Matplotlib: Plotting Basics.

---

Common mistakes and troubleshooting

The three that cause the most silent damage, in prose:

The ddof default mismatch is the assassin, because it’s never loud. Both numpy and pandas spell it .std(); neither warns you they mean different things. On the big datasets where you first learn the tools, the difference is a rounding error, so you internalise “they’re the same.” Then you write a monitoring rule on small rolling windows (n=20), compute the control limits in numpy (population), and the alerting service recomputes them in pandas (sample) — now your ±3σ limits are ~2.5% apart and alerts fire slightly wrong forever. There is no traceback. The only defence is to make ddof explicit in every std/var call in code that matters, so the intent is on the page.

NaN poisoning vs NaN silently dropped — both are traps, in opposite directions. numpy reductions propagate NaN: one missing value turns np.mean of a million-element array into nan, which at least fails loudly. pandas, by contrast, silently skips NaN in .mean(), .std(), etc. — friendlier, but it means your “mean of 10,000 rows” might secretly be the mean of 9,300 because 700 were missing, and the denominator quietly changed under you. Neither behaviour is wrong; you just have to know which library you’re in. Count your non-nulls (s.count() vs len(s)) whenever missingness is possible — and handle the gaps deliberately (drop, fill, or interpolate) before you summarise, which is exactly what Pandas: GroupBy, Merge & Missing Data is for.

Applying normal-distribution tools to non-normal data is the confident-nonsense generator. The empirical rule, the |z|>3 outlier rule, “mean ± 2σ” error bars — every one assumes a bell shape, and every one is silently wrong on skewed data, as the 97.96%-in-1σ result showed. The habit that saves you: describe the shape before you trust a shape-dependent tool. One sps.skew() and a glance at mean-vs-median tells you in two seconds whether the normal-world rules even apply.

---

Cheat-sheet

---

Interview and exam questions

Q: When should you report the median instead of the mean, and how do you decide quickly? A: Report the median for skewed data or data with outliers — income, latency, file sizes, session lengths — because the mean is dragged toward the tail while the median is robust. The quick test is to compute both: if mean > median the data is right-skewed (a long high tail), if mean < median it’s left-skewed, and if they’re close it’s roughly symmetric and the mean is fine. On the latency example the mean (32.1 ms) was 59% above the median (20.2 ms) and 81% of requests were faster than the mean — so the median is the honest “typical” value.

Q: numpy.std([1,2,3,4,5]) and pandas.Series([1,2,3,4,5]).std() give different answers. Why? A: Different default ddof (delta degrees of freedom). numpy divides by n (ddof=0, the population std), pandas divides by n−1 (ddof=1, the sample std with Bessel’s correction). The sample version is larger by a factor of √(n/(n−1)) — for n=5 that’s √(5/4) ≈ 1.118, so pandas is ~11.8% higher. Fix by passing ddof explicitly. Use ddof=1 when your data is a sample from a larger population (almost always), ddof=0 when it is the whole population.

Q: What is Bessel’s correction and why n−1? A: Dividing the sample variance by n−1 instead of n. A sample’s own mean is, by construction, the point the sample sits closest to, so deviations measured from it are slightly too small — using n would systematically underestimate the true population variance. Dividing by the smaller n−1 inflates the estimate just enough to make it unbiased. The effect vanishes as n grows (which is why it only matters on small samples).

Q: Explain the 5-number summary and how it maps to a boxplot. A: Minimum, Q1 (25th percentile), median (50th), Q3 (75th), maximum. On a boxplot the box spans Q1–Q3 (its height is the IQR), the line inside is the median, the whiskers extend to the furthest points within 1.5×IQR of the box, and anything beyond is drawn as an individual outlier point. It describes location and spread without assuming any distribution shape, which is why it’s robust and why df.describe() leads with it.

Q: Two ways to detect outliers — when does each fail? A: The 1.5×IQR rule (fences at Q1−1.5·IQR and Q3+1.5·IQR) is built from quartiles, so it’s robust and works on skewed data — the default choice. The z-score rule (|z|>3) is built from the mean and std, which are themselves inflated by the outliers, so on skewed data it under-detects: on the latency set IQR flagged 533 points but z-score only 177, because the fat tail bloated σ and pushed the 3σ fence out. Use IQR (or the median/MAD “modified z-score”) on non-normal data; the z-score rule assumes roughly normal.

Q: State the Central Limit Theorem and why it matters. A: For any population with finite variance, the distribution of sample means approaches a normal distribution as the sample size n grows — regardless of the population’s shape — centred on the true mean with standard deviation σ/√n. It matters because it’s what lets us do inference: confidence intervals, t-tests, and A/B testing all rely on sample means being normal, and the CLT guarantees that without our needing to know or assume the population’s real distribution. It’s why the normal distribution is central to statistics even though real data rarely is normal.

Q: The CLT says averages become normal as n grows. Demonstrate you understand “how fast.” A: It depends on the population’s skew — the folk rule is n ≥ 30, but that’s only a rule of thumb. Simulating from a skewed exponential (population skew ≈ 2), the skew of the sample means fell to 1.38 at n=2, 0.88 at n=5, 0.33 at n=30, and 0.21 at n=100 — still not perfectly symmetric at 30 because the source is heavily skewed. A near-symmetric population (like a uniform) produces normal-looking means at much smaller n. Also, the spread of the means shrinks as σ/√n, so quadrupling the sample halves the uncertainty.

Q: What’s the difference between .pdf, .cdf, .ppf, and .sf on a scipy distribution? A: .pdf(x) (or .pmf(k) for discrete) is the density/probability at a point — for continuous distributions it’s a density, not a probability, and can exceed 1. .cdf(x) is P(X ≤ x), the cumulative probability to the left. .sf(x) is the survival function P(X > x) = 1 − cdf, more numerically accurate than 1 - cdf in the far tail. .ppf(q) is the inverse cdf — give it a probability, get the value at that percentile (e.g. norm.ppf(0.975) ≈ 1.96).

Q: Why standardise data, and what does a z-score of 1.5 mean? A: Standardising (z = (x−μ)/σ) converts values to a unitless “number of standard deviations from the mean,” so you can compare across different scales and feed scale-sensitive ML models. A z-score of 1.5 means the value is 1.5 standard deviations above its mean; via the standard normal, that’s roughly the 93rd percentile. It’s how you fairly compare a 78-on-100 (z=1.0, 84th percentile) with a 150-on-200 (z=1.5, 93rd) — the second is the stronger relative result despite the lower raw fraction.

Q: A colleague finds Pearson r = 0.2 between two variables and concludes they’re unrelated. What’s the flaw? A: Pearson only measures linear association. A strong but curved monotonic relationship (say y = x²) or a non-monotonic one can have low Pearson r while the variables are clearly related — check Spearman’s ρ (rank-based, catches any monotonic link) and, better, plot the data. Also, a single outlier can wreck or fabricate a Pearson r (one bad point flipped a perfect +1 to −0.73 in the lesson). Low Pearson r means “not linearly correlated,” not “unrelated.”

Q (coding): Write a function that summarises a numeric array robustly and flags outliers with the IQR rule. A:

import numpy as np
def robust_summary(x):
    x = np.asarray(x, dtype=float)
    x = x[~np.isnan(x)]                       # drop NaN explicitly
    q1, med, q3 = np.percentile(x, [25, 50, 75])
    iqr = q3 - q1
    lo, hi = q1 - 1.5*iqr, q3 + 1.5*iqr
    return {
        "n": x.size, "median": med, "iqr": iqr,
        "mean": x.mean(), "std_sample": x.std(ddof=1),
        "p95": np.percentile(x, 95), "p99": np.percentile(x, 99),
        "outliers": int(((x < lo) | (x > hi)).sum()),
    }

It leads with median/IQR (robust), includes the mean/std for reference, exposes p95/p99 for the tail, uses ddof=1 explicitly, and handles NaN — every habit from this lesson in one function.

Q: Why prefer np.random.default_rng() over np.random.seed()? A: default_rng() returns an explicit Generator object you pass around, so randomness is local and reproducible without hidden global state. np.random.seed() sets a global legacy generator shared by all np.random.* calls, so any library that reseeds breaks your reproducibility, and the two produce different streams (same seed, different numbers). default_rng is also statistically better (PCG64) and faster. Legacy functions still work, but all new code should use default_rng.

---

Key takeaways

• Every summary number throws away data on purpose — know which one loses the part that matters. The mean loses the shape (and lies on skewed data); the median loses the tail (and is robust); the std assumes a shape your data may not have.
• mean vs median is a free skew detector. If they diverge, the data is skewed — report the median (plus p95/p99), not the mean.
• ddof silently splits your tools: numpy defaults to population (÷n), pandas to sample (÷n−1). The gap is √(n/(n−1)) — invisible on big data, 12% at n=5. Pass ddof explicitly in code that matters.
• On skewed data, reach for robust tools: median over mean, IQR over std, 1.5×IQR over z-score. The z-score and empirical rules assume normality and fail loudly-but-silently when it’s absent (98% in ±1σ instead of 68%).
• scipy.stats gives every distribution one interface — .pdf/.pmf, .cdf, .sf, .ppf, .rvs — and lets you freeze parameters. Learn the five workhorses: normal, uniform, binomial, Poisson, exponential.
• Standardise (z-scores) to compare across scales and to feed scale-sensitive models; standardised data has mean 0 and std 1 by construction.
• The Central Limit Theorem is the keystone: the distribution of sample means becomes normal regardless of the population’s shape, with spread σ/√n. It’s why inference works — and the whole of Part 2 stands on it.
• Seed with np.random.default_rng(seed), never the legacy global np.random.seed — and remember correlation is never, by itself, causation.