Z-Score

A z-score, also called a standard score, says how many standard deviations a value sits from the mean. It turns raw numbers into a common scale so that measurements on different units can be compared.

\[ z = \frac{x - \mu}{\sigma} \qquad \text{or, from a sample,} \qquad z = \frac{x - \bar{x}}{s} \]

The reading is direct. \(z = 0\) is exactly at the mean, \(z = +1\) is one standard deviation above it, \(z = -2\) is two below. A student scoring 85 where the mean is 70 and the standard deviation is 10 has \(z = 1.5\). A person 150 cm tall where the mean is 170 and the standard deviation is 8 has \(z = -2.5\).

Two things are then usually said: that a large \(|z|\) marks an outlier, and that z-scores correspond to probabilities. Both are conditional, and the conditions are not small.

A z-score computed from the sample has a ceiling

When \(\mu\) and \(\sigma\) are known, \(z\) can be any size. When they are replaced by \(\bar{x}\) and \(s\) from the same \(n\) values, the extreme point is inside both of them: it pulls the mean toward itself and inflates the standard deviation it is divided by. The result is a hard bound that depends only on \(n\).

import numpy as np
from scipy.stats import norm, t, lognorm, expon, uniform

def max_z(n):
    x = np.zeros(n)
    x[0] = 1e12
    return float(abs((x - x.mean()) / x.std(ddof=1)).max())
"""One value pushed arbitrarily far from the rest. The z-score it earns
   stops growing, because the same value also drags the mean and the sd."""

print(f"{'n':>8}{'largest |z| reachable':>24}{'(n-1)/sqrt(n)':>16}")
for n in (5, 10, 11, 20, 30, 100, 1000):
    print(f"{n:>8}{max_z(n):>24.4f}{(n - 1) / np.sqrt(n):>16.4f}")

rng = np.random.default_rng(0)
m = 2_000_000
fams = (("normal", norm.rvs(size=m, random_state=1)),
        ("uniform", uniform.rvs(size=m, random_state=2)),
        ("t, 3 df", t.rvs(3, size=m, random_state=3)),
        ("exponential", expon.rvs(size=m, random_state=4)),
        ("lognormal, s=1", lognorm.rvs(1.0, size=m, random_state=5)))
print(f"{'distribution':>16}{'|z| > 2':>11}{'|z| > 3':>11}{'z > 3':>10}{'z < -3':>10}")
for name, x in fams:
    zz = (x - x.mean()) / x.std()
    print(f"{name:>16}{np.mean(np.abs(zz) > 2):>11.4f}{np.mean(np.abs(zz) > 3):>11.4f}"
          f"{np.mean(zz > 3):>10.4f}{np.mean(zz < -3):>10.4f}")

x = np.array([2, 3, 3, 4, 4, 4, 5, 5, 6, 40.0])
zz = (x - x.mean()) / x.std(ddof=1)
print(f"    data   {np.array2string(x, precision=1)}")
print(f"    mean {x.mean():.4f}, sd {x.std(ddof=1):.4f}, so the 40 scores z = {zz[-1]:.4f}")
med = np.median(x)
mad = np.median(np.abs(x - med))
print(f"    that is {zz[-1] / max_z(10) * 100:.1f}% of the largest z-score ten values can produce,")
print(f"    puts it at {0.6745 * (x[-1] - med) / mad:.4f}, which is the answer the eye gives")
#        n   largest |z| reachable   (n-1)/sqrt(n)
#        5                  1.7889          1.7889
#       10                  2.8460          2.8460
#       11                  3.0151          3.0151
#       20                  4.2485          4.2485
#       30                  5.2947          5.2947
#      100                  9.9000          9.9000
#     1000                 31.5912         31.5912
#     distribution    |z| > 2    |z| > 3     z > 3    z < -3
#           normal     0.0455     0.0027    0.0013    0.0014
#          uniform     0.0000     0.0000    0.0000    0.0000
#          t, 3 df     0.0412     0.0141    0.0070    0.0071
#      exponential     0.0497     0.0183    0.0183    0.0000
#   lognormal, s=1     0.0370     0.0180    0.0180    0.0000
#     data   [ 2.  3.  3.  4.  4.  4.  5.  5.  6. 40.]
#     mean 7.6000, sd 11.4426, so the 40 scores z = 2.8315
#     that is 99.5% of the largest z-score ten values can produce,
#     puts it at 24.2820, which is the answer the eye gives

What the two limits mean in practice

For z-scores computed from the same nonconstant sample, the maximum absolute value is \((n-1)/\sqrt{n}\) with ddof=1, and \(\sqrt{n-1}\) with ddof=0. The bound is attained when one value differs and all the others are equal. Therefore a strict |z| > 3 rule cannot flag an observation below n=11 under either convention; the exact bounds differ.

The standard repair is to compute the location and spread from statistics less sensitive to an isolated outlier: the median and the median absolute deviation. The median and MAD are robust, not immovable; if MAD is zero, the score needs a separate handling rule. The same 40 scores 24.28 that way, which is the answer anyone looking at the numbers would give. The alternative is to compute \(z\) with the suspect point held out, which costs \(n\) passes but keeps the familiar scale.

The second limit is the one people rarely state out loud. “Beyond two standard deviations” describes a distance, and it only becomes a rarity under the normal distribution. The table shows the share past \(|z| = 2\) running from 0.0000 for a uniform — which cannot reach 2 at all, since its own ceiling is 1.7321 — to 0.0497 for an exponential. Past \(|z| = 3\) the normal gives 0.0027 while a \(t\) with three degrees of freedom gives 0.0141, five times as many. And skew decides where they land: for the lognormal and the exponential, every value past 3 is above the mean and none is below.

What it is genuinely good for

Standardizing makes different scales commensurable, and that use survives everything above. Comparing a math mark with an English mark, or a latency in milliseconds with a payload in bytes, works because \(z\) removes the units — and it needs no distributional assumption, because it is only claiming a distance. Most machine-learning preprocessing uses it for exactly this: gradient-based methods converge badly when features differ in scale by orders of magnitude, and distance-based methods let the largest-scaled feature dominate.

Two habits keep it honest. When standardizing a training set, store \(\bar{x}\) and \(s\) and apply those same values to later data, rather than recomputing them — otherwise the transform shifts under the model. And when the sample standard deviation comes from few observations it is itself noisy, so a \(z\) built on it inherits that noise on top of everything else.

References

• SciPy documentation, scipy.stats.norm, scipy.stats.t and scipy.stats.lognorm. docs.scipy.org