Counts are whole numbers, never negative, and they come with their own shapes. How many tries, how rare the event, and how much the units differ decide which one.
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes
from scipy import stats
from statly import style
style.use()
rng = np.random.default_rng(1)
visits = (
sm.datasets.randhie.load_pandas().data
) # doctor visits per person per year (RAND study)A coin toss, a patient who recovers or not, a customer who buys or not: one try, two outcomes. Repeat the try a fixed number of times and count the yeses. That count follows a binomial distribution, with two parameters: the number of tries and the chance of a yes in each.
A Galton board shows it. A ball falls through rows of pins and at each pin goes left or right with equal chance. Where it lands is the number of times it went right.
ROWS, BALLS = 12, 10_000
rights = rng.integers(0, 2, (BALLS, ROWS)) # 1 = the ball goes right at that pin
landed = rights.sum(axis=1)
def plot_board(ax: Axes, shown: int = 40) -> None:
"""The paths of a few balls: each pin sends them left or right."""
position = np.hstack([np.zeros((shown, 1)), (2 * rights[:shown] - 1).cumsum(axis=1)])
for path in position:
ax.plot(path / 2 + ROWS / 2, np.arange(ROWS + 1), color=style.BLUE, lw=0.8, alpha=0.35)
ax.invert_yaxis()
ax.set(title="The board", xlabel="where the ball lands (rights)", ylabel="row of pins")
def plot_landing(ax: Axes) -> None:
"""Where the balls land, against the binomial distribution."""
k = np.arange(ROWS + 1)
counts = np.bincount(landed, minlength=ROWS + 1)
ax.bar(k, counts, color=style.BLUE, alpha=0.55)
ax.plot(k, stats.binom.pmf(k, ROWS, 0.5) * BALLS, "o", ms=4, color=style.CORAL)
ax.set_xticks(range(0, ROWS + 1, 2))
ax.set(title="10,000 balls", xlabel="where the ball lands (rights)", ylabel="balls")
def plot_chances(ax: Axes, tries: int = 20) -> None:
"""Binomial distributions for the same number of tries and different chances."""
k = np.arange(tries + 1)
for p, color in ((0.1, style.CORAL), (0.5, style.BLUE), (0.8, style.GREEN)):
ax.plot(
k, stats.binom.pmf(k, tries, p), "o-", ms=3, lw=1.2, color=color, label=f"{p:.0%}"
)
ax.yaxis.set_major_formatter(style.PERCENT)
ax.legend(title="chance of yes", fontsize=8, title_fontsize=8)
ax.set(title=f"{tries} tries", xlabel="number of yeses", ylabel="probability")
fig, axes = style.subplots(3, height=3, width=10.5)
for ax, plot in zip(axes, [plot_board, plot_landing, plot_chances]):
plot(ax)
fig.tight_layout()
plt.show()
Most balls land near the middle, because most paths mix lefts and rights. It is the same story as the heights in Chapter 1.1: many small pushes, and a bell appears. With a chance far from 50% the bell is lopsided, squeezed against 0 or against the number of tries, because a count cannot go below 0 or above the number of tries.
The probability of exactly k yeses in n tries, each with chance p, is
Read it from right to left. One particular path with k yeses and n − k noes has probability : multiply the chances of each step. The binomial coefficient , “n choose k”, counts how many paths end in the same place. On the Galton board there is only one path to the far left (all lefts) but 924 paths to the middle of 12 rows; that is why the middle fills up.
On average the count is : 20 tries at 10% give 2 yeses. Its variance is , largest at p = 50%, when each try is least predictable. A Bernoulli variable is the binomial with one try: mean p, variance p(1 − p).
Many counts have no fixed number of tries: emails per hour, accidents per day, typos per page. There are endless chances for the event, and each one is tiny. Count such events in a fixed window and you get a Poisson distribution. It has one parameter: the average count.
The customers of Chapter 2.1, arriving about every two minutes, show the link. The gaps between them were exponential. Count them in ten-minute windows instead:
WINDOW, AVERAGE = 10, 5 # minutes per window, arrivals per window on average
arrivals = np.sort(rng.uniform(0, 60_000, 30_000)) # 1,000 hours, one arrival every 2 minutes
counts = np.histogram(arrivals, bins=np.arange(0, 60_001, WINDOW))[0]
def plot_windows(ax: Axes) -> None:
"""The first hour of arrivals, cut into windows, with the count in each."""
hour = arrivals[arrivals < 60]
ax.vlines(hour, 0, 1, color=style.BLUE, lw=1.2)
for start in range(0, 60, WINDOW):
if start // WINDOW % 2 == 0:
ax.axvspan(start, start + WINDOW, color=style.GREY, alpha=0.12, lw=0)
n = int(((hour >= start) & (hour < start + WINDOW)).sum())
ax.text(start + WINDOW / 2, 1.15, str(n), ha="center", fontsize=11, color=style.CORAL)
ax.set(title="Counted per 10 minutes", xlabel="minutes", xlim=(0, 60), ylim=(0, 1.4))
style.clean(ax, y=False)
def plot_counts(ax: Axes) -> None:
"""How often each count turns up, against the Poisson distribution."""
k = np.arange(16)
ax.bar(k, np.bincount(counts, minlength=16)[:16], color=style.BLUE, alpha=0.55)
ax.plot(k, stats.poisson.pmf(k, AVERAGE) * len(counts), "o", ms=4, color=style.CORAL)
ax.set(title="6,000 windows", xlabel="arrivals in 10 minutes", ylabel="windows")
def plot_averages(ax: Axes) -> None:
"""Poisson distributions: small averages are lopsided, large ones look like a bell."""
k = np.arange(31)
for mean, color in ((1, style.CORAL), (5, style.BLUE), (15, style.GREEN)):
ax.plot(k, stats.poisson.pmf(k, mean), "o-", ms=3, lw=1.2, color=color, label=str(mean))
ax.yaxis.set_major_formatter(style.PERCENT)
ax.legend(title="average", fontsize=8, title_fontsize=8)
ax.set(title="Poisson shapes", xlabel="count", ylabel="probability")
fig, axes = style.subplots(3, height=3, width=10.5)
for ax, plot in zip(axes, [plot_windows, plot_counts, plot_averages]):
plot(ax)
fig.tight_layout()
plt.show()
The Poisson probability of k events, when λ (“lambda”) are expected, is
where . It is the binomial with a huge number of tries and a tiny chance each, keeping the expected count fixed:
def plot_limit(ax: Axes) -> None:
"""Binomials with n·p = 5 for growing n, against the Poisson."""
k = np.arange(16)
for n, color in ((10, style.GOLD), (25, style.GREEN), (1000, style.BLUE)):
ax.plot(
k,
stats.binom.pmf(k, n, 5 / n),
"o-",
ms=3,
lw=1.1,
color=color,
label=f"binomial, {n} tries at {5 / n:.1%}",
)
ax.plot(
k,
stats.poisson.pmf(k, 5),
"s",
ms=7,
mfc="none",
mec=style.CORAL,
mew=1.5,
label="Poisson, average 5",
)
ax.yaxis.set_major_formatter(style.PERCENT)
ax.legend(fontsize=8)
ax.set(xlabel="count", ylabel="probability")
fig, ax = plt.subplots(figsize=(8, 3.2))
plot_limit(ax)
plt.show()
The Poisson has one fingerprint worth remembering: its variance equals its mean, λ. A count that averages 5 has a variance of about 5, a standard deviation of about 2.2.
How many tries until the first success? If each try succeeds with chance p, the first success comes on try k when k − 1 failures come first:
That is the geometric distribution, the count version of the exponential waiting time of Chapter 2.1, and like it, it has no memory: after ten failures, the next try still succeeds with chance p.
def plot_geometric(ax: Axes) -> None:
"""Geometric probabilities for three chances of success."""
k = np.arange(1, 26)
for p, color in ((0.5, style.CORAL), (0.2, style.BLUE), (0.08, style.GREEN)):
ax.plot(
k,
stats.geom.pmf(k, p),
"o-",
ms=3,
lw=1.1,
color=color,
label=f"chance {p:.0%}: on average {1 / p:.1f} tries",
)
ax.yaxis.set_major_formatter(style.PERCENT)
ax.legend(fontsize=8)
ax.set(xlabel="try on which the first success comes", ylabel="probability")
fig, ax = plt.subplots(figsize=(8, 3.1))
plot_geometric(ax)
plt.show()
Counting the failures before the r-th success instead of the first gives the negative binomial, the distribution of the next section, though there it earns its place for another reason.
The Poisson assumes every unit has the same rate: every window, every person, every page. Real units are rarely that alike. Some people see a doctor once a year, others every month. Mix many different rates and the counts spread out more than a Poisson allows: more zeros, and a longer tail. This is called overdispersion, and the negative binomial distribution is made for it.
per_person = visits.mdvis.astype(int)
MOST = 20
k = np.arange(MOST + 1)
people = np.bincount(per_person, minlength=MOST + 1)[: MOST + 1]
mean, variance = per_person.mean(), per_person.var()
def plot_fit(ax: Axes, probabilities: np.ndarray, title: str) -> None:
"""The counts of people per number of visits, with one distribution over them."""
ax.bar(k, people, color=style.GOLD, alpha=0.7)
ax.plot(k, probabilities * len(per_person), "o-", ms=3, lw=1.2, color=style.CORAL)
ax.set(title=title, xlabel="doctor visits per year", ylabel="people")
fitted = sm.NegativeBinomial(per_person, np.ones(len(per_person))).fit(disp=0)
n = 1 / fitted.params.iloc[-1] # the fitted shape: small n means rates differ a lot
p = n / (n + mean)
fig, axes = style.subplots(2, height=3.1, width=9)
plot_fit(axes[0], stats.poisson.pmf(k, mean), "Poisson: too tidy")
plot_fit(axes[1], stats.nbinom.pmf(k, n, p), "Negative binomial")
style.note(axes[0], f"mean {mean:.1f}\nvariance {variance:.1f}", x=0.6, color=style.INK)
fig.tight_layout()
plt.show()
The mechanism is easy to simulate. Give each person their own yearly rate, drawn from a gamma distribution (Chapter 2.1), then draw each person’s visits from a Poisson with that rate:
rates = rng.gamma(n, mean / n, len(per_person)) # one rate per simulated person
simulated = rng.poisson(rates)
def plot_rates(ax: Axes) -> None:
"""The spread of personal rates."""
style.hist(ax, rates[rates < 20], bins=50, color=style.PURPLE)
ax.axvline(mean, color=style.INK, ls="--", lw=1.2)
ax.set(title="Each person's own rate", xlabel="expected visits per year", ylabel="people")
def plot_simulated(ax: Axes) -> None:
"""Counts from the mixed rates against the negative binomial and a single Poisson."""
ax.bar(
k, np.bincount(simulated, minlength=MOST + 1)[: MOST + 1], color=style.PURPLE, alpha=0.5
)
ax.plot(
k, stats.nbinom.pmf(k, n, p) * len(simulated), "o-", ms=3, lw=1.2, color=style.CORAL
)
ax.plot(k, stats.poisson.pmf(k, mean) * len(simulated), "o-", ms=3, lw=1, color=style.GREY)
ax.set(title="Their visits", xlabel="visits per year", ylabel="people")
fig, axes = style.subplots(2, height=3.1, width=9)
plot_rates(axes[0])
plot_simulated(axes[1])
fig.tight_layout()
plt.show()
A Poisson mixed over gamma-distributed rates is the negative binomial. Its variance is , where (the shape) measures how alike the units are: as grows, the rates become equal and the negative binomial turns into the Poisson.
The check is quick: compare the variance with the mean. For doctor visits the average is under 3 visits but the variance is several times larger, far from the “variance = mean” of a Poisson. A model that assumes Poisson here would be far too sure of itself; Chapter 4.6 returns to this.
The binomial counts two outcomes. With more categories, such as the faces of a die or the parties in an election, the counts in all categories together follow a multinomial distribution. Each single category on its own is still binomial (that category against all others), and the counts are tied: if one category gets more, the others must get fewer, so their covariances are negative.
rolls = rng.multinomial(60, [1 / 6] * 6, size=500)
def plot_faces(ax: Axes) -> None:
"""Counts of each face across experiments."""
for face in range(6):
ax.scatter(
face + 1 + rng.uniform(-0.2, 0.2, len(rolls)),
rolls[:, face],
s=4,
color=style.BLUE,
alpha=0.3,
)
ax.axhline(10, color=style.CORAL, lw=1.4)
ax.set(
title="Six counts", xlabel="face", ylabel="times rolled in 60 rolls", xticks=range(1, 7)
)
def plot_pair(ax: Axes) -> None:
"""Ones against twos."""
jitter = rng.uniform(-0.25, 0.25, (len(rolls), 2))
ax.scatter(
rolls[:, 0] + jitter[:, 0], rolls[:, 1] + jitter[:, 1], s=6, color=style.BLUE, alpha=0.4
)
r = np.corrcoef(rolls[:, 0], rolls[:, 1])[0, 1]
style.note(ax, f"correlation {r:.2f}", x=0.03, y=0.98, color=style.INK)
ax.set(title="Tied together", xlabel="ones", ylabel="twos")
fig, axes = style.subplots(2, height=3.1, width=9)
plot_faces(axes[0])
plot_pair(axes[1])
fig.tight_layout()
plt.show()
The correlation between any two of k equally likely categories is −1/(k − 1): −0.2 for a die. Chapter 3.5 tests whether observed category counts fit expected ones.
| Mechanism | Distribution | Mean | Variance | |
|---|---|---|---|---|
| one try, yes or no | Bernoulli (p) | p or 1 − p | p | p(1 − p) |
| a fixed number of tries, count the yeses | binomial (n, p) | np | np(1 − p) | |
| many chances, each rare, same rate | Poisson (λ) | λ | λ | |
| tries until the first success | geometric (p) | 1/p | (1 − p)/p² | |
| rare events, rates differ between units | negative binomial (μ, r) | μ | μ + μ²/r | |
| counts in several categories | multinomial (n, p₁ … p_k) | np_j | np_j(1 − p_j) |
Next: how to find the distribution that fits, and how to check it.