Normal distribution and z-scores

#Why it matters

The normal distribution and the z score are the vocabulary of half the statistics you will use: confidence intervals, A/B test calculators, control charts, anomaly thresholds and standardised features in models all lean on them. The skill an employer wants is not the formula, it is judgement: knowing when data is close enough to normal for these tools to work, and what to do when it is not (most money and time data is not). Statistics appears in 58% of the Data Analyst and 82% of the Data Scientist job descriptions we read, and "what share of values lies within two standard deviations, and when is that wrong?" is a standard screening question.

#The concept

#The normal distribution

A normal distribution $N(\mu ,\sigma )$ is the symmetric bell curve set by two numbers: the mean $\mu$ (its centre) and the standard deviation $\sigma$ (its width). Mean, median and mode coincide. It shows up in nature when a quantity is the sum of many small independent effects (adult height, measurement error, the weight of animals of one species), and in statistics because averages of many observations tend towards it, which is the central limit theorem lesson that comes next.

The 68-95-99.7 rule: about 68.3% of values fall within 1σ of the mean, 95.4% within 2σ (exactly 95% within 1.96σ) and 99.7% within 3σ. It holds for normal data only.

#z-scores

A z-score says how many standard deviations a value sits from the mean:

Standardising puts different scales on one ruler: a z of 1.5 in exam marks and a z of 1.5 in salary are equally unusual, if both are normal. z-scores are also how most software turns a normal question into a standard-normal one, $N(0,1)$.

#Areas with scipy.stats.norm

Printed z tables (the grid of "area to the left of z" at the back of textbooks) are only for exams without laptops. In work, norm.cdf and norm.ppf give more digits and fewer lookup mistakes.

#Checking normality

1. Histogram with the fitted curve: quick, but bin choice can hide or invent shapes.

2. Q-Q plot (quantile-quantile): sort the data and plot it against the quantiles a normal distribution would have. Normal data lies on a straight line. A curve bending up at the right end means a long right tail; an S-shape means heavier tails on both sides. scipy.stats.probplot returns the points and the correlation of the line, which is a handy single number.

3. Formal tests such as Shapiro-Wilk: the null hypothesis is "the data is normal". With small samples the test has little power to detect anything; with large samples it rejects for deviations far too small to matter. At n in the thousands, almost all real data "fails". Use the Q-Q plot and ask whether the deviation matters for what you are about to do.

#Log-normal data

If the log of a variable is normal, the variable is log-normal: right-skewed, always positive, with a long tail. Incomes, order values, session lengths, house prices, taxi fares and file sizes are often roughly log-normal, because they come from multiplicative effects (a price times a quantity times a discount). Fit a normal to the logs, and convert back: the median of a log-normal is $e^{{\mu }_{\log }}$, the geometric mean.

#z-scores as outlier flags, and why they fail on skew

A common rule flags $|z|>3$. For normal data that flags about 0.27% of values. On right-skewed data it fails twice: the long tail inflates σ, so real oddities hide, and the threshold is symmetric, so the low side can never be flagged (mean minus 3σ is often negative). Use percentiles, a robust z based on median and MAD, or z-scores on the log scale.

#Worked example: penguin body mass

Palmer Penguins (CC0) again: body mass in grams for three species.

#Step 1: fit a normal per species and check it

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

url = "https://raw.githubusercontent.com/allisonhorst/palmerpenguins/main/inst/extdata/penguins.csv"
penguins = pd.read_csv(url).dropna(subset=["body_mass_g"])

rows = []
for species, grp in penguins.groupby("species"):
    mass = grp["body_mass_g"]
    (_, _), (_, _, qq_r) = stats.probplot(mass, dist="norm")
    rows.append({"species": species, "n": len(mass), "mean": mass.mean(), "sd": mass.std(),
                 "skew": stats.skew(mass), "qq_r": qq_r, "shapiro_p": stats.shapiro(mass).pvalue})
all_mass = penguins["body_mass_g"]
rows.append({"species": "all three", "n": len(all_mass), "mean": all_mass.mean(), "sd": all_mass.std(),
             "skew": stats.skew(all_mass), "qq_r": stats.probplot(all_mass)[1][2],
             "shapiro_p": stats.shapiro(all_mass).pvalue})
print(pd.DataFrame(rows).round(4).to_string(index=False))

  species   n      mean       sd   skew   qq_r  shapiro_p
   Adelie 151 3700.6623 458.5661 0.2825 0.9926     0.0324
Chinstrap  68 3733.0882 384.3351 0.2419 0.9899     0.5605
   Gentoo 123 5076.0163 504.1162 0.0688 0.9946     0.2336
all three 342 4201.7544 801.9545 0.4683 0.9806     0.0000

Within each species, body mass is close to normal: skew near zero and Q-Q correlations of 0.99. Shapiro-Wilk mildly rejects for Adelie (p = 0.03) on 151 birds, a small deviation you should look at, not panic about. Pooling the three species gives a clearly non-normal mix: two overlapping bells (Adelie and Chinstrap around 3.7 kg, Gentoo around 5.1 kg). Mixtures of groups are one of the most common reasons data "is not normal".

Draw the Q-Q plots to see it:

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 4, figsize=(15, 3.5))
groups = list(penguins.groupby("species")) + [("all three", penguins)]
for ax, (name, grp) in zip(axes, groups):
    stats.probplot(grp["body_mass_g"], dist="norm", plot=ax)
    ax.set_title(name)
fig.tight_layout()
fig.savefig("penguin_qq.png", dpi=120)

#Step 2: z-scores, areas and percentile ranks

A Gentoo penguin weighs 5,800 g. How unusual is that?

gentoo = penguins.loc[penguins["species"] == "Gentoo", "body_mass_g"]
mu, sigma = gentoo.mean(), gentoo.std()

z = (5800 - mu) / sigma
print(f"Gentoo: mean {mu:.0f} g, sd {sigma:.0f} g")
print(f"z for 5,800 g = {z:.2f}")
print(f"percentile rank from the normal model: {stats.norm.cdf(z):.1%}")
print(f"actual share of Gentoo at or below 5,800 g: {(gentoo <= 5800).mean():.1%}")
print(f"normal model: share between 4,500 and 5,500 g = "
      f"{stats.norm.cdf(5500, mu, sigma) - stats.norm.cdf(4500, mu, sigma):.1%}, "
      f"actual = {gentoo.between(4500, 5500).mean():.1%}")
print(f"heaviest 5% start at {stats.norm.ppf(0.95, mu, sigma):.0f} g (normal), "
      f"{gentoo.quantile(0.95):.0f} g (data)")
for k in (1, 2, 3):
    print(f"within {k} sd: {((gentoo - mu).abs() <= k * sigma).mean():.1%}  (normal: {stats.norm.cdf(k) - stats.norm.cdf(-k):.1%})")

Gentoo: mean 5076 g, sd 504 g
z for 5,800 g = 1.44
percentile rank from the normal model: 92.5%
actual share of Gentoo at or below 5,800 g: 92.7%
normal model: share between 4,500 and 5,500 g = 67.3%, actual = 64.2%
heaviest 5% start at 5905 g (normal), 5850 g (data)
within 1 sd: 67.5%  (normal: 68.3%)
within 2 sd: 98.4%  (normal: 95.4%)
within 3 sd: 100.0%  (normal: 99.7%)

When the data is close to normal, the model's answers match the counts within about three percentage points, and the 68-95-99.7 rule roughly holds (98.4% within 2 SD here, a sign of slightly light tails on 123 birds). That is the condition for using these tools.

#Worked example: taxi fares are not normal

Now the January 2019 NYC TLC yellow taxi fares from Mean vs median: central tendency and outliers, keeping only valid fares (above 0 and below 500).

#Step 3: normal against log-normal

import duckdb

fares = duckdb.sql("""
    select fare_amount from 'yellow_tripdata_2019-01.csv.gz'
    where tpep_pickup_datetime >= '2019-01-01' and tpep_pickup_datetime < '2019-02-01'
      and fare_amount > 0 and fare_amount < 500
""").df()["fare_amount"]

mu, sigma = fares.mean(), fares.std()
log_mu, log_sigma = np.log(fares).mean(), np.log(fares).std()
lognormal = stats.lognorm(s=log_sigma, scale=np.exp(log_mu))

print(f"n = {len(fares):,}, mean {mu:.2f}, sd {sigma:.2f}")
print(f"normal model: share of fares below 0 = {stats.norm.cdf(0, mu, sigma):.1%}")
table = pd.DataFrame({
    "data": fares.quantile([0.01, 0.10, 0.50, 0.90, 0.99]),
    "normal": stats.norm.ppf([0.01, 0.10, 0.50, 0.90, 0.99], mu, sigma),
    "lognormal": lognormal.ppf([0.01, 0.10, 0.50, 0.90, 0.99]),
})
print(table.round(2))
print(f"share above 50: data {(fares > 50).mean():.4f}, normal {stats.norm.sf(50, mu, sigma):.4f}, "
      f"lognormal {lognormal.sf(50):.4f}")

n = 7,657,663, mean 12.25, sd 11.13
normal model: share of fares below 0 = 13.6%
      data  normal  lognormal
0.01   3.5  -13.65       2.20
0.10   5.0   -2.02       4.27
0.50   8.8   12.25       9.66
0.90  23.5   26.52      21.83
0.99  52.0   38.15      42.45
share above 50: data 0.0298, normal 0.0003, lognormal 0.0049

The normal model predicts that 13.6% of taxi fares are negative, puts the 1st percentile at minus 13.65, and says fares above 50 almost never happen (0.03%), when in reality 3% of trips cost more than that. The log-normal is not perfect (it under-states the airport cluster around 52 and over-states the low end), but its median and spread are in the right place. For money data, start on the log scale.

#Step 4: Shapiro-Wilk at large n

rng = np.random.default_rng(3)
for n in [30, 300, 3000]:
    rejected = np.mean([stats.shapiro(rng.standard_t(10, size=n)).pvalue < 0.05 for _ in range(500)])
    print(f"t-distribution with 10 df, n = {n:>4}: Shapiro-Wilk rejects normality in {rejected:.0%} of samples")

t-distribution with 10 df, n =   30: Shapiro-Wilk rejects normality in 9% of samples
t-distribution with 10 df, n =  300: Shapiro-Wilk rejects normality in 49% of samples
t-distribution with 10 df, n = 3000: Shapiro-Wilk rejects normality in 100% of samples

A t-distribution with 10 degrees of freedom is very close to normal (slightly heavier tails), close enough for any t-test or confidence interval. Shapiro-Wilk barely notices at n = 30 and rejects every time at n = 3,000. The test answers "is this exactly normal?", which is never true of real data. The useful question is "is the deviation big enough to break the method I am using?" and you answer it with a Q-Q plot and knowledge of the method (t-tests on means, for example, are robust at large n).

#Step 5: z-score outlier flags on skewed data

z = (fares - mu) / sigma
print(f"|z| > 3 flags {(z.abs() > 3).mean():.2%} of fares; normal theory says 0.27%")
print(f"flagged on the low side: {(z < -3).sum()}, smallest flagged fare: {fares[z > 3].min():.2f}")

log_z = (np.log(fares) - log_mu) / log_sigma
print(f"log-scale |z| > 3 flags {(log_z.abs() > 3).mean():.2%}: "
      f"low side {(log_z < -3).sum():,} fares (max {fares[log_z < -3].max():.2f}), "
      f"high side {(log_z > 3).sum():,} fares (min {fares[log_z > 3].min():.2f})")

med = fares.median()
mad = stats.median_abs_deviation(fares, scale="normal")
robust_z = (fares - med) / mad
print(f"robust z (median, MAD) > 3 flags {(robust_z.abs() > 3).mean():.2%}")

|z| > 3 flags 3.27% of fares; normal theory says 0.27%
flagged on the low side: 0, smallest flagged fare: 45.66
log-scale |z| > 3 flags 0.32%: low side 1,060 fares (max 1.25), high side 23,629 fares (min 65.20)
robust z (median, MAD) > 3 flags 10.14%

On raw fares the |z| > 3 rule flags 3.3%, twelve times what it promises. Every flag is on the high side, and the cut-off of 45.66 sits just below the 52 JFK flat fare, so it labels airport trips as anomalies. The robust z (median and MAD) is worse here, flagging 10%: it measures distance in units of the tight middle of the data, so on skewed data it flags the whole tail. On the log scale the flags are rare (0.32%) and two-sided: 1,060 fares of 1.25 or less, which are almost certainly errors or cancelled trips, and 23,629 fares of 65.20 or more, which are long trips worth a look. No rule replaces looking at the rows, but on skewed data the log scale is the place to start.

#Common mistakes

• Quoting 68-95-99.7 for any data. It only holds for roughly normal data; for skewed delivery times you saw 84% within 1 SD in the variance lesson.

• Running Shapiro-Wilk on 10,000 rows and declaring the data unusable. At that size every test rejects. Look at the Q-Q plot and think about the method's robustness.

• Fitting one normal to a mixture of groups. Pooled penguin species are bimodal. Split by the obvious grouping first.

• Modelling money with a normal. It predicts negative prices. Use logs.

• Using |z| > 3 (or a robust z) on raw skewed data. It flags the profitable tail and never flags the low side. Move to the log scale first.

• Using z-score tables in a notebook. norm.cdf and norm.ppf are exact and fast.

#Interview angle

#Practice task

1. For each penguin species, fit a normal to flipper_length_mm and draw Q-Q plots. Which species looks least normal, and why might that be?

2. For Adelie penguins, use norm.cdf to estimate the share heavier than 4,500 g and compare with the actual share.

3. Using the Olist delivery days from Variance and standard deviation explained, compare the normal and log-normal predictions for the share of orders taking more than 30 days with the actual share.

4. Flag unusual Olist delivery times with |z| > 3 on raw days and on log days. Which flags would you send to the logistics team?

Done when each answer has the model's number next to the data's number and one sentence saying which model you trust.

for species, grp in penguins.dropna(subset=["flipper_length_mm"]).groupby("species"):
    f = grp["flipper_length_mm"]
    print(species, round(stats.probplot(f)[1][2], 4), round(stats.skew(f), 3))

adelie = penguins.loc[penguins["species"] == "Adelie", "body_mass_g"]
print("model:", round(stats.norm.sf(4500, adelie.mean(), adelie.std()), 3),
      "actual:", round((adelie > 4500).mean(), 3))

base = "https://raw.githubusercontent.com/olist/work-at-olist-data/master/datasets/"
o = pd.read_csv(base + "olist_orders_dataset.csv",
                parse_dates=["order_purchase_timestamp", "order_delivered_customer_date"])
days = ((o["order_delivered_customer_date"] - o["order_purchase_timestamp"]).dt.total_seconds() / 86400).dropna()
days = days[days > 0]
print("actual:", round((days > 30).mean(), 4),
      "normal:", round(stats.norm.sf(30, days.mean(), days.std()), 4),
      "lognormal:", round(stats.lognorm(s=np.log(days).std(), scale=np.exp(np.log(days).mean())).sf(30), 4))

#Further reading

• SciPy norm and lognorm: note that lognorm takes s as the log-scale SD and scale as exp of the log-scale mean.

• SciPy probplot: Q-Q plots and the fitted-line correlation.

• NIST/SEMATECH e-Handbook: normal probability plot: how to read the shapes.

• Limpert, Stahel and Abbt, Log-normal Distributions across the Sciences (BioScience, 2001): why so many real quantities are log-normal.