1.1 Data

Data is a collection of anything observable such as numbers, categories, words, images, sound, or other information that describes the world. Data becomes useful when it is organized and can be used for understanding the world, for making predictions or even for generating new data.

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from matplotlib.axes import Axes
from scipy import stats
from sklearn import datasets
from sklearn.cluster import KMeans
from statsmodels.regression.linear_model import RegressionResults

from statly import style

style.use()
rng = np.random.default_rng(1)
SPECIES = ["setosa", "versicolor", "virginica"]

# Real data that ship with statsmodels and scikit-learn (no download needed).
engel = (
    sm.datasets.engel.load_pandas().data
)  # income and food spending, Belgian households (1857)
iris = datasets.load_iris(as_frame=True).frame  # flowers of three iris species
visits = (
    sm.datasets.randhie.load_pandas().data
)  # doctor visits per person per year (RAND study)
sun = sm.datasets.sunspots.load_pandas().data  # sunspots per year, 1700 to 2008
co2 = sm.datasets.co2.load_pandas().data.interpolate()  # weekly CO2 at Mauna Loa, 1958 to 2001
cells = datasets.load_breast_cancer(
    as_frame=True
).frame  # cell nuclei measured in 569 breast tumour samples

1.1.1 How is data organised?

To be used, data is usually organised as a table. Each row is one thing we observed: a flower, a household, a week. Each column is one property we measured on it: a length, an income, a count. Statisticians call the rows observations and the columns variables.

Words, images and sound fit in a table too, once they are turned into numbers: a text becomes counts of its words, an image a grid of pixel values, a sound a series of measurements over time.

SPACE_INVADER = [
    "..X.....X..",
    "...X...X...",
    "..XXXXXXX..",
    ".XX.XXX.XX.",
    "XXXXXXXXXXX",
    "X.XXXXXXX.X",
    "X.X.....X.X",
    "...XX.XX...",
]


def plot_image_as_numbers(ax: Axes) -> None:
    """A space invader, 11 by 8 pixels, with each pixel's number written on it."""
    pixels = np.array([[char == "X" for char in line] for line in SPACE_INVADER], dtype=int)
    ax.imshow(pixels, cmap="Greys", vmin=0, vmax=1)
    for (row, col), value in np.ndenumerate(pixels):
        color = "white" if value == 1 else style.INK
        ax.text(col, row, f"{value:.0f}", ha="center", va="center", fontsize=7.5, color=color)
    ax.set(title="An image: a grid of numbers", xticks=[], yticks=[])
    ax.grid(False)


def plot_sound_as_numbers(ax: Axes, rate: int = 8000, milliseconds: float = 6) -> None:
    """A short tone, measured `rate` times a second: one number per measurement."""
    t = np.arange(0, milliseconds / 1000, 1 / rate)
    pressure = np.sin(2 * np.pi * 440 * t) + 0.4 * np.sin(2 * np.pi * 1320 * t)
    ax.plot(t * 1000, pressure, color=style.BLUE, lw=1.2)
    ax.plot(t * 1000, pressure, "o", ms=3.5, color=style.CORAL)
    first = ", ".join(f"{value:.2f}" for value in pressure[:5])
    style.note(ax, f"the numbers: {first}, …", y=0.04, va="bottom")
    ax.set_ylim(-2.0, 1.3)  # room at the bottom for the numbers
    ax.set(
        title="A sound: numbers over time", xlabel="time (milliseconds)", ylabel="air pressure"
    )


fig, (left, right) = plt.subplots(1, 2, figsize=(9.4, 3.4), width_ratios=[1, 1.9])
plot_image_as_numbers(left)
plot_sound_as_numbers(right)
fig.tight_layout()
plt.show()

Figure 1.1.1: An image and a sound are numbers too. Left: a space invader, 11 × 8 pixels, each one a number for its colour (0 is white, 1 is black). Right: a tone (simulated), measured 8,000 times a second; each dot is one number.

The 88 numbers of the image become one row of a table, with one column per pixel; a collection of images becomes a table with one row per image. A sound becomes a long column of measurements. This book mostly works with tables of numbers and categories.

def iris_sample(n: int = 6) -> pd.DataFrame:
    """A few rows of the iris table, with short column names and the species spelled out."""
    table = iris.sample(n, random_state=4).rename(
        columns=lambda name: name.replace(" (cm)", "")
    )
    table["species"] = table.pop("target").map(dict(enumerate(SPECIES)))
    return table


print(iris_sample().to_markdown(index=False))

Table 1.1.1: Six of the 150 iris flowers: one row per flower, one column per measurement (cm).

sepal length	sepal width	petal length	petal width	species
6.4	2.8	5.6	2.1	virginica
5.7	3.8	1.7	0.3	setosa
7.4	2.8	6.1	1.9	virginica
7.6	3	6.6	2.1	virginica
7.3	2.9	6.3	1.8	virginica
6	2.9	4.5	1.5	versicolor

The kind of values in a column decides what we can do with it:

def plot_measurement(ax: Axes) -> None:
    """Petal lengths: a measurement can take any value in a range."""
    style.hist(ax, iris["petal length (cm)"], bins=30)
    ax.set(title="Measurement", xlabel="petal length (cm): any value", ylabel="flowers")


def plot_count(ax: Axes, most: int = 10) -> None:
    """Doctor visits per year: a count is a whole number, never negative."""
    per_person = visits.mdvis.astype(int)
    people = per_person[per_person <= most].value_counts().sort_index()
    ax.bar(people.index, people.values, color=style.GOLD, alpha=0.8)
    ax.set(title="Count", xlabel="doctor visits: 0, 1, 2, …", ylabel="people")


def plot_category(ax: Axes) -> None:
    """Flowers per species: categories are groups without an order."""
    flowers = iris.target.map(dict(enumerate(SPECIES))).value_counts()
    ax.bar(
        flowers.index, flowers.values, color=[style.BLUE, style.GREEN, style.CORAL], alpha=0.8
    )
    ax.set(title="Category", xlabel="species: groups, no order", ylabel="flowers")
    ax.tick_params(axis="x", labelsize=8)


def plot_in_time(ax: Axes) -> None:
    """CO₂ through one year: values in time, where the order matters."""
    year = co2.co2.loc["2000"]
    ax.plot(year.index, year.values, "o-", ms=3, color=style.CORAL, lw=1)
    ax.set(title="Ordered in time", xlabel="through the year 2000", ylabel="CO₂ (ppm)")
    ax.xaxis.set_major_locator(mdates.MonthLocator(bymonth=[1, 7]))
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%b"))


fig, axes = style.subplots(4, cols=4, height=2.6, width=10.5)
for ax, plot in zip(axes, [plot_measurement, plot_count, plot_category, plot_in_time]):
    plot(ax)
fig.tight_layout()
plt.show()

Figure 1.1.2: Four kinds of variables. Each one asks for its own pictures, summaries and models.

• Measurements can take any value in a range (length, weight, income).
• Counts are whole numbers, never negative, often with many zeros.
• Categories are groups without an order (species, country). With two groups (yes/no, sick/healthy) they are called binary.
• Ordered in time (or space): the same thing measured again and again. The order carries information, and neighbouring values are related.

One more thing about every dataset: it is a sample. The 150 flowers stand in for all iris flowers, the 235 households for all households like them. Statistics is about what the sample tells us about the rest, and how sure we can be.

1.1.2 What can data be used for?

Almost every use of data is one of these:

• Describe: what is typical, how much things vary, what the shape is.
• Predict: a household earns 1,000. How much will it spend on food?
• Understand (statisticians say infer): how much does food spending rise with income, and is the relationship real, or could it be chance?
• Forecast: what comes next in a series, such as next month’s CO₂ or next year’s sunspots?
• Find structure: are there groups nobody told us about?
• Generate: make new data that behaves like the real thing, for simulations or for generative models that write text and draw images. This book simulates data all the time to show how statistics works.

def plot_describe(ax: Axes) -> None:
    """Household incomes: what is typical and how much they vary."""
    style.hist(ax, engel.income, bins=30)
    middle = engel.income.median()
    ax.axvline(middle, color=style.CORAL, lw=1.6)
    style.note(ax, f"typical: {middle:,.0f}", x=0.97, ha="right")
    ax.set(title="Describe", xlabel="household income", ylabel="households")


def food_line() -> RegressionResults:
    """A straight line through food spending against income."""
    return sm.OLS(engel.foodexp, sm.add_constant(engel.income)).fit()


def plot_food(ax: Axes, alpha: float) -> None:
    """Food spending against income, with the fitted line."""
    grid = np.linspace(0, 2500, 50)
    ax.scatter(engel.income, engel.foodexp, s=6, color=style.BLUE, alpha=alpha)
    ax.plot(grid, food_line().predict(sm.add_constant(grid)), color=style.GREY, lw=1.2)
    ax.set(xlabel="household income", ylabel="food spending", xlim=(0, 2500), ylim=(0, 1600))


def plot_predict(ax: Axes, income: float = 1000) -> None:
    """A new household: from its income, a guess at its food spending."""
    plot_food(ax, alpha=0.25)
    guess = food_line().predict([[1, income]])[0]
    ax.plot([income, income, 0], [0, guess, guess], color=style.CORAL, lw=1.2, ls="--")
    ax.plot(income, guess, "o", color=style.CORAL)
    style.note(ax, f"income {income:,.0f}\n→ food about {guess:,.0f}")
    ax.set_title("Predict")


def plot_understand(ax: Axes, step: float = 1000, x0: float = 1200) -> None:
    """How much food spending rises with income."""
    plot_food(ax, alpha=0.15)
    fit = food_line()
    rise = fit.params["income"] * step
    y0 = fit.params["const"] + fit.params["income"] * x0
    ax.plot([x0, x0 + step, x0 + step], [y0, y0, y0 + rise], color=style.CORAL, lw=1.5)
    style.note(ax, f"+{step:,.0f} income\n→ +{rise:,.0f} food")
    ax.set_title("Understand")


def plot_forecast(ax: Axes, last: str = "1998", years: int = 3) -> None:
    """CO₂ ahead: last year's pattern, raised by the average yearly increase."""
    weekly = co2.co2.resample("W").mean().interpolate()
    known = weekly[:last]
    one_year = 52
    rise = (known.iloc[-1] - known.iloc[-one_year * 10 - 1]) / 10
    pattern = known.iloc[-one_year:].to_numpy()
    ahead = np.concatenate([pattern + rise * (k + 1) for k in range(years)])
    dates = pd.date_range(known.index[-1], periods=len(ahead) + 1, freq="W")[1:]
    shown = known["1993":]
    ax.plot(shown.index, shown.values, color=style.BLUE, lw=1)
    ax.plot(dates, ahead, color=style.CORAL, lw=1, ls="--")
    style.note(ax, "forecast", x=0.97, y=0.1, ha="right", va="bottom")
    ax.set(title="Forecast", xlabel="year", ylabel="CO₂ (ppm)")
    ax.xaxis.set_major_locator(mdates.YearLocator(3))
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))


def plot_find_structure(ax: Axes, groups: int = 3) -> None:
    """Flowers sorted into groups from their petals alone, without being told the species."""
    petals = iris[["petal length (cm)", "petal width (cm)"]]
    found = KMeans(n_clusters=groups, n_init=10, random_state=0).fit_predict(petals)
    colors = np.array([style.BLUE, style.GREEN, style.GOLD])[found]
    ax.scatter(petals.iloc[:, 0], petals.iloc[:, 1], s=8, c=colors, alpha=0.7)
    style.note(ax, f"{groups} groups found", color=style.INK)
    ax.set(title="Find structure", xlabel="petal length (cm)", ylabel="petal width (cm)")


def plot_generate(ax: Axes, n: int = 235, seed: int = 2) -> None:
    """New, made-up households that behave like the real ones."""
    rng = np.random.default_rng(seed)
    fit = food_line()
    log_income = np.log(engel.income)
    income = np.exp(rng.normal(log_income.mean(), log_income.std(), n))
    noise = rng.normal(0, np.sqrt(fit.scale), n) * income / engel.income.mean()
    food = fit.params["const"] + fit.params["income"] * income + noise
    ax.scatter(engel.income, engel.foodexp, s=6, color=style.BLUE, alpha=0.3, label="real")
    ax.scatter(income, food, s=6, color=style.CORAL, alpha=0.5, label="made up")
    ax.legend(loc="upper left", fontsize=8, frameon=False)
    ax.set(
        title="Generate",
        xlabel="household income",
        ylabel="food spending",
        xlim=(0, 2500),
        ylim=(0, 1600),
    )


uses = [
    plot_describe,
    plot_predict,
    plot_understand,
    plot_forecast,
    plot_find_structure,
    plot_generate,
]
fig, axes = style.subplots(6, cols=3, height=2.8, width=10.5)
for ax, plot in zip(axes, uses):
    plot(ax)
fig.tight_layout()
plt.show()

Figure 1.1.3: Six uses of data. The made-up households (bottom right) are simulated; everything else is real data.

Later chapters take these one at a time.

1.1.3 Real data comes in shapes

Measure the same thing many times (the income of many households, the length of many flowers, the number of doctor visits of many people) and the values pile up in a pattern. That pattern is the data’s shape, and it is the first thing to look at.

The usual way to see a shape is a histogram. Put every value on a line, cut the line into ranges of the same width (called bins), and count how many values land in each range. Each count becomes a bar. Tall bars show where values are common, low bars where they are rare.

sepal = iris.loc[iris.target == 0, "sepal length (cm)"]
EDGES = np.arange(4.25, 6.0, 0.2)  # bin edges sit between values (measured to 0.1 cm)


def plot_dots(ax: Axes) -> None:
    """Every flower as one dot on the line (shifted up or down a little so none hide)."""
    ax.scatter(sepal, rng.uniform(-1, 1, len(sepal)), s=14, color=style.BLUE, alpha=0.7)
    ax.set(title="The values", ylim=(-4, 4))
    style.clean(ax, y=False)


def plot_stacks(ax: Axes) -> None:
    """The dots stacked in their bin, one above the other."""
    bin_of = np.digitize(sepal, EDGES)
    for b in np.unique(bin_of):
        middle = (EDGES[b - 1] + EDGES[b]) / 2
        n = int((bin_of == b).sum())
        ax.scatter([middle] * n, np.arange(1, n + 1), s=14, color=style.BLUE, alpha=0.7)
    for edge in EDGES:
        ax.axvline(edge, color=style.INK, lw=0.6, alpha=0.5)
    ax.set(title="Cut into bins, stacked", ylabel="flowers")


def plot_bars(ax: Axes) -> None:
    """One bar per bin, as tall as the number of flowers in it."""
    style.hist(ax, sepal, bins=EDGES)
    ax.set(title="The histogram", ylabel="flowers")


fig, axes = style.subplots(3, cols=3, height=2.6, width=10.5, sharex=True)
for ax, plot in zip(axes, [plot_dots, plot_stacks, plot_bars]):
    plot(ax)
    ax.set_xlabel("sepal length (cm)")
for ax in axes[1:]:
    ax.set_ylim(0, 13)
    ax.yaxis.set_major_locator(plt.MaxNLocator(integer=True))
fig.tight_layout()
plt.show()

Figure 1.1.4: How a histogram is made, from the sepal lengths of the 50 setosa flowers. Left: one dot per flower. Middle: the line cut into bins of 0.2 cm, with the dots stacked in their bin. Right: each stack becomes a bar as tall as its count.

The width of the bins is a choice. Very wide bins blur the shape into a few blocks; very narrow bins break it into noise. Somewhere in between, the shape shows.

def plot_visits(ax: Axes, most: int = 20) -> None:
    """People with 0, 1, 2, … doctor visits in a year (up to `most`; a few had more)."""
    per_person = visits.mdvis.astype(int)
    people = per_person[per_person <= most].value_counts().sort_index()
    ax.bar(people.index, people.values, color=style.GOLD, alpha=0.8)


def plot_co2(ax: Axes) -> None:
    """Weekly CO₂ at Mauna Loa, with a year label every ten years."""
    ax.plot(co2.index, co2.co2, color=style.CORAL, lw=1)
    ax.xaxis.set_major_locator(mdates.YearLocator(10))
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))


fig, ax = style.subplots(6, cols=3, height=2.8, width=10)
style.hist(ax[0], cells["mean smoothness"], bins=20, color=style.GREY)
style.hist(ax[1], engel.income, bins=40)
style.hist(ax[2], iris["petal length (cm)"], bins=30, color=style.GREEN)
plot_visits(ax[3])
style.mean_median(ax[0], cells["mean smoothness"], fmt="{:.3f}")
style.mean_median(ax[1], engel.income)
style.mean_median(ax[2], iris["petal length (cm)"], fmt="{:.1f}")
style.mean_median(ax[3], visits.mdvis, fmt="{:.1f}")
ax[4].plot(sun.YEAR, sun.SUNACTIVITY, color=style.PURPLE, lw=1.2)
plot_co2(ax[5])

labels = [  # (title, x axis, y axis) for each panel
    ("Cell smoothness: a bell", "smoothness of cell nuclei", "samples"),
    ("Income: long right tail", "household income", "households"),
    ("Petal length: two humps", "petal length (cm)", "flowers"),
    ("Doctor visits: counts, many 0s", "visits per year", "people"),
    ("Sunspots: cycles in time", "year", "sunspots"),
    ("CO₂: trend + season", "year", "CO₂ (ppm)"),
]
for a, (title, x, y) in zip(ax, labels):
    a.set(title=title, xlabel=x, ylabel=y)
fig.tight_layout()
plt.show()

Figure 1.1.5: Six real datasets, six different shapes.

The solid red line in each histogram is the mean: add all values and divide by how many there are. The dashed line is the median: the middle value, with half the values below it and half above. For a balanced shape they sit together. A long tail pulls the mean towards it, while the median stays with the typical value.

Each shape tells you something before any model does:

• A bell (cell smoothness): most values lie near the middle, fewer and fewer further out, about the same on both sides. The mean and the median sit together. This is the most common shape of all; the next section gives it its name, the normal distribution.
• A long right tail (income): most values are modest, a few are huge. The few pull the mean up, well above the median: the median is the more typical value.
• Two humps (petal length): there are probably two (or more) groups mixed together, here the species. The mean and the median both land in the gap: one average describes nobody.
• Whole numbers with many zeros (doctor visits): counts behave differently from measurements, and “zero” is often its own story.
• Ordered in time (sunspots, CO₂): the order matters. Today’s value depends on yesterday’s, and patterns like trends and seasons only exist because of the order.

1.1.4 Shapes have names: distributions

A distribution is a description of how values spread out: which values are common, which are rare. A few shapes turn up again and again, because they come from the same few mechanisms. When a named distribution fits the data well, a couple of numbers (its parameters) describe thousands of values.

def plot_normal_fit(ax: Axes) -> None:
    """Sepal length of one species, with a fitted normal distribution."""
    x = iris.loc[iris.target == 0, "sepal length (cm)"]
    width = style.hist(ax, x, bins=12)
    grid = np.linspace(x.min() - 0.4, x.max() + 0.4, 200)
    ax.plot(grid, stats.norm.pdf(grid, *stats.norm.fit(x)) * len(x) * width, color=style.CORAL)
    ax.set(title="Normal", xlabel="sepal length, one species (cm)", ylabel="flowers")


def plot_lognormal_fit(ax: Axes) -> None:
    """Household income, with a fitted log-normal distribution."""
    x = engel.income
    width = style.hist(ax, x, bins=40)
    grid = np.linspace(1, x.max(), 300)
    curve = stats.lognorm.pdf(grid, *stats.lognorm.fit(x, floc=0)) * len(x) * width
    ax.plot(grid, curve, color=style.CORAL)
    ax.set(title="Log-normal", xlabel="household income", ylabel="households")


def plot_count_fit(ax: Axes, most: int = 15) -> None:
    """Doctor visits, with a negative binomial: a distribution made for counts."""
    per_person = visits.mdvis.astype(int)
    people = per_person[per_person <= most].value_counts().sort_index()
    ax.bar(people.index, people.values, color=style.GOLD, alpha=0.7)
    k = np.arange(0, most + 1)
    mean, var = per_person.mean(), per_person.var()
    p = mean / var  # the negative binomial with the data's mean and variance
    n = mean * p / (1 - p)
    ax.plot(
        k, stats.nbinom.pmf(k, n, p) * len(per_person), "o-", ms=3, color=style.CORAL, lw=1.5
    )
    ax.set(title="Negative binomial", xlabel="doctor visits per year", ylabel="people")


def plot_mixture_fit(ax: Axes) -> None:
    """Petal length, with a mix of two normal distributions (two groups)."""
    x = iris["petal length (cm)"]
    width = style.hist(ax, x, bins=30, color=style.GREEN)
    grid = np.linspace(0.5, 7.5, 300)
    small, large = x[iris.target == 0], x[iris.target > 0]
    mix = sum(len(g) * stats.norm.pdf(grid, g.mean(), g.std()) for g in (small, large))
    ax.plot(grid, mix * width, color=style.CORAL)
    ax.set(title="Mixture of two", xlabel="petal length (cm)", ylabel="flowers")


fig, ax = style.subplots(4, cols=4, height=2.9, width=11)
for a, plot in zip(ax, [plot_normal_fit, plot_lognormal_fit, plot_count_fit, plot_mixture_fit]):
    plot(a)
fig.tight_layout()
plt.show()

Figure 1.1.6: The same kind of data with a named distribution drawn over it.

Each panel pairs a kind of data with a shape that suits it: measurements around a typical value (normal), positive values with a long tail (log-normal), counts (the negative binomial, a shape made for counts), and two groups mixed together (a mixture of two normals).

Try the shapes yourself. Each distribution below comes with the mechanism that produces it:

Which distribution is the right one?

No distribution fits real data exactly. The question is which one is close enough, and three checks settle most cases:

1. What kind of values are they? Only positive (incomes, waiting times)? Whole numbers (counts)? Between 0 and 1 (proportions)? A distribution that allows values the data can never take is a poor choice.
2. How did the values come about? Many small effects added together give a normal shape, multiplied together a log-normal, counted events a count distribution. The explorer above lists the mechanism for each one.
3. Draw it over the data. Fit the curve and look: does it follow the bars, tails included?

Household income shows all three at work:

def plot_candidate(ax: Axes, normal: bool) -> None:
    """Household income with one fitted curve over the bars: a normal or a log-normal."""
    x = engel.income
    width = style.hist(ax, x, bins=40)
    grid = np.linspace(-1200, x.max(), 600)
    if normal:
        params = stats.norm.fit(x)
        curve = stats.norm.pdf(grid, *params) * len(x) * width
        below_zero = grid < 0
        ax.fill_between(grid[below_zero], curve[below_zero], color=style.GREY, alpha=0.4)
        share = stats.norm.cdf(0, *params)
        note = f"expects {share:.0%} of households\nto earn less than 0"
        style.note(ax, note, x=0.45, color=style.INK)
    else:
        curve = stats.lognorm.pdf(grid, *stats.lognorm.fit(x, floc=0)) * len(x) * width
        style.note(ax, "never below 0,\nfollows the long tail", x=0.45)
    ax.plot(grid, curve, color=style.GREY if normal else style.CORAL, lw=2)
    ax.axvline(0, color=style.INK, lw=0.8)
    ax.set(
        title="Normal" if normal else "Log-normal",
        xlabel="household income",
        ylabel="households",
    )


fig, (left, right) = style.subplots(2, height=3.2)
plot_candidate(left, normal=True)
plot_candidate(right, normal=False)
fig.tight_layout()
plt.show()

Figure 1.1.7: Two candidates for household income. The normal curve (left) spills below zero and misses the long tail. The log-normal curve (right) stays above zero and follows the bars.

Income can never be negative (check 1), it comes from many effects that multiply, such as a raise being a percentage of what you already earn (check 2), and the log-normal follows the bars while the normal does not (check 3). Chapter 2.3 adds sharper tools for the third check: QQ plots, and tests that measure how far the data are from a curve. Chapters 2.1 and 2.2 go through each distribution in detail.

1.1.5 Why are distributions important?

Because a distribution turns how common into how likely.

The area under the curve is a probability

Take the flowers of one species and draw a histogram where each bar is the share of the flowers in its range. All the bars together hold 100% of the flowers. The fitted curve is a smooth version of the same bars, so the whole area under it is also 100%: every possible outcome.

Then any part of that area is the share of outcomes that land there. And the share of outcomes in a range is the probability that one new flower lands in that range.

BIN = 0.2  # bar width in cm: every height below is "share of flowers per 0.2 cm"


def setosa_sepals() -> pd.Series:
    """Sepal lengths of one species (setosa), in cm."""
    return iris.loc[iris.target == 0, "sepal length (cm)"]


def share_axis(ax: Axes) -> None:
    """Label the y axis as a share of the flowers, in percent."""
    ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda value, _: f"{value:.0%}"))
    ax.set(ylabel="share of flowers", xlabel="sepal length (cm)")


def plot_shares(ax: Axes) -> None:
    """Each bar is the share of flowers in its range; together the bars make 100%."""
    x = setosa_sepals()
    edges = np.arange(4.2, 6.01, BIN)
    ax.hist(x, bins=edges, weights=np.full(len(x), 1 / len(x)), color=style.BLUE, alpha=0.55)
    style.note(ax, "all bars together:\n100% of the flowers", x=0.6)
    ax.set_title("Bars: shares")
    share_axis(ax)


def curve(grid: np.ndarray) -> np.ndarray:
    """The fitted normal curve, scaled to the bars (share of flowers per bar width)."""
    return stats.norm.pdf(grid, *stats.norm.fit(setosa_sepals())) * BIN


def plot_whole_area(ax: Axes) -> None:
    """The fitted curve with all of its area shaded: 100% of outcomes."""
    grid = np.linspace(3.9, 6.1, 300)
    ax.fill_between(grid, curve(grid), color=style.BLUE, alpha=0.25)
    ax.plot(grid, curve(grid), color=style.BLUE)
    style.note(ax, "whole area:\n100% of outcomes", x=0.62, color=style.BLUE)
    ax.set_title("Curve: the whole area")
    share_axis(ax)


def plot_part_area(ax: Axes, cut: float = 5.5) -> None:
    """The area beyond `cut`: the probability that a flower is longer than that."""
    x = setosa_sepals()
    grid = np.linspace(3.9, 6.1, 300)
    beyond = grid >= cut
    ax.plot(grid, curve(grid), color=style.BLUE)
    ax.fill_between(grid[beyond], curve(grid[beyond]), color=style.CORAL, alpha=0.5)
    probability = stats.norm.sf(cut, *stats.norm.fit(x))
    note = (
        f"area beyond {cut} cm:\n{probability:.1%} of outcomes\n(data: {(x > cut).mean():.0%})"
    )
    style.note(ax, note, x=0.6)
    ax.set_title("Part of it: a probability")
    share_axis(ax)


fig, ax = style.subplots(3, height=3.1, width=10.5)
for a, plot in zip(ax, [plot_shares, plot_whole_area, plot_part_area]):
    plot(a)
fig.tight_layout()
plt.show()

Figure 1.1.8: From bars to probabilities. Left: each bar is a share of the flowers, and together they make 100%. Middle: the curve smooths the bars; its whole area is 100% of outcomes. Right: the part of the area beyond 5.5 cm is the probability that a flower is longer than 5.5 cm.

So the probability of a range is the area above it, not the height of the curve. The height only shows where outcomes are more or less common.

What follows from that

• Distributions compress. Two numbers, a log-normal’s centre and spread, describe all 235 incomes.
• They tell us what is unusual. A value far out in the tail sits over a tiny area: it is rare under that shape. That is the idea behind every statistical test.
• Models are built on them. A regression assumes how the noise is spread; counts call for a count distribution, yes/no outcomes for a binomial one.

def plot_unusual(ax: Axes, income: float = 4000) -> None:
    """How rare an income above `income` is under a fitted log-normal."""
    shape, loc, scale = stats.lognorm.fit(engel.income, floc=0)
    grid = np.linspace(1, 5200, 400)
    ax.plot(grid, stats.lognorm.pdf(grid, shape, loc, scale), color=style.BLUE)
    ax.axvline(income, color=style.CORAL, ls="--", lw=1.5)
    share = stats.lognorm.sf(income, shape, loc, scale)
    style.note(ax, f"above {income:,.0f}:\n{share:.2%} of\nhouseholds", x=0.3)
    ax.set(title="How unusual?", xlabel="household income", ylabel="how common")


def plot_chance(ax: Axes, flips: int = 20, heads: int = 15) -> None:
    """How often a fair coin gives `heads` or more heads in `flips` flips."""
    k = np.arange(0, flips + 1)
    colors = np.where(k >= heads, style.CORAL, style.BLUE)
    ax.bar(k, stats.binom.pmf(k, flips, 0.5), color=colors, alpha=0.7)
    share = stats.binom.sf(heads - 1, flips, 0.5)
    style.note(ax, f"{heads}+ heads in\n{flips} fair flips:\n{share:.1%}", x=0.66)
    ax.set(
        title="What chance looks like", xlabel=f"heads in {flips} flips", ylabel="probability"
    )
    ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda value, _: f"{value:.0%}"))


fig, ax = style.subplots(2, height=3.1, width=9)
plot_unusual(ax[0])
plot_chance(ax[1])
ax[0].set_yticks([])  # a curve's height only means something compared with itself
fig.tight_layout()
plt.show()

Figure 1.1.9: Two more things a distribution gives us: how unusual a value is (left), and what chance alone produces (right).

The right panel is a statistical test in miniature. If someone gets 15 heads in 20 flips, is the coin fair? A fair coin does that in about 2% of tries: a small area of the bars. That is rare, so we doubt the coin. Every test in this book works this way: work out what chance alone would produce, and see where the data fall.

1.1.6 The central limit theorem

Many small influences make a bell

Why is the bell the most common shape? Think of how tall a person grows. No single thing decides it. Hundreds of genes each add or take away a little, and so do food, sleep, illness and much more. A person’s height is the sum of all these small pushes, up and down.

Most people get a mix of ups and downs that largely cancel, so they end up near the middle. To end up very tall, nearly all the pushes must go up at once, and that is as rare as a long run of heads with a coin. So values near the middle are common, values far out are rare, and the further out, the rarer. That is the bell.

INFLUENCES, PEOPLE, START = 300, 10_000, 170
pushes = rng.uniform(-0.7, 0.7, (PEOPLE, INFLUENCES))
paths = START + pushes.cumsum(axis=1)  # each person's height after each influence
heights = paths[:, -1]


def plot_one_person(ax: Axes) -> None:
    """One person's height, built up one small push at a time."""
    ax.plot(paths[0], color=style.CORAL, lw=1.2)
    ax.axhline(START, color=style.INK, lw=0.8, ls="--")
    ax.set(title="One person", xlabel="influences added", ylabel="height (cm)")


def plot_many_people(ax: Axes, shown: int = 200) -> None:
    """Many people's paths: most stay near the middle."""
    ax.plot(paths[:shown].T, color=style.BLUE, lw=0.5, alpha=0.15)
    ax.plot(paths[0], color=style.CORAL, lw=1.2)
    ax.set(title=f"{shown} people", xlabel="influences added")


def plot_end_heights(ax: Axes) -> None:
    """Where everyone ends up: a bell, lying on its side."""
    width = style.hist(ax, heights, bins=40, orientation="horizontal")
    grid = np.linspace(heights.min(), heights.max(), 200)
    curve = stats.norm.pdf(grid, heights.mean(), heights.std()) * PEOPLE * width
    ax.plot(curve, grid, color=style.CORAL)
    ax.set(title="Final heights", xlabel="people")


fig, axes = style.subplots(3, cols=3, height=2.8, width=10.5, sharey=True)
for ax, plot in zip(axes, [plot_one_person, plot_many_people, plot_end_heights]):
    plot(ax)
fig.tight_layout()
plt.show()

Figure 1.1.10: Heights made of 300 small influences, each pushing up or down by up to 0.7 cm from 170 cm. Left: one person, influence by influence. Middle: 200 people; the paths spread out, but few reach the edges. Right: where 10,000 people end up, with a normal curve (red).

No single influence is normal: each one is equally likely to be anything between −0.7 and +0.7 cm. The bell comes from adding many of them. Many things we measure work the same way (blood pressure, the weight of a loaf of bread, errors in a measurement), and that is why the normal distribution is everywhere.

Averages make a bell too

An average is a sum divided by a count, so it gets the same bell. Draw a sample from any population, with replacement, and take its average. Repeat that thousands of times: the averages pile up in a bell around the population’s mean, and the larger each sample, the narrower and more normal the bell. This is the central limit theorem.

waiting_times = rng.exponential(1, 10_000)  # a population of waits, 1 minute on average
POPULATIONS = [
    ("Waiting times", waiting_times, "waiting time (minutes)"),
    ("Incomes", engel.income.to_numpy(), "household income"),
]


def plot_averages(
    ax: Axes, name: str, population: np.ndarray, n: int, samples: int = 20_000
) -> None:
    """Averages of `samples` samples of `n` values, with the normal curve the theorem predicts."""
    averages = rng.choice(population, (samples, n), replace=True).mean(axis=1)
    width = style.hist(ax, averages, bins=50, color=style.GREY if n == 1 else style.BLUE)
    if n > 1:
        mean, spread = population.mean(), population.std() / np.sqrt(n)
        grid = np.linspace(max(0, mean - 4 * spread), mean + 4 * spread, 300)
        ax.plot(grid, stats.norm.pdf(grid, mean, spread) * samples * width, color=style.CORAL)
    ax.set(
        title=f"{name}: " + ("single values" if n == 1 else f"averages of {n}"),
        ylabel="samples",
    )


fig, ax = style.subplots(6, cols=3, height=2.5, width=10)
for row, (name, population, xlabel) in enumerate(POPULATIONS):
    for a, n in zip(ax[3 * row : 3 * row + 3], (1, 5, 30)):
        plot_averages(a, name, population, n)
        a.set(xlabel=xlabel, xlim=(0, np.quantile(population, 0.99)))
fig.tight_layout()
plt.show()

Figure 1.1.11: Averages of samples drawn with replacement, from 10,000 simulated waiting times (top) and from the 235 household incomes (bottom). Single values are skewed. Averages of 5 are narrower but still lean, with a tail to the right. Averages of 30 are close to the normal curve (red).

Try it with other populations and sample sizes:

Why it matters:

• Averages are predictable even when single values are not. That is why we can put a confidence interval around an average, and why t-tests work on data that are not normal.
• The spread shrinks with √n. Four times as much data halves the uncertainty of an average.
• It has conditions. The values must be independent (not 30 measurements of the same person), and with very long tails, like incomes or city sizes, “large enough” can mean hundreds.

Without a distribution: the bootstrap

Not every question has a named distribution to lean on. What share of their income do the 235 households spend on food, and how sure can we be of that share? No simple bell answers that. The bootstrap uses the trick above, with the sample itself as the population: draw 235 households from the 235, with replacement (some come twice, some not at all), and work out the share again. Repeat that thousands of times. How much these shares vary shows how much the share could vary if we had drawn other households.

RESAMPLES = 10_000
food, income = engel.foodexp.to_numpy(), engel.income.to_numpy()
picks = rng.integers(0, len(income), (RESAMPLES, len(income)))  # households in each resample
shares = food[picks].sum(axis=1) / income[picks].sum(axis=1)
low, high = np.quantile(shares, [0.025, 0.975])
percent = plt.FuncFormatter(lambda value, _: f"{value:.0%}")


def plot_sample(ax: Axes) -> None:
    """The one sample we have: each household's food share, and the overall share."""
    style.hist(ax, food / income, bins=30, color=style.GREY)
    ax.axvline(food.sum() / income.sum(), color=style.INK, lw=1.6, ls="--")
    ax.xaxis.set_major_formatter(percent)
    ax.set(title="The sample", xlabel="share of income spent on food", ylabel="households")


def plot_shares(ax: Axes) -> None:
    """The overall share in each resample, with the middle 95% marked."""
    style.hist(ax, shares, bins=40)
    for edge in (low, high):
        ax.axvline(edge, color=style.CORAL, lw=1.6)
    ax.set_xlim(right=high + 0.8 * (high - low))  # room for the note
    style.note(ax, f"95% between\n{low:.1%} and {high:.1%}", x=0.98, y=0.6, ha="right")
    ax.xaxis.set_major_formatter(percent)
    ax.set(title="Resampled shares", xlabel="share of income spent on food", ylabel="resamples")


fig, axes = style.subplots(2, height=3.1, width=9)
plot_sample(axes[0])
plot_shares(axes[1])
fig.tight_layout()
plt.show()

Figure 1.1.12: The bootstrap. Left: the share of income each of the 235 households spends on food, and the share for all of them together (dashed). Right: that overall share in 10,000 resamples of the same 235 households, drawn with replacement. The middle 95% of them (red lines) show how far off the share could be.

No formula and no named shape: only the data, resampled. The price is that the sample must stand in well for the population, since the bootstrap can only reshuffle what the data holds.

1.1.7 In one picture

• Data is anything we observe, organised as a table: rows are observations, columns are variables (measurements, counts, categories, or values in time).
• Every dataset is a sample standing in for something bigger.
• We use data to describe, predict, understand, forecast, find structure and generate.
• Look at the shape first: tails, humps, counts, order in time.
• A distribution names a shape. The area under it is probability: all of it is 100%, a part of it is the chance of landing there.
• The central limit theorem: averages are normal, whatever the data. Most tests rest on this.
• The bootstrap needs no named distribution: resample the data with replacement and see how much the result varies.

Next: what a model is, and how flexible it should be.