1.2 Describing data · Visual Statistics

1.2.1 The centre

The mean adds all values and divides by how many there are:

$\bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i$

The bar over the x is the usual sign for “average”, and $\sum$ means “add up”. The mean is the point where the histogram would balance on a finger. The median is the middle value: half the values lie below it, half above. The mode is the most common value, the top of the histogram.

Code

income = engel.income.to_numpy()


def plot_balance(ax: Axes) -> None:
    """The histogram on a balance point at the mean, with the median for comparison."""
    style.hist(ax, income, bins=40)
    ax.plot(
        income.mean(),
        0,
        "^",
        ms=14,
        color=style.CORAL,
        clip_on=False,
        transform=ax.get_xaxis_transform(),
    )
    ax.axvline(np.median(income), color=style.INK, ls="--", lw=1.4)
    style.note(ax, f"mean {income.mean():,.0f}", x=0.55, y=0.9)
    style.note(ax, f"median {np.median(income):,.0f}", x=0.55, y=0.78, color=style.INK)
    ax.set(title="The balance point", xlabel="household income", ylabel="households")


def plot_one_value(ax: Axes) -> None:
    """Mean and median before and after making the largest income much larger."""
    changed = income.copy()
    changed[np.argmax(changed)] = 20_000
    rows = [("as measured", income), ("richest × 4", changed)]
    for i, (name, values) in enumerate(rows):
        ax.plot(
            values.mean(), i, "o", ms=9, color=style.CORAL, label="mean" if i == 0 else None
        )
        ax.plot(
            np.median(values), i, "s", ms=8, color=style.INK, label="median" if i == 0 else None
        )
    ax.set_yticks([0, 1], [name for name, _ in rows])
    ax.legend(fontsize=8, loc="lower right")
    ax.set(
        title="One extreme value", xlabel="household income", ylim=(-0.6, 1.6), xlim=(850, 1100)
    )


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_balance(axes[0])
plot_one_value(axes[1])
fig.tight_layout()
plt.show()

Figure 1.2.1: Left: the incomes of the 235 households balance at the mean (triangle); the long right tail pulls the balance point to the right of the median (dashed). Right: what one value does. Moving the richest household from 4,957 to 20,000 shifts the mean by 64 and leaves the median exactly where it was.

The mean uses every value, which makes it efficient and fragile: one wrong or extreme value moves it. The median uses only the order, which makes it robust. For a symmetric shape they agree; for a skewed one, report the median, or both.

The mean and the median as best guesses

There is a second way to find both. Pick a guess and add up how far every value is from it. If distances are counted as they are, the guess with the smallest total is the median. If they are squared first, the guess with the smallest total is the mean:

$\text{median} = \underset{c}{\arg\min} \sum_i |x_i - c|, \qquad \text{mean} = \underset{c}{\arg\min} \sum_i (x_i - c)^2$

( $\arg\min$ means “the c that makes the sum smallest”.)

Code

values = np.array([1.0, 1.6, 2.1, 2.4, 2.9, 3.3, 3.8, 4.5, 9.0])
guesses = np.linspace(0, 10, 500)


def plot_total(ax: Axes, power: int, centre: float, name: str, color: str) -> None:
    """Total distance (power 1) or squared distance (power 2) from each guess to the values."""
    total = (np.abs(values[:, None] - guesses) ** power).sum(axis=0)
    ax.plot(guesses, total, color=color, lw=2)
    ax.plot(values, np.zeros_like(values), "o", color=style.INK, ms=6, clip_on=False)
    ax.axvline(centre, color=color, ls="--", lw=1)
    style.note(ax, f"lowest at the {name}, {centre:.1f}", x=0.45, y=0.92, color=color)
    label = "total distance" if power == 1 else "total squared distance"
    ax.set(xlabel="guess c", ylabel=label, ylim=(0, None))


fig, axes = style.subplots(2, height=3.1, width=10)
plot_total(axes[0], 1, np.median(values), "median", style.BLUE)
plot_total(axes[1], 2, values.mean(), "mean", style.CORAL)
fig.tight_layout()
plt.show()

Figure 1.2.2: Nine values (dots) and every possible guess c across the bottom. Left: the total distance from c to all values, counted as is; it falls in straight pieces and is lowest at the median. Right: the total of the squared distances, a smooth bowl lowest at the mean. The large value on the right pulls the bowl’s bottom towards it, but not the median.

Squares punish large distances heavily, so the mean moves towards any value far away to shorten that one distance. This view returns in Chapter 4.1: a regression line is a mean drawn as a line, because it too makes the total squared distance smallest.

1.2.2 The spread: variance and standard deviation

Every value lies some distance from the mean. Some distances are positive, some negative, and they always add up to exactly zero, so their average is useless. Squaring makes them all positive. The variance is the average squared distance from the mean:

$s^2 = \frac{1}{n - 1} \sum_{i=1}^{n} (x_i - \bar{x})^2$

and the standard deviation is its square root, $s = \sqrt{s^2}$ , which brings it back to the units of the data.

Each squared distance is literally a square: a value 3 cm from the mean contributes a square of 3 × 3 cm². The variance is the size of the average square, and the standard deviation is the side of that average square.

Code

ten = iris.loc[iris.target == 0, "sepal length (cm)"].to_numpy()[:10]
mean = ten.mean()
sd = ten.std(ddof=1)


def plot_squares(ax: Axes, gap: float = 0.08) -> None:
    """Every deviation drawn as a square on the mean line, then the average square."""
    x0 = 0.0
    for value in ten:
        side = value - mean
        color = style.CORAL if side > 0 else style.BLUE
        ax.add_patch(
            plt.Rectangle((x0, min(mean, value)), abs(side), abs(side), color=color, alpha=0.3)
        )
        ax.plot(x0 + abs(side), value, "o", ms=5, color=color)
        x0 += abs(side) + gap
    ax.axhline(mean, color=style.INK, lw=1.2)
    x0 += 0.25
    ax.add_patch(plt.Rectangle((x0, mean), sd, sd, fill=False, ec=style.INK, lw=1.5, ls="--"))
    ax.text(
        x0 + sd / 2,
        mean + sd + 0.03,
        f"average square:\nside s = {sd:.2f} cm",
        ha="center",
        fontsize=9,
    )
    ax.text(0, mean - 0.05, f"mean {mean:.2f} cm", fontsize=9, va="top")
    ax.set(
        xlim=(-0.05, x0 + sd + 0.3), ylim=(mean - 0.55, mean + 0.65), aspect="equal", xticks=[]
    )
    ax.set_ylabel("sepal length (cm)")


fig, ax = plt.subplots(figsize=(8, 3.6))
plot_squares(ax)
plt.show()

Figure 1.2.3: Ten setosa flowers’ sepal lengths (dots) and their mean (line). Each value’s distance from the mean is the side of a square (red: above the mean, blue: below), drawn one after another. The variance is the area of the average square (dashed); the standard deviation is its side.

Why divide by n − 1 and not n? The distances are measured from the sample’s own mean, which sits in the middle of the sample by construction, so they come out a little too small. Dividing by n − 1 corrects that exactly on average. For large samples it hardly matters.

Code

samples = rng.normal(0, 1, (20_000, 5))


def plot_estimates(ax: Axes, ddof: int, title: str) -> None:
    """Histogram of variance estimates with their average and the truth."""
    estimates = samples.var(axis=1, ddof=ddof)
    style.hist(
        ax, estimates[estimates < 5], bins=50, color=style.GREY if ddof == 0 else style.BLUE
    )
    ax.axvline(1, color=style.INK, ls="--", lw=1.4)
    ax.axvline(estimates.mean(), color=style.CORAL, lw=1.6)
    style.note(ax, f"average estimate\n{estimates.mean():.2f}", x=0.55)
    ax.set(title=title, xlabel="estimated variance", ylabel="samples")


fig, axes = style.subplots(2, height=3, width=9.5, sharey=True)
plot_estimates(axes[0], 0, "Divide by n")
plot_estimates(axes[1], 1, "Divide by n − 1")
fig.tight_layout()
plt.show()

Figure 1.2.4: Variances of 20,000 samples of 5 values from a population with variance 1 (dashed). Dividing by n gives estimates that are too small on average (left); dividing by n − 1 is right on average (right).

1.2.3 What a standard deviation means

The standard deviation is roughly the typical distance of a value from the mean. Two groups can share a mean and differ completely in spread:

Code

def plot_spread(ax: Axes, values: np.ndarray, title: str, color: str) -> None:
    """A centred histogram with ± one standard deviation marked."""
    centred = values - values.mean()
    style.hist(ax, centred, bins=np.linspace(-1.6, 1.6, 33), color=color)
    s = centred.std(ddof=1)
    top = ax.get_ylim()[1]
    ax.annotate(
        "",
        xy=(-s, top * 0.92),
        xytext=(s, top * 0.92),
        arrowprops=dict(arrowstyle="<->", color=style.CORAL, lw=1.8),
    )
    within = (abs(centred) <= s).mean()
    style.note(ax, f"s = {s:.2f} cm\n{within:.0%} within ± s", x=0.03, color=style.INK)
    ax.set(title=title, xlabel="distance from the mean (cm)", ylabel="flowers")


fig, axes = style.subplots(2, height=3.1, width=9.5, sharex=True)
plot_spread(
    axes[0],
    iris.loc[iris.target == 0, "sepal width (cm)"].to_numpy(),
    "Setosa sepal width",
    style.BLUE,
)
plot_spread(
    axes[1],
    iris.loc[iris.target == 1, "petal length (cm)"].to_numpy(),
    "Versicolor petal length",
    style.GREEN,
)
fig.tight_layout()
plt.show()

Figure 1.2.5: Sepal width of setosa flowers and petal length of versicolor flowers, shifted to the same mean of 0 for comparison. Same centre, different spread: the red arrow is one standard deviation either side of the mean.

How many values lie within one or two standard deviations depends on the shape. For a bell, 68% and 95% (Chapter 2.1). For any shape at all, at least 75% lie within two standard deviations and at least 89% within three (Chebyshev’s inequality). So “more than three standard deviations from the mean” is unusual whatever the data look like.

To compare spreads across different units, divide by the mean: the coefficient of variation s / x̄. A spread of 10 kg means a lot for cats and little for cars.

1.2.4 Quantiles and the boxplot

The median splits the data in halves. Quartiles split it in quarters: 25% of the values lie below the first quartile, 75% below the third. The distance between them, the interquartile range (IQR), holds the middle half of the data. Like the median, it ignores the extremes.

The boxplot draws these five numbers at once:

Code

def plot_anatomy(ax: Axes) -> None:
    """A labelled boxplot of incomes."""
    ax.boxplot(
        income,
        vert=False,
        widths=0.5,
        patch_artist=True,
        boxprops=dict(facecolor=style.BLUE, alpha=0.3, color=style.INK),
        medianprops=dict(color=style.CORAL, lw=2),
        whiskerprops=dict(color=style.INK),
        capprops=dict(color=style.INK),
        flierprops=dict(markeredgecolor=style.INK),
    )
    q1, med, q3 = np.percentile(income, [25, 50, 75])
    labels = [
        (q1, q1 - 350, "right", "first quartile"),
        (med, med, "center", "median"),
        (q3, q3 + 350, "left", "third quartile"),
    ]
    for value, text_x, align, label in labels:
        ax.annotate(
            label,
            (value, 1.27),
            xytext=(text_x, 1.62 if align == "center" else 1.48),
            ha=align,
            fontsize=8.5,
            arrowprops=dict(arrowstyle="-", color=style.INK, lw=0.6),
        )
    top_whisker = income[income <= q3 + 1.5 * (q3 - q1)].max()
    ax.annotate(
        "whisker",
        (top_whisker, 1.02),
        xytext=(top_whisker, 0.55),
        ha="center",
        fontsize=8.5,
        arrowprops=dict(arrowstyle="-", color=style.INK, lw=0.6),
    )
    ax.annotate(
        "single values",
        (income.max(), 1.02),
        xytext=(income.max() - 400, 0.55),
        ha="center",
        fontsize=8.5,
        arrowprops=dict(arrowstyle="-", color=style.INK, lw=0.6),
    )
    ax.set(title="Anatomy", xlabel="household income", yticks=[], ylim=(0.3, 1.9))


def plot_hidden(ax_hist: Axes, ax_box: Axes) -> None:
    """A two-humped shape: clear in a histogram, hidden in a boxplot."""
    petal = iris["petal length (cm)"].to_numpy()
    style.hist(ax_hist, petal, bins=30, color=style.GREEN)
    ax_hist.set(title="Two humps", ylabel="flowers")
    ax_hist.tick_params(labelbottom=False)
    ax_box.boxplot(
        petal,
        vert=False,
        widths=0.6,
        patch_artist=True,
        boxprops=dict(facecolor=style.GREEN, alpha=0.3, color=style.INK),
        medianprops=dict(color=style.CORAL, lw=2),
        whiskerprops=dict(color=style.INK),
        capprops=dict(color=style.INK),
    )
    ax_box.set(xlabel="petal length (cm)", yticks=[])


fig = plt.figure(figsize=(10, 3.6))
grid = fig.add_gridspec(2, 2, height_ratios=[2, 1])
plot_anatomy(fig.add_subplot(grid[:, 0]))
top = fig.add_subplot(grid[0, 1])
plot_hidden(top, fig.add_subplot(grid[1, 1], sharex=top))
fig.tight_layout()
plt.show()

Figure 1.2.6: Left: the parts of a boxplot, on the household incomes. The box holds the middle half, the line in it is the median, the whiskers reach to the furthest values within 1.5 IQR of the box, and points beyond are drawn one by one. Right: boxplots compress, and can hide: the petal lengths of all iris flowers have two clear humps in the histogram (top), invisible in the boxplot (bottom).

Boxplots are best for comparing many groups side by side. For one group, a histogram shows more.

1.2.5 Two variables: covariance

Do two measurements move together? Put a cross at the two means. For every point, draw the rectangle from the cross to the point. Points up and to the right, or down and to the left, give rectangles with a positive area (both distances have the same sign); the other two corners give negative areas. The covariance is the average area:

$\text{cov}(x, y) = \frac{1}{n - 1} \sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})$

Code

def plot_rectangles(
    ax: Axes, x: np.ndarray, y: np.ndarray, title: str, xlabel: str, ylabel: str
) -> None:
    """Rectangles from the means to every point, coloured by the sign of their area."""
    mx, my = x.mean(), y.mean()
    for a, b in zip(x, y):
        color = style.CORAL if (a - mx) * (b - my) > 0 else style.BLUE
        ax.add_patch(
            plt.Rectangle(
                (min(mx, a), min(my, b)),
                abs(a - mx),
                abs(b - my),
                color=color,
                alpha=0.08,
                lw=0,
            )
        )
    ax.scatter(x, y, s=12, color=style.INK, alpha=0.8, zorder=3)
    ax.plot(mx, my, "+", ms=16, mew=2, color=style.INK)
    cov = np.cov(x, y)[0, 1]
    style.note(
        ax,
        f"covariance {cov:+.3f}" if abs(cov) < 1 else f"covariance {cov:+.1f}",
        x=0.03,
        y=0.98,
        color=style.INK,
    )
    ax.set(title=title, xlabel=xlabel, ylabel=ylabel)


versicolor = iris[iris.target == 1]
fig, axes = style.subplots(2, height=3.4, width=9.5)
plot_rectangles(
    axes[0],
    versicolor["petal length (cm)"].to_numpy(),
    versicolor["petal width (cm)"].to_numpy(),
    "Versicolor petals",
    "petal length (cm)",
    "petal width (cm)",
)
plot_rectangles(
    axes[1],
    states.poverty.to_numpy(),
    states.hs_grad.to_numpy(),
    "US states",
    "poverty (%)",
    "high school graduates (%)",
)
fig.tight_layout()
plt.show()

Figure 1.2.7: Covariance as rectangles from the means (cross) to each point; red rectangles count positive, blue negative. Left: petal length and width of versicolor flowers, mostly red: a positive covariance. Right: poverty and high school graduation across the US states, mostly blue: a negative covariance.

The covariance has a sign that means something, but its size depends on the units: measure the petals in millimetres and it grows 100 times. Divide by both standard deviations and the units cancel. That is the correlation:

$r = \frac{\text{cov}(x, y)}{s_x \, s_y}$

always between −1 and +1. It is the covariance of the two columns after turning both into z-scores. Chapter 3.6 shows how to test it and what it cannot see.

1.2.6 Many variables: the correlation matrix

With many columns, compute the correlation of every pair and show them as a grid of colours. Patterns jump out: blocks of measurements that move together, and ones that move against each other.

Code

names = [
    c.replace("_", " ").replace("od280/od315 of diluted wines", "od280")
    for c in wine.columns[:-1]
]
matrix = np.corrcoef(wine.iloc[:, :-1].to_numpy().T)


def plot_matrix(ax: Axes) -> None:
    """The correlation matrix as coloured cells with the values written in."""
    image = ax.imshow(matrix, cmap="RdBu_r", vmin=-1, vmax=1)
    for i in range(len(names)):
        for j in range(len(names)):
            ax.text(
                j,
                i,
                f"{matrix[i, j]:.1f}".replace("0.", ".").replace("-.", "−."),
                ha="center",
                va="center",
                fontsize=6.5,
                color="black" if abs(matrix[i, j]) < 0.6 else "white",
            )
    ax.set_xticks(range(len(names)), names, rotation=60, ha="right", fontsize=8)
    ax.set_yticks(range(len(names)), names, fontsize=8)
    ax.grid(False)
    plt.colorbar(image, ax=ax, shrink=0.8, label="correlation")


fig, ax = plt.subplots(figsize=(7.5, 6.2))
plot_matrix(ax)
plt.show()

Figure 1.2.8: Correlations between the 13 measurements of 178 wines. Red: move together, blue: move against each other. The phenol measurements form a block; colour intensity and hue go opposite ways.

The same table with covariances instead of correlations is the covariance matrix. It is the starting point of principal components (Chapter 7.1) and of the multivariate normal distribution (Chapter 2.4).

1.2.7 In one picture

What	Formula	Robust to extremes?
mean	$\bar x = \frac{1}{n}\sum x_i$	no
median	middle value	yes
variance	$s^2 = \frac{1}{n-1}\sum (x_i - \bar x)^2$	no
standard deviation	$s = \sqrt{s^2}$	no
interquartile range	third quartile − first quartile	yes
covariance	$\frac{1}{n-1}\sum (x_i - \bar x)(y_i - \bar y)$	no
correlation	$r = \text{cov} / (s_x s_y)$	no (Spearman’s is)

The mean is the balance point; the median the middle. Skewed data pull them apart.
The standard deviation is the side of the average square: the typical distance from the mean.
Whatever the shape, at least 75% of values lie within two standard deviations.
Boxplots compare groups; histograms show shape.
Covariance is the average rectangle; correlation is covariance without units.

Next: what a model is.